COALIA: A Computational Model of Human EEG for
Consciousness Research

Siouar Bensaid, Julien Modolo, Isabelle Merlet, Fabrice Wendling, Pascal Benquet

► To cite this version:

Siouar Bensaid, Julien Modolo, Isabelle Merlet, Fabrice Wendling, Pascal Benquet. COALIA: A Computational Model of Human EEG for Consciousness Research. Frontiers in Systems Neuroscience, 2019, 13, pp.59. 10.3389/fnsys.2019.00059. hal-02430077

HAL Id: hal-02430077

https://univ-rennes.hal.science/hal-02430077

Submitted on 7 Jan 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for
updates

# COALIA: A Computational Model of Human EEG for Consciousness Research

Siouar Bensaid†, Julien Modolo†, Isabelle Merlet, Fabrice Wendling* and Pascal Benquet

INSERM, Laboratoire Traitement du Signal et de l'Image (LTSI)—U1099, University of Rennes, Rennes, France

Understanding the origin of the main physiological processes involved in consciousness is a major challenge of contemporary neuroscience, with crucial implications for the study of Disorders of Consciousness (DOC). The difficulties in achieving this task include the considerable quantity of experimental data in this field, along with the non-intuitive, nonlinear nature of neuronal dynamics. One possibility of integrating the main results from the experimental literature into a cohesive framework, while accounting for nonlinear brain dynamics, is the use of physiologically-inspired computational models. In this study, we present a physiologically-grounded computational model, attempting to account for the main micro-circuits identified in the human cortex, while including the specificities of each neuronal type. More specifically, the model accounts for thalamo-cortical (vertical) regulation of cortico-cortical (horizontal) connectivity, which is a central mechanism for brain information integration and processing. The distinct neuronal assemblies communicate through feedforward and feedback excitatory and inhibitory synaptic connections implemented in a template brain accounting for long-range connectome. The EEG generated by this physiologically-based simulated brain is validated through comparison with brain rhythms recorded in humans in two states of consciousness (wakefulness, sleep). Using the model, it is possible to reproduce the local disynaptic disinhibition of basket cells (fast GABAergic inhibition) and glutamatergic pyramidal neurons through long-range activation of vasoactive intestinal-peptide (VIP) interneurons that induced inhibition of somatostatin positive (SST) interneurons. The model (COALIA) predicts that the strength and dynamics of the thalamic output on the cortex control the local and long-range cortical processing of information. Furthermore, the model reproduces and explains clinical results regarding the complexity of transcranial magnetic stimulation TMS-evoked EEG responses in DOC patients and healthy volunteers, through a modulation of thalamo-cortical connectivity that governs the level of cortico-cortical communication. This new model provides a quantitative framework to accelerate the study of the physiological mechanisms involved in the emergence, maintenance and disruption (sleep, anesthesia, DOC) of consciousness.

**Keywords:** computational modeling, brain connectivity, feedforward inhibition, GABA, disinhibition, TMS-EEG, disorders of consciousness (DOC)

**OPEN ACCESS**

**Edited by:**

Igor Timofeev,
Laval University, Canada

**Reviewed by:**

Marco Atzori,
Universidad Autónoma de San Luis
Potosà, Mexico

Jamie Sleigh,
The University of Auckland,
New Zealand

*Correspondence:

Fabrice Wendling
fabrice.wendling@inserm.fr

†Co-first authors

**Received:** 08 March 2019

**Accepted:** 07 October 2019

**Published:** 13 November 2019

**Citation:**

Bensaid S, Modolo J, Merlet I, Wendling F and Benquet P (2019) COALIA: A Computational Model of Human EEG for Consciousness Research. *Front. Syst. Neurosci.* 13:59. doi: 10.3389/fnsys.2019.00059
# INTRODUCTION

The characterization and understanding of the mechanisms underlying consciousness is one, if not the greatest challenge that contemporary neuroscience is facing. Beyond the purely fundamental interest of this question, a major clinical issue is also at stake: evaluating residual consciousness in patients suffering from Disorders of Consciousness (DOC), which can be extremely difficult, and have crucial implications in terms of clinical care. For example, motor imagery paradigms can reveal covert consciousness in coma patients, using functional magnetic resonance imaging (fMRI; Owen et al., 2006) for instance. This illustrates the pressing need for an improved characterization of the mechanisms that underlie consciousness, which could be exploited to propose novel quantified measures, or *metrics*, of consciousness.

Many theories have attempted, at various levels of description, to integrate the multifaceted aspects of consciousness. One of the first theories that found a significant echo in the neuroscience community is the Dynamic Core Hypothesis (DCH; Tononi and Edelman, 1998), which was the first to relate the concept of information with consciousness. In this theory, functional clusters in the thalamocortical system are central, and involve fast re-entrant interactions, as well as a high level of integration and differentiation giving rise to complex patterns of neuronal activity. Another popular theory that has been gradually expanded over the years, and that has solid ties with neurophysiology, is the Global Workspace Theory (GWT; Dehaene et al., 1998; Dehaene and Changeux, 2011). In short, GWT states that conscious information is globally available within the brain, and that the “ignition” of large-scale networks, i.e., the sudden communication between distant brain regions to engage into the processing of information, enables a stimulus to reach the global workspace, hence consciousness. Ignition is thought to involve long-range glutamatergic fibers that enable long-distance communication between cortical regions. Several experiments have supported GWT, for example that non-masked words involve the activation of much wider networks as compared to masked words (Dehaene et al., 2001), with a similar result found for sub-liminal vs. supra-threshold visual stimuli (Modolo et al., 2018; van Vugt et al., 2018). The Integrated Information Theory (IIT) is based on a different approach. Instead of beginning from the large-scale structure of the thalamo-cortical system as in the DCH and GWT, IIT introduces several axioms to derive general principles of consciousness. One of the leading ideas of IIT is that consciousness involves the *integration* of information between distant areas (reminiscent of ignition in GWT), which increases the complexity of the processed information. Segregation, or *differentiation* of information, is also key in IIT: for example, large-scale synchronization of several regions with the same activity is indeed integrated, however with low complexity (Koch et al., 2016). Therefore, *integration* and *differentiation* appear as the two concepts leading to increasing the *complexity* of the information conveyed by brain-scale networks. A more recent theory named algorithmic information theory of consciousness a.k.a.

Kolmogorov Theory (KT; Ruffini, 2017) is also based on the idea that conscious states are associated with higher levels of complexity, and that subjective experience occurs following processes of information compression.

In contrast with the aforementioned theoretical studies of consciousness, only few studies have actually attempted to simulate brain activity associated with consciousness states using neurophysiologically-plausible computational models. Obviously, the *in silico* implementation of neural mechanisms that underlie the emergence and maintenance of consciousness represents a considerable challenge. However, capturing the main features of the most significant common principles from the main theories of consciousness using a computational neuroscience framework appears at reach. For example, a computational model exploring how cortico-cortical connectivity is functionally impaired during sleep (“connectivity breakdown”) has been proposed (Esser et al., 2009). However, most of these models are limited in terms of spatial scale and the represented micro-circuitry. This limitation hinders bridging the micro-circuit scale with the brain-scale, which is of interest in consciousness. The present study proposes to fill this gap, and provides new links between different levels of description (from local neuronal population to whole-brain scale).

Using a bottom-up approach, we developed COALIA, a new computational model of brain-scale electrophysiological activity. This is neither the intention nor objective of the present article to investigate the philosophical question of “qualia” that represents the perception properties. Instead, the model name “COALIA” has been chosen to ground this brain-scale computational model in the field of consciousness research. The model starts from neuronal micro-circuits involving different cellular subtypes that have been reliably identified through neurobiological studies. The basic unit of the model is the neural mass, representing a local population of a few thousands of neurons, which has proven its ability to capture the dynamics of actual neuronal assemblies (Wendling et al., 2002). At the local level, the model includes subsets of pyramidal neurons (glutamatergic), and three different types of interneurons (GABAergic) with appropriate physiologically-based kinetics (fast vs. slow). At the global level, the large-scale model is then constructed on the basis of a standard 66-region brain atlas (Desikan et al., 2006), with one neural mass representing the local field activity of one atlas region. Neural masses are spatially distributed over the cortex, using the template brain morphology (Colin). As they account for distinct cortical regions, neural masses are synaptically connected through long-range glutamatergic projections among pyramidal neurons and GABAergic interneurons, Connectivity is derived from DTI (Diffusion Tensor Imaging) data. Results show that the model captures the large-scale structure of brain connectivity between regions while accounting for the properties of local micro-circuits. It can accurately reproduce EEG activity for different conscious states (e.g., sleep vs. wake), and the breakdown of functional connectivity during sleep as assessed through the replication of TMS-EEG experiments.

In this article, we first describe basic concepts of cortico-cortical and thalamo-cortical networks involved in theories of
consciousness. Then, we present the neural mass modeling approach used to develop the brain-scale model, along with the various microcircuits considered and their functional role. A toy model involving a limited number of neuronal populations is then investigated in order to validate the implemented micro-circuits, before an extension to the whole-brain model. Simulated responses to TMS are generated and quantified in two consciousness states, namely awake and asleep. Results are discussed according to the novelty, performance and limitations of the model, along with its usefulness in consciousness studies. Future extensions are described.

# BACKGROUND: ROLE OF CORTICO- AND THALAMO-CORTICAL NETWORKS IN CONSCIOUSNESS

Consciousness is a global functional state of the brain that is intrinsically linked with neuronal oscillations generated by large-scale cortico-cortical and thalamo-cortical networks (Llinás et al., 1998). More specifically, wakefulness is determined by widespread thalamocortical projections (Timofeev and Steriade, 1996; Laureys, 2004) while awareness requires the activation of a wide cortico-cortical network, involving lateral and medial frontal regions parieto-temporal and posterior parietal areas, bilaterally (Laureys et al., 1999). In the model, we included these two key components of consciousness that are briefly reviewed below.

## "Horizontal" Cortico-Cortical Connectivity

Functional connectivity studies have shed light on the functional networks involved in various conscious states (Jin and Chung, 2012). During general anesthesia-induced loss of consciousness, there is a breakdown in cortical effective connectivity (Ferrarelli et al., 2010; Hudetz, 2012; Gómez et al., 2013). As a reminder, effective connectivity is defined as the ability of a neuronal group to causally affect the firing of other neuronal groups (Friston, 2011). In unresponsive patients, impaired consciousness was associated with altered effective connectivity (Varotto et al., 2014; Crone et al., 2015). A protocol of TMS triggered a simple local EEG response indicating a breakdown of effective connectivity at the cortical level, similar to the one previously observed in unconscious sleeping or anesthetized subjects (Casali et al., 2013). Sleep stages have a drastic impact on consciousness and also on functional connectivity. For example, there is a strong reduction of both wakefulness and awareness components of consciousness during NREM sleep, associated with thalamic up-and-down state and cortical slow wave sleep (SWS). However, brain effective connectivity changes significantly (Tononi and Sporns, 2003; Esser et al., 2009). More specifically, cortical activations become more local and stereotypical upon falling into NREM sleep, which indicates an impaired effective cortical connectivity (Massimini et al., 2010). "Horizontal" communication through coherence that involves high frequency oscillations synchronization is a fundamental mechanism in cortical function and perception (Fries, 2005, 2009). Overall, the "awareness" component of consciousness depends on large-scale synchronized communication among

distant neuronal populations distributed over the neocortex (see Figure 1A, left panel).

## "Vertical" Thalamo-Cortical Connectivity

As reported in IIT (Tononi, 2004, 2012), consciousness depends on the brain's ability to integrate information, which relies on the effective connectivity among functionally specialized regions (or clusters) of the thalamocortical system, and on the segregation of information. One important modulator of cortical connectivity is the activity pattern of thalamocortical cells, tonic vs. up- and down, which is able to modify the excitability level of cortical neuron subpopulations. The thalamic-mediated synchronization of distant cortical areas may coordinate the large-scale integration of information across multiple cortical circuits, consequently influencing the level of arousal and consciousness (Saalmann, 2014). Conversely, during sleep or anesthesia-induced transitions in consciousness, both thalamo-cortical and intra-thalamic functional connectivity are modified (Kim et al., 2012; Hale et al., 2016). In addition, thalamic input to neocortex modifies cortico-cortical connectivity. Upon falling into NREM sleep (when rhythmic thalamo-cortical up-and-down activity develops), cortical activations become more local and stereotypical, indicating a significant decrease of cortico-cortical connectivity (Esser et al., 2009; Massimini et al., 2010; Usami et al., 2015). Recently, it was shown that direct and tonic optogenetic activation of thalamic reticular nuclei (TRN) GABAergic interneurons induces a spatially restricted cortical slow-wave activity (Lewis et al., 2015). This activity was reminiscent of sleep rhythms, and animals exhibited behavioral changes that were consistent with a decrease of arousal.

Overall, this brief literature review suggests that both components of consciousness, namely awareness and wakefulness, are impaired when large-scale cortico-cortical functional connectivity mediated through the binding of synchronized high-frequency oscillations in the beta-gamma band (Schoenberg et al., 2018) and regulated through the thalamus (Nakajima and Halassa, 2017). Meanwhile, during this decrease of awareness and wakefulness, an increase of "vertical" thalamo-cortical connectivity is observed, along with a stronger synchronization of delta oscillations between thalamic cell (TC) assemblies and isolated groups of neocortical neurons (Hill and Tononi, 2005; Figure 1A, right panel).

# MATERIALS AND METHODS

## Modeling of Micro- and Macro-Circuits: Neural Mass Model Approach

Neural mass models (NMMs) are a mathematical description of neural dynamics at a mesoscopic scale (from a millimeter to several centimeters of the cortex). This class of models was proposed in the 1970s as an alternative to detailed microscopic models that require a more extensive computational cost (Wilson and Cowan, 1973; Nunez, 1974; Lopes da Silva et al., 1976; Freeman, 1978). NMMs can indeed model the local field potential (LFP) of an entire cortical region using only few state variables (Breakspear, 2017), whereas in detailed
A

«Horizontal» cortico-cortical connectivity

<figure><img src="image_3.png" /></figure>

«Vertical» thalamo-cortical connectivity

<figure><img src="image_4.png" /></figure>

B

<figure><img src="image_5.png" /></figure>

**FIGURE 1** | (A) Illustration of “horizontal” and “vertical” connectivity. Left: Horizontal connectivity refers to the cortico-cortical connections which are functionally effective during wakefulness, with a weak level of thalamo-cortical coupling. Right: Vertical connectivity refers to thalamo-cortical projections that functionally impair cortico-cortical connectivity during sleep. (B) General architecture of the micro- and macro-circuits implemented in the computational model. The local neural mass model (NMM) of cortical activity is composed of a PC exciting two GABAergic interneurons, namely, the somatic-projecting BC and the dendritic-projecting SST, responsible for the generation of fast and slow oscillations, respectively; the VIP were introduced as they play a crucial role in cortical column communication through disinhibition of SST. The subcortical module consisted in TC sending excitatory glutamatergic projections to the TRN block composed of fast and slow GABAergic interneurons TRN₁ and TRN₂, respectively. VIP, vasoactive intestinal peptide positive GABAergic interneurons; SST, Somatostatin-positive GABAergic interneurons; BC, Basket-type GABAergic interneurons; PC, Glutamatergic Pyramidal Cells.

models this activity is meticulously described at the level of spatially distributed and interconnected neuron models, each including the properties of ionic channels, axons and dendrites (Wang and Buzsáki, 1996; Whittington et al., 2000; Maex and De Schutter, 2007).

Despite their simplicity, NMMs are neurophysiologically grounded, since they include the connectivity, synaptic kinetics

and firing rates of neuronal sub-types present in the region of interest. The reduced complexity and performance (in term of reproducing actual LFPs) of the NMM approach made it a powerful tool to investigate various cerebral mechanisms, such as the generation of brain rhythms (Jansen and Rit, 1995; David and Friston, 2003; Ursino et al., 2010). NMMs have also been extensively used to study pathological dynamics
such as in epilepsy (Wendling et al., 2002; Traub et al., 2005; Molaee-Ardekani et al., 2010; for a review, see Wendling et al., 2016), Alzheimer's disease (AD; Bhattacharya et al., 2011) and Parkinson's disease (PD; Liu et al., 2016, 2017).

Designing a NMM involves identifying the main subsets of neurons implied in the modeled brain tissue and describing their synaptic interconnections. Based on an exhaustive literature review, we developed a model consisting of coupled NMMs able to simulate both cortical and thalamic activity, as described hereafter.

## A Local Neural Mass Model of Neocortical Activity

Until recently, the classification of GABAergic interneurons was a highly challenging task, regarding the electrophysiological properties, morphology, biochemistry markers and connectivity. However, all recent studies converge towards much simpler functional categories, involving three main classes of GABAergic interneurons in the neocortex (Rudy et al., 2011). As briefly reviewed below, these classes account for: (1) somatic-targeting parvalbumine positive (PV+) basket cells (BC); (2) dendritic-targeting somatostatin positive (SST) interneurons; and (3) vasoactive intestinal-peptide (VIP) expressing interneurons (Tremblay et al., 2016; **Figure 1B**, left panel).

### Basket Cells and Fast Oscillations

PV+ participate to the generation of cortical gamma oscillations through: (1) thalamocortical feedforward inhibition in layer 4; (2) feedback inhibition in layer 2/3; and (3) *via* direct PV+/PV+ coupling through electrical gap-junctions (Povysheva et al., 2008; Buzsáki and Wang, 2012; Lewis et al., 2012; Varga et al., 2014; Womelsdorf et al., 2014; Chen et al., 2017). Recent optogenetic experiments have demonstrated the causal role of somatic-targeting interneurons BC in mediating fast oscillations (>20 Hz; Chen et al., 2017).

### Surround Inhibition by SST GABAergic Interneurons and Slow Oscillations

Neocortical SST neurons can exhibit high levels of spontaneous slow oscillations, and their tonic activity might facilitate fine scale up- and down-regulation of global inhibition levels in the neocortex (Urban-Ciecko and Barth, 2016). Lateral inhibition has been shown to be a fundamental principle in neuronal networks (Harris and Mrsic-Flogel, 2013; Karnani et al., 2014; Harris and Gordon, 2015). Specifically, lateral inhibition between nearby pyramidal cells (PCs) is thought to work through SST interneurons (Kapfer et al., 2007; Silberberg and Markram, 2007; Adesnik et al., 2012). Dendritic inhibition is more effective than perisomatic inhibition in regulating excitatory synaptic integration, therefore SSTs are the key regulator of input-output transformations (Lovett-Barron et al., 2012). Furthermore, the slow kinetics of dendritic-targeting inhibitory postsynaptic membrane potentials (IPSPs) are particularly suited to maximize the localized shunting inhibition effect (Gidon and Segev, 2012; Paulus and Rothwell, 2016). Finally, it was shown that this cell type is especially involved in slow oscillations (Womelsdorf et al., 2014; Urban-Ciecko and Barth, 2016; Funk et al., 2017).

## Cortical Column Communication Through VIP-Controlled Disinhibition

The disinhibition of cortical PCs gates information flow through and between cortical columns (Walker et al., 2016). One of the factors underlying this PC disinhibition is the inhibition of SST cells during active cortical processing, which enhances distal dendritic excitability (Gentet et al., 2012). The activation of VIP neurons strongly inhibits dendritic-targeting SST interneurons mediating surround inhibition, which leads to PCs disinhibition (suppress the inhibition on PCs; Lee et al., 2013; Pi et al., 2013; Fu et al., 2014; Pfeffer, 2014; Yang et al., 2016). This disynaptic mechanism of disinhibition seems to be a generic motif able to suppress the blanket of inhibition mediated by SST neurons (Fino and Yuste, 2011; Karnani et al., 2014). It has been demonstrated indeed in motor, sensory and associative neocortical areas that transient SST activity suppression by VIP activation occurs during visual processing (Lee and Mihalas, 2017), somatosensory integration (Lee et al., 2013; Sohn et al., 2016), locomotion (Dipoppa et al., 2018), top-down modulation (Ayzenshat et al., 2016) and plasticity during perceptual learning (Williams and Holtmaat, 2019). Since VIP neurons are targeted by long-range cortical glutamatergic projections, they represent a key factor for distal cortico-cortical activation through disynaptic disinhibition.

### Formal Description of Neocortical Model

Based on the above information, we have developed a neocortical module involving PCs and three types of inhibitory subpopulations, namely, BC, SST and VIP (see **Figure 1B**, left panel). BC and SST receive excitatory inputs from PCs that are reciprocally inhibited by both of them. Pyramidal collateral excitation was also implemented *via* an excitatory feedback loop passed by a supplementary excitatory population (PC) analogous to PC, except that it projects only from and to subpopulation PC. The electrical gap-junction mentioned in “Basket Cells and Fast Oscillations” section was implemented through an inhibitory feedback loop characterized by a connectivity constant $C_{BC}^n$, where *n* is the index of the NMM (e.g., *n* = 23th neural mass). Communication through disinhibition mentioned in “Cortical Column Communication Through VIP-Controlled Disinhibition” section was modeled by inhibitory projections; first from VIP to SST, and second, from the latter to BC. The nonspecific influence from neighboring and distant populations was modeled by a Gaussian input noise corresponding to an excitatory input $p_C^n(t)$ that globally describes the average density of afferent action potentials. The set of ordinary differential equations (ODEs) modeling the neocortical module is included as a **Supplementary Material** (see Supplementary section “ODEs of the Neocortical Module”).

## A Local Neural Mass Model of Thalamic Activity

The thalamus is considered as a complex relay extensively connected with the cortex, as well as most subcortical areas. This central position underlies its key role in several cognitive functions including perception, attention, memory and consciousness. Importantly, even a limited damage in
the thalamus can have major consequences on all the aforementioned functions (Ward, 2013). The thalamo-cortical circuitry as a neural correlate of consciousness has been mentioned by consciousness theories. More evidence for this role was provided by a study analyzing the metabolic activity of posterior midline cortical areas driven by the thalamic nuclei, across different altered consciousness states (Laureys et al., 2006). The authors reported that unresponsive wakefulness syndrome (UWS) patients can be differentiated from minimally conscious state (MCS) patients by a difference in glucose metabolism in these areas. The pivotal role of thalamo-cortical loop in the generation of slow waves and the up-and-down state was appraised in the review of Crunelli et al. (2015) where they enumerated the main lines of evidence supporting this assertion, namely, the strong interconnection between thalamic cells (TCs) and neocortical layers involved in slow waves suggesting that thalamic nuclei can control up-and-down state dynamics in neocortical circuits, the early TC firing in relation to the initiation of cortical UP states, the rhythmic up-and-down state generated by TC neurons and TRN in isolated conditions, and the neocortical UP states readily induced in head-restrained mice by selective optogenetic activation of TC neurons proving the stimuli role of the latter.

The importance of thalamo-cortical connectivity has motivated the development of thalamo-cortical models simulating the interactions between the cortex and thalamus at a mesoscopic level (Suffczynski et al., 2004; Sotero et al., 2007; Bhattacharya et al., 2011; Roberts and Robinson, 2012; Mina et al., 2013; Sen Bhattacharya et al., 2013; Cona et al., 2014). While some models were developed to generate alpha activity (8–12 Hz; Sotero et al., 2007; Bhattacharya et al., 2011; Sen Bhattacharya et al., 2013), others were used to simulate the sleep-wake cycle (Suffczynski et al., 2004; Roberts and Robinson, 2012; Cona et al., 2014).

Since our aim is to develop a computational model able to reproduce brain rhythms corresponding to different consciousness states (e.g., sleep-wake cycle), the inclusion of the thalamus in the model is crucial. A description of the thalamus model is provided in the right panel of **Figure 1B**. The thalamic module includes one population of excitatory glutamatergic neurons TCs, and two inhibitory interneurons from the TRN, TRN₁ and TRN₂ accounting for fast and slow GABAergic IPSPs, respectively. TCs receive GABAergic IPSPs with slow and fast kinetics from the TRNs, whereas the latter receive excitatory inputs from the former. Similarly to the cortical module, a Gaussian input noise corresponding to excitatory input $p_{Th}^n(t)$ was used to represent nonspecific inputs on TCs. The set of ODEs modeling the thalamic module is provided as a **Supplementary Material** (see Supplementary section “ODEs of the Thalamic Module”).

# Modeling of Large-Scale Cortico-Cortical and Thalamo-Cortical Connectivity

## Cortico-Cortical Connections

PCs originating from a single cortical column target several cell types in distant cortical columns. Glutamatergic PCs target

not only remote PCs by common feedforward excitation but also GABAergic cells by disynaptic cortico-cortical feedforward inhibition (FFI; see **Figure 2**). Even if $PV^+$ BCs have been shown to receive stronger thalamocortical and intracortical excitatory inputs than SST neurons, it appears that cortico-cortical FFI could be mediated by both types of interneurons (Ma et al., 2010; Tremblay et al., 2016). Importantly, cortico-cortical glutamatergic axons also target VIP neurons (Sohn et al., 2016).

Therefore, feedforward excitation was included in the model by means of a connectivity constant $K_{P,P}^{ij}$ modeling the average strength of glutamatergic projections from pyramidal subpopulations in NMM “i” to their counterpart in NMM “j.” Disynaptic cortico-cortical feedforward inhibition was similarly integrated via the connectivity constants $K_{P,BC}^{ij}$, $K_{P,SST}^{ij}$ and $K_{P,VIP}^{ij}$ denoting glutamatergic projections from PC subpopulations in NMM “i” to BC, SST and VIP in NMM “j,” respectively. The set of ODEs modeling the cortico-cortical connections is included as a **Supplementary Material** [see Supplementary section “ODEs of Large Scale Cortico- and Thalamo-Cortical Connectivity,” equations (9)–(12)].

In long-range cortico-cortical connections, the time-delay between NMMs “i” and “j” was controlled by a distance parameter $D^{ij}$ that reports the Cartesian distance in centimeters (cm) between the two populations. By setting the conduction velocity of action potentials in the brain to $c$ [$c \in (10,100)$ cm/s], the time-delay taken by NMM “j” to receive a firing rate is straightforwardly deduced as $D^{ij}/c$.

## Thalamo-Cortical Connections

The main connections between the thalamus and neocortex were taken into consideration in the model (see **Figure 2**). As in classical thalamocortical models, TCs receive glutamatergic excitatory postsynaptic potentials (EPSPs) from PCs. In turn, PCs receive excitatory input from TCs. Similarly, TRNs also receive excitatory cortical projections. In terms of GABAergic cortical targets, thalamic projections mainly target $PV^+$ basket cells (Cruikshank et al., 2007; Yang et al., 2013). In the adult brain, thalamic projections onto SST neurons are present but are much weaker as compared to the projections onto PCs and $PV^+$ cells (Ji et al., 2016). However, robust thalamocortical activation of non-Martinotti dendritic-targeting GABAergic interneurons has been demonstrated (Tan et al., 2008). Long-range connections from cortical areas and/or thalamic nucleus areas can activate VIP neurons, which in turn inhibit SST neurons, and disinhibit PCs dendrites. Such dendritic disinhibitory circuit has been proposed to gate excitatory inputs targeting pyramidal dendrites (Yang et al., 2016; Williams and Holtmaat, 2019). Note that the thalamic compartment implements FFI, since it induces first an EPSP (PCs activation) followed later on by an IPSP (cortical interneurons activation).

Since only one thalamic population was considered in this model, the TC projections to—and from—PCs were included in the model. Connectivity constants $K_{Th,P}^n$ and $K_{P,Th}^n$ allow to adjust the strength of efferent TC projections to PCs of NMM “n” and efferent pyramidal projections from the NMM “n” to TC, respectively. Projections from PCs of NMM “n” onto TRN₁ and TRN₂ were also integrated and adjusted with
<figure><img src="image_6.png" /><figcaption>**FIGURE 2** | Large-scale architecture of the model. Illustration of the synaptic projections between cortical modules, and between thalamic and cortical modules. Note the presence of long-range thalamocortical and cortico-cortical feedforward inhibition. For the sake of clarity, long-range connections between cortical NMMs and NMMs are unidirectional, whereas in the model, PCs of NMMs also project on neurons of NMMs. The strength of the synaptic input of PC onto distant VIP cells is larger than the input onto the other distant interneurons (SST, Basket).</figcaption></figure>

the connectivity constants $K_{P,TRN1}^n$ and $K_{P,TRN2}^n$, respectively. Likewise, projections from TC to GABAergic interneurons in NMM "n," namely BCs, SSTs and VIPs, were included through the connectivity constants $K_{Th,BC}^n$, $K_{Th,SST}^n$ and $K_{Th,VIP}^n$, respectively. The set of ODEs modeling the thalamo-cortical connections is described in Supplementary section "ODEs of Large Scale Cortico- and Thalamo-Cortical Connectivity" [equations (13)–(19)]. It is noteworthy that time-delays between the thalamus and the cortical NMMs were included as in the cortico-cortical long-range connections.

## Implementation and Parameter Tuning

An important step in NMM approaches consists in tuning the model parameters. Some of these parameters, namely time constants in the three modules, were set close to the "standard values" used in neuronal population models, while other parameters such as connectivity constants were adjusted according to the target EEG activities. Two classical well-known examples of conscious and unconscious states are deep sleep- a.k.a. SWS characterized by delta oscillations (0–4 Hz), and wakefulness (background activity). The objective was to reproduce these manifestations of consciousness modulation in the model.

In **Supplementary Table S1**, we provide physiological interpretation and values of model parameters. Regarding the model output, the signal simulated at the level of PCs in the cortical compartment was chosen as the model output since it corresponds to the sum of PSPs, which is the main contribution to LFPs recorded in the neocortex.

When several cortical NMMs are interconnected, a simple way to handle large-scale connectivity is to arrange the connectivity constants in arrays where line and column indices refer to source and target NMMs, respectively. Based on "Modeling of Large-Scale Cortico-Cortical and Thalamo-Cortical Connectivity" section, there are two categories of connectivity array, "excitatory to excitatory" and "excitatory to inhibitory" arrays. In the first category, we consider the matrix $K_{EXC}$ whose (i, j)th element represents the glutamatergic projections from NMM "i" onto an excitatory subpopulation in NMM "j." Hence, if both of them belongs to the cortical module, the (i, j)th element would correspond to $K_{P,P}^{ij}$. However, if NMM "i" coincides with a thalamic population, the (i, j)th element would be then equal to $K_{Th,P}^j$. In the second category, a similar scheme is considered with the difference that the target subpopulation is always an inhibitory interneuron (BC, SST, VIP or TRNs). Consequently, five arrays are considered, namely, $K_{BC}$,
$K_{SST}$, $K_{VIP}$, $K_{TRN1}$ and $K_{TRN2}$ (see Supplementary Figure S1). Note that when only one thalamic NMM is considered, the two last matrices are reduced to vector arrays. Cortico- and thalamo-cortical time-delays were similarly arranged in matrix $D$ whose $(i, j)$th element coincides with the aforementioned $D^{ij}$.

The set of second order stochastic nonlinear ODEs obtained for all synaptic interactions present in the model was numerically solved using a fixed step ($\Delta t = 1$ ms) 4th-order Runge-Kutta method. The model was implemented using an object-oriented language (Objective-C).

To study the model behavior, we followed a two-step approach. We first implemented a “toy model” with a small number of coupled neural masses. This reduced-complexity model allowed to assess the effects of cortico- and thalamo-cortical connectivity matrices on cortical activities associated with different consciousness states (see “A Local Neural Mass Model of Neocortical Activity,” “A Local Neural Mass Model of Thalamic Activity” and “Modeling of Large-Scale Cortico-Cortical and Thalamo-Cortical Connectivity” sections). After validation in the toy model, all mechanisms were implemented in an extended, more realistic model, where the whole brain was considered.

## Toy Model of Cortical Activity

The toy model was composed of one thalamic population connected to four cortical populations ($N = 4$) that were identically tuned. The default values of NMMs intrinsic parameters are listed in Supplementary Table S1. Unless explicitly mentioned, these parameters were kept unchanged for all neural activities generated afterward (SWS and background activity). The objective was to verify, in a simplified model, the hypothesis that cortical activity is modulated from deep sleep to wakefulness by tuning only thalamo- and cortico-cortical connectivity, so that, when thalamocortical connectivity increases, the model goes deeper into sleep. This connectivity mechanism was implemented in all connectivity matrices in the model (see Supplementary Figure S2).

## Whole Brain Model: From Cortical Activity to EEG

The pipeline used for the simulation of scalp EEG data is described in Figure 3. In order to obtain a “realistic” activity during wakefulness and SWS over the entire neocortex, we considered one thalamic population and 66 cortical populations ($N + 1 = 67$). Each time-course at the output of these populations represented the activity of one macro-region of the anatomical parcellation described in Desikan et al. (2006), in which the activity was assumed to be homogeneous. In order to generate 67 time-courses from the coupled NMMs (66 cortical regions plus the thalamus), we used a combination of connectivity arrays (Figure 3B). First, a matrix of connection weight $K_{DTI}$ representing a density of fibers between all pairs of 66 cortical regions was used to set the structural connectivity. This matrix, obtained from DTI, is provided in Hagmann et al. (2008). Second, we considered additional functional horizontal (i.e., cortico-cortical) matrices $K_{Hx}$, reproducing the coefficients weights used for both wakefulness and sleep in the toy model, taking also into account the new number of populations. In order to apply the connectivity weights defined earlier in the toy model

only to pairs of NMMs that are structurally connected. Structural and horizontal functional matrices were combined using the Hadamard product. Finally the vertical (i.e., thalamo-cortical) connectivity $K_{Vx}$ was added to this product to obtain a set of anatomo-functional connectivity matrices $K_x$ such that:

$$K_x = K_{DTI} \odot K_{Hx} + K_{Vx} \quad (1)$$

with $x = \{EXC, SST, BC, VIP, TRN1, TRN2\}$. All connectivity matrices are described in Supplementary Figure S3.

Using this combination of connectivity weights, and specific parameters defined in Supplementary Material for delta and background activity, we built a spatio-temporal source matrix S containing the time-varying activities of the thalamus and of all cortical macro-regions.

To reconstruct simulated scalp EEG data, we first solved the forward problem using the Boundary Element Method (BEM, OpenMEEG, Gramfort et al., 2010). To this end, a realistic head model was built in Brainstorm (Tadel et al., 2011) from the segmentation of a template T1 magnetic resonance imaging (MRI; Colin 27 template brain, Holmes et al., 1998) previously obtained using the Freesurfer image analysis suite<fnref n="1" /> (Dale et al., 1999; Fischl et al., 1999), as illustrated in Figure 3A. The head model consisted in three nested homogeneous mesh surfaces shaping the cortical surface (642 vertices), the skull (642 vertices) and the scalp (1,082 vertices) with conductivity values of $0.33 \text{ Sm}^{-1}$, $0.0082 \text{ Sm}^{-1}$ and $0.33 \text{ Sm}^{-1}$, respectively (Gonçalves et al., 2003). The forward problem was then numerically calculated for each vertex of the source mesh obtained from the segmented white matter/gray matter interface of the same template brain. As a result, the leadfield matrix A represented the contribution of each dipole of the source mesh at the level of 257 scalp electrode positions (high-density EEG), placed over the scalp according to the geodesic convention (EGI®), Eugene, USA). All leadfield vectors of A belonging to a common region of Desikan Atlas were added to obtain a simplified 66 × 257 matrix G. The spatio-temporal matrix X of simulated EEG data was given by the matrix product:

$$X = GS \quad (2)$$

The entire pipeline enabling the simulation of scalp EEG is provided in Figure 3. Finally, as detailed in the Supplementary Material, a total number of 1,060 ODE’s must be solved to run the model. To give an idea of the required computing time, 60 s of simulated EEG could be simulated in 49 s on a 3.5 GHz 6-Core Intel Xeon EG with 64 GB of 1,866 MHz RAM (OsX Mojave) and in 72 s on a standard Intel 2.5 GHz 2-Core Xeon EG with 8 GB RAM (Windows 10). Therefore, the model performed almost real-time simulation on a standard PC.

## Cortical Response to TMS and PCI

In order to simulate TMS-evoked EEG responses that can be compared to those recorded experimentally (Casali et al., 2013), we included the effect of an exogenous, TMS-induced, stimulation in the whole-brain model. The model can simulate

<footnote n="1">¹http://surfer.nmr.mgh.harvard.edu/</footnote>
<figure><img src="image_7.png" /><figcaption>**FIGURE 3** | Full processing pipeline leading to the simulation of scalp EEGs. The pipeline to simulate EEG is a two-step process. (A) First, the forward problem is solved at the level of 257 scalp electrodes from a dipole layer constrained to the surface of a cortical mesh (15,000 vertices) midway between the Gray/White matter interface and the cortex surface. The boundary Element Method is used for the calculation within a realistic head model that accounts for the conductivity properties and the geometry of brain, skull and scalp. This step provides a 257 × 15,000 leadfield matrix **A** representing the contribution of each individual cortical dipole at each of the 257 scalp electrodes. In this matrix, leadfield vectors belonging to a common region of the Desikan Atlas are added to obtain a simplified 66 × 257 matrix **G**. Second, the time courses **S** at the whole brain level are obtained in the mean-field model from a set of 66 cortical and one thalamic coupled NMMs. (B) Coupling between these 67 NMMs is done using combination of connection weight matrices. Pairs of structurally connected cortical NMMs are first defined from a matrix of connection weights representing a density of fibers between all pairs of 66 cortical regions of the Desikan Atlas. This matrix is provided in Hagmann et al. (2008). Using an element-wise multiplication, this matrix is combined with a set of horizontal (i.e., cortico-cortical) functional connectivity matrices that reproduce the coefficients weights used for wakefulness and sleep in the toy model. Vertical (i.e., thalamo-cortical) connectivity matrices are added to each of these products to obtain connectivity weight matrices that account for anatomical connections as well as cortico-cortical and thalamo-cortical connectivity matrices. Cortico- and thalamo-cortical time-delays were similarly organized in the form of matrices where the elements represent the Cartesian distance between cortical NMMs divided by the mean velocity of traveling for action potentials. (C) The mean-field model includes explicitly the contribution of an external stimulus term that represents the effect of TMS. At the output of the pipeline, scalp EEG signals at the level of 257 channels are obtained as the product of leadfield **G** and source time courses **S**.</figcaption></figure>

not only the activity of the 67 pre-defined anatomical regions but also the EEG signals recorded by 257 scalp electrodes. Furthermore, since the model includes, as a connectivity matrix between the 67 regions, a DTI-derived connectivity matrix (Hagmann et al., 2008), it is possible to track the spatio-temporal dynamics of the stimulation-evoked network, i.e., the activated regions, along with the peak latency for each region. In addition, since the model can simulate “wake” and “sleep” states, it provides the opportunity to compare TMS-evoked responses in these two states of consciousness, which have been experimentally recorded in humans. Therefore, the structure of our computational model provides a unique framework to interpret TMS-evoked EEG responses obtained in humans.

In principle, TMS involves a high-intensity current flowing through the stimulation coil, thereby generating a magnetic field penetrating without attenuation through the head. Since

the stimulation pulse is very short (≈0.1 ms), the magnetic field gradient dB/dt is extremely high (>30,000 T/s), resulting per Maxwell-Ampere's law into an electric field at the level of brain tissue. Since the electric field induced in brain tissue is high (>V/m, Miranda et al., 2003), this induces neuronal firing, presumably at the bending point of cortical axons, triggering a series of complex activations within the stimulated area (Di Lazzaro and Ziemann, 2013).

We used a simple approach to represent the effect of TMS, consisting in simulating an afferent volley of action potentials (in terms of pulses/s) to the stimulated cortical region. Since the 1 ms time step used to numerically solve the equations of the model was higher than the duration of an actual TMS pulse, the simulated length of this volley of incoming action potential was adapted and fixed to 5 ms. The amplitude of the simulated evoked volley of action potentials was fixed to 1,000 pulses/s and was applied to each cellular subtype of the stimulated region (i.e., PCs
and all types of GABAergic interneurons). For the purpose of this artilce, we chose to simulate the stimulation of the right motor cortex in both conditions (“wake” and “sleep”), which is known to have a number of anatomical connections that should result in the propagation of the TMS-evoked responses in several cortical structures. The repetition rate of the TMS protocol was fixed to 2 s, and the total simulation duration was a full minute, resulting in a total of 30 TMS-evoked responses available at the source level. The EEG activity at the level of each scalp electrode was then computed from the simulated source activity using our EEG forward problem pipeline, as described previously. From the simulated TMS-evoked responses, the studied outcomes were the anatomical regions activated by TMS, and also their latency from the onset of the TMS pulse.

Finally, simulated EEG responses were used to compute the Perturbational Complexity Index (PCI; Casali et al., 2013). PCI is a measure of TMS-evoked responses complexity, based on the Lampel-Ziv compression algorithm. The basic idea is that, if a sequence of information is complex, then it can be only marginally compressed (low compression rate); and on the opposite, a very simple sequence can be described by a very limited amount of information (high compression rate). To derive PCI, a process similar to the one proposed by Casali et al. (2013) was used:

* The sources activity is stored in a 2D matrix (number of lines equal to the number of channels, number of columns equal to the number of time points).

* A threshold value for the simulated sources activity was set by preserving the highest 20% (proportional threshold) of values once fixing to 0 values below it.

* Each value of the sources activity matrix that is above or equal to the threshold value is set to 1, while all other values are set to 0 (binarization process).

* The Lempel-Ziv algorithm is then applied to the resulting binary matrix.

PCI values were computed in the two different scenarios, corresponding to the “wake” and “sleep” state, respectively.

Overall, our implementation of TMS-evoked responses enables a meaningful comparison with human data, since it results in similar experimentally measurable quantities: impacted anatomical regions, latency of the TMS-evoked response within specific regions, and complexity of the brain-scale response through PCI.

# RESULTS

## Toy Model: The Impact of Cortico- and Thalamo-Cortical Connectivity on the Cortical Rhythms During Sleep and Wakefulness

During the deep sleep state, the thalamocortical connectivity is meant to be strong compared to the cortico-cortical one (see Supplementary Figure S2, first column). Our strategy consisted in progressively and simultaneously decreasing and increasing the thalamo- and cortico-cortical connectivity, respectively, in

order to switch to wakefulness (see Supplementary Figure S2, column 2). Note that this connectivity process was mostly reflected by arrays $K_{EXC}$ and $K_{BC}$, to simulate the strong thalamic projections onto PCs and BCs as reported in “Introduction” section. It is also noteworthy that no time delays were injected in the toy model example.

In Figure 4, we provide an example of simulated LFPs a.k.a. intracerebral EEG (iEEG) in comparison with real (human) ones. The left column depicts the scheme of high (and low) thalamo- (and cortico-) cortical connectivity. As depicted, the reinforcement of the thalamocortical loop (TC→PC→TRNs and TC→BC) resulted in the generation of delta oscillations that characterize SWS. In this example, the simulated delta was around 2-3 Hz (slightly faster than real iEEG). Conversely, by reducing the thalamocortical loop (right column), delta waves disappeared and were replaced by background activity, indicative of oscillatory changes observed during the switch from sleep to wake. We emphasize that only the large scale connectivity was tuned while all remaining parameters were kept unchanged, confirming the crucial role of thalamo- and cortico-cortical connectivity in modulating consciousness.

## Whole Brain Model: The Impact of Cortico- and Thalamo-Cortical Connectivity on Scalp EEG Rhythms During Sleep and Wakefulness

The morphology of simulated intracerebral signals was not modified when connectivity matrices were adapted to a larger number of NMMs in order to account for whole-brain activity. In the case of low thalamo-cortical connectivity (wakefulness condition), background activity was similar to signals obtained in the toy model and resembled real intracerebral background activity recorded during wakefulness in humans. Similarly, by increasing thalamo-cortical connectivity, the whole-brain model generated delta activity at a mean frequency of 3.8 Hz that was consistent with the morphology and spectral content of delta activity obtained both from the toy model and from real intracerebral recordings during SWS.

Signals obtained with the whole-brain model in both conditions of wakefulness and SWS were used in the forward calculation to generate simulated scalp EEG data at the level of 257 electrodes (Figure 5). In the low-thalamocortical connectivity condition, simulated scalp EEG (Figure 5A) resembled scalp EEG background activity as recorded in humans during wakefulness (Figure 5C). The spectral analysis disclosed similar sub-band distribution in the simulated vs. real case, although simulated signals contained more beta frequency than real background activity. In the high-thalamo-cortical connectivity condition, simulated scalp EEG (Figure 5B) was comparable to scalp EEG signals recorded in humans during SWS (Figure 5D). In the simulation case, the peak frequency was 3.8 Hz, thus slightly higher than in the real case (2 Hz). Topographical voltage maps at the peak of delta waves showed analogous distribution over the vertex, the activity in the simulated case being slightly more posterior than on the real case example.
<figure><img src="image_8.png" /><figcaption>**FIGURE 4** | Comparison of real and simulated intracerebral EEG (iEEG; toy model, *N* = 4). Left column: in the condition of high thalamo-cortical connectivity (i.e., low cortico-cortical connectivity), signals generated by the mean-field model are characterized by delta waves (~4 Hz). These simulated signals are similar (although slightly faster) to delta waves recorded by iEEG during slow wave sleep (SWS) in non-epileptic cortical regions of one patient undergoing an invasive EEG exploration. Right column: in the condition of low thalamo-cortical connectivity, signals generated by the mean-field model are similar to background activity recorded by iEEG in real conditions during wakefulness. Note that these iEEG recordings are performed in patients who are candidates to epilepsy surgery. For the sake of this study, only iEEG signals that do not show epileptic activity were retained.</figcaption></figure>

## Bridging Brain Circuits, TMS-Evoked EEG Responses and Complexity Metrics

In **Figure 6**, we describe the process used to estimate the complexity of TMS-evoked EEG responses within our brain-scale model and in humans (Casali et al., 2013). As depicted, the simulated TMS-evoked EEG response (**Figure 6B**) was very similar to the human response (**Figure 6A**), not only in terms of length (250 and 350 ms, respectively) but also in terms of regions distant from the right motor cortex (the stimulated area) that are activated post-stimulation. Indeed, in the simulated and experimental data, a first activation of the right pre-central gyrus occurred within 15 ms from the TMS pulse, followed by activity evoked notably in the right precuneus within 40 ms and a common propagation in the contralateral left precuneus at about 60 ms from the stimulus. In both simulated and experimental TMS-evoked responses, activity was evoked in the right frontal lobe within 110-120 ms, with activation of the right superior parietal cortex within 175 ms (human data) and 150 ms (model data). The most important difference is that the human TMS-evoked response was notably longer as compared to the model.

We then compared the simulated and experimental TMS-evoked responses in wakefulness and in sleep, along with their PCI value, which is presented in **Figure 7**. In the two cases, the right motor area is stimulated. In the case of wakefulness, as also illustrated in **Figure 6**, there was a satisfactory agreement in the duration and global shape of the response, and also in terms of the sequence of activated brain regions. More precisely, by comparing **Figure 7A** (upper) and **Figure 7B** (upper), we can conclude that the tracking of the propagation of the TMS-evoked activity revealed that this propagation occurred,

as expected, along documented inter-regional connections from the DTI-derived connectivity matrix used in the model. PCI values obtained were also extremely similar (0.52 and 0.51 for the model and for humans, respectively). Conversely, in the “Sleep” condition, the time course of the TMS-evoked response was significantly shorter (less than 200 ms), which is also comparable to TMS-EEG human recordings. In addition, another striking similarity with human data is that the TMS-evoked activity remained confined to the stimulated area, i.e., the right motor area. The PCI value in this condition was 0.19, which is very similar to the value obtained in humans during sleep (0.23, see **Figure 7A**, lower panel), notably lower than the “Wakefulness” condition, which is also consistent with human data (Casali et al., 2013). Therefore, despite using the exact same TMS pulse characteristics within the two conditions (“Wakefulness” and “Sleep”), the simulated response was drastically different within the model: in wakefulness, the TMS-evoked response resulted in a complex sequence of successive activations within distant, anatomically connected areas in ipsilateral and contralateral regions; while during sleep the TMS-evoked activity remained confined to the stimulation site, even when the anatomical connections were present as in the “sleep” condition.

## DISCUSSION AND PERSPECTIVE

In this article, we have developed the first brain-scale computational model that can reproduce neuronal activity patterns associated with various conscious states, while accounting for key microcircuits at the cellular type scale. A major asset of the model is its strong link with recent neurophysiological and neuroanatomical data: the main cellular
<figure><img src="image_9.png" /><figcaption>**FIGURE 5** | Simulated vs. Real EEG during wakefulness and SWS. Signals simulated with the whole-brain model using weak thalamo-cortical connectivity parameters (A) display background activity. The morphology and spectral content of these simulated signals are similar to scalp EEG recorded in a human subject during wakefulness in humans (C), except for a higher power spectral density in the beta sub-band. Signals simulated with the whole-brain model using strong thalamo-cortical connectivity parameters (B) display delta waves similar to the activity recorded in real condition during SWS in humans (D). The spectral content of signals as well as the topographical voltage distribution at the peak of delta waves were similar in the simulated and real conditions.</figcaption></figure>

types are included (PCs and different types of interneurons), along with their recently elucidated connections that underlie the selective disinhibition of distant neural populations (through

VIP to SST projections), realistic synaptic kinetics, large-scale structural connectivity obtained in humans through DTI and propagation delays between regions based on their spatial
**A** Real, wakefulness +TMS stimulation (modified from Casali et al, 2013)

<figure><img src="image_10.png" /></figure>

**B** Virtual Brain, « wakefulness mode » + Simulated TMS

<figure><img src="image_11.png" /></figure>

**FIGURE 6** | Comparison of TMS-EEG evoked responses in silico and in humans. (A) Time course of a human TMS-EEG response (modified from Casali et al., 2013) following stimulation of the motor cortex during wakefulness. Once that cortical sources have been computed from EEG recordings, a spatio-temporal matrix of significant sources was built and the Lempel-Ziv compression algorithm was used to compute the complexity of the evoked response (Perturbational Complexity Index, PCI). (B) Time course of a simulated TMS-EEG response using our brain-scale, following stimulation of the motor area in the wakefulness mode. Cortical sources were reconstructed from the simulated EEG, and a similar procedure was used to compute PCI. A similar PCI value was obtained in the simulated and experimental TMS-evoked EEG responses in the awake state.

distance. Furthermore, since the model features detailed micro-circuits, DTI-derived brain connectivity matrix and large-scale macro-circuits results also provided a bottom-up description of TMS-evoked responses in humans.

The model offers novel, key insights into the maintenance of the neuronal activity associated with conscious states. Interestingly, this model including a variety of cellular subtypes with accurate synaptic kinetics provides realistic
<figure><img src="image_12.png" /><figcaption>Real, NREM sleep</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**B** Virtual brain, wakefulness mode</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Virtual brain, NREM sleep mode</figcaption></figure>

<figure><img src="image_15.png" /></figure>

**FIGURE 7** | Comparison of TMS-evoked EEG responses in wakefulness and sleep. (A) TMS-evoked EEG responses obtained through TMS of the motor cortex in humans (upper panel, during wakefulness; lower panel, during sleep) with associated PCI values. (B) TMS-evoked EEG responses obtained through simulated TMS of the motor region (upper panel, in the wakefulness mode; lower panel, during the sleep mode) with associated PCI values.

electrophysiological signals both at the level of cortical sources and of the EEG. First, it explains how thalamo-cortical (vertical) connectivity is critically involved in the gating of cortico-cortical (horizontal) information propagation. If thalamo-cortical activity is indeed rhythmically patterned, the communication between cortical areas is disrupted due to the resulting rhythmic inhibition. This result is therefore in line with the “connectivity breakdown” observed during sleep (Esser et al., 2009; Ferrarelli et al., 2010; Casali et al., 2013).

Crucially, the various conscious states simulated with the model can be tuned by modulating only the thalamo-cortical input, which supports our hypothesis that this represents the crucial control parameter for consciousness, i.e., the dimension of wakefulness. It can be seen essentially as the possibility for information processing to take place between spatially distant cortical areas. Since the propagation of activity is impaired at the cortical level during sleep (absence of wakefulness), no cortical processing can take place, which explains the absence of consciousness. Second, regarding cortico-cortical activity (horizontal connectivity), since the model reproduces with a satisfactory qualitative and quantitative agreement the TMS-evoked EEG responses observed in awake humans (Casali et al., 2013), this suggests that: (1) our modeling hypotheses and choices appear sufficient to capture the essence of TMS-EEG responses; (2) TMS-EEG responses are mainly driven by the underlying structural connectome, which is in

line with recent research pointing at the tight links between structural and functional connectivity (Avena-Koenigsberger et al., 2017); and (3) this modeling approach could be used to assist in the interpretation of TMS-evoked EEG responses in DOC patients, since our model links explicitly the underlying connectivity with the observed TMS-evoked response. Finally, the model validates that, at the brain scale, the disynaptic disinhibition of distant PCs through the activation of VIP neurons is indeed an effective mechanism enabling the transmission of activity between a few cortical regions. The model therefore confirms the role of a cellular-scale micro-circuit that regulates brain-scale propagation of activity within the cortex.

Among the possibilities to improve the realism and predictive power of our brain-scale model, the most immediate would be the use of structural connectivity matrices averaged among a large number of healthy participants, such as those from the Human Connectome Project<fnref n="2" />. Furthermore, no regional specificities were accounted for between the 66 cortical regions included within the model, since we used standard parameter values for the synaptic gains (e.g., A, B and G) and the same within-population connectivity parameters. By doing so, we have assumed that the large-scale anatomic structure of brain connectivity and the cellular-scale micro-circuits included

<footnote n="2">²http://www.humanconnectomeproject.org</footnote>
are the main factors explaining the resulting simulated EEG
signals. Another limitation of COALIA under its current form
is the simplification of the physiological processes involved
in the transition from wakefulness to sleep. The brainstem
is indeed involved in the modulation of thalamic activity,
and the regulation of cortico-cortical connectivity (Crunelli
et al., 2015). In the model, the brainstem input onto the
thalamus is not explicitly included. Instead, the modulation of
thalamic inhibition between the sleep/wake mode was simulated
by changing the excitatory input onto the RT-interneurons
or TC cells.

It should be mentioned that in the case of deep sleep
simulated signals, the median frequency of delta activity was
higher in the model (2–3 Hz) as compared to human data
(~1 Hz). This discrepancy is explained by the fact that model
does not implement the mechanisms underlying slow oscillations
(~1 Hz) generated in cortical and thalamic networks, among
which: (i) the sequence of depolarizing periods followed by
silent periods (up to 2 s) during up-and-down states as well as;
and (ii) GABA<sub>B</sub>-mediated pre-synaptic slow inhibition that also
appears to play some role (see review in Neske, 2015). Other
limitations likely explain the moderate discrepancies between
the simulated and experimental TMS-EEG responses in terms
of latencies and localizations, such as the lack of asymmetry
in connectivity weights (i.e., all connections were assumed
bidirectional and identical).

Let us emphasize that the results presented in the present
article are only a first step towards understanding the transition
between wake and sleep and that COALIA provides an
appropriate framework to achieve that objective. The future
prospects regarding our brain-scale EEG model are numerous: in
terms of consciousness studies, the model could be used to study
the mechanisms underlying the so-called “slow-wave activity
saturation,” delta-band activity that appears when the blood
concentration of anesthetics is increased (Ni Mhuircheartaigh
et al., 2013), and that constitutes a solid marker of the conscious
state. Furthermore, this study improves our understanding of
active probing paradigms of brain circuits in DOCs, such as
the PCI and paves the way toward the design of optimized

stimulation-based metrics to measure consciousness. The model
indeed allows to test *in silico* novel neuromodulation protocols
based on TMS, transcranial direct current stimulation (tDCS)
and transcranial alternating current stimulation (tACS), aiming
at quantifying the level of residual consciousness in DOC
patients. Beyond applications for consciousness, the model could
be exploited to understand the detailed dynamics of TMS-EEG
responses and their underlying mechanisms, or to shed light on
the mechanisms underlying the generation and propagation of
epileptiform activity.

DATA AVAILABILITY STATEMENT

The raw data supporting the conclusions of this manuscript will
be made available by the authors, without undue reservation, to
any qualified researcher.

AUTHOR CONTRIBUTIONS

SB: model implementation, simulations, data analysis, manuscript write-up. JM: model design, simulations, data analysis, manuscript write-up. IM: simulations, data analysis, manuscript write-up. FW and PB: model design, data analysis, manuscript write-up.

FUNDING

This work has been fully funded by the LUMINOUS Project. This
project has received funding from the European Union's Horizon
2020 research and innovation program H2020-FETOPEN-
2014-2015-RIA under agreement No. 686764. This study was
supported in part by the Institut National de la Santé et de la
Recherche Médicale (INSERM, France).

SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online
at: https://www.frontiersin.org/articles/10.3389/fnsys.2019.000
59/full#supplementary-material.

REFERENCES

Adesnik, H., Bruns, W., Taniguchi, H., Huang, Z. J., and Scanziani, M. (2012). A neural circuit for spatial summation in visual cortex. *Nature* 490, 226–231. doi: 10.1038/nature11526

Avena-Koenigsberger, A., Misic, B., and Sporns, O. (2017). Communication dynamics in complex brain networks. *Nat. Rev. Neurosci.* 19, 17–33. doi: 10.1038/nrn.2017.149

Ayzenshtat, I., Karnani, M. M., Jackson, J., and Yuste, R. (2016). Cortical control of spatial resolution by VIP⁺ interneurons. *J. Neurosci.* 36, 11498–11509. doi: 10.1523/jneurosci.1920-16.2016

Bhattacharya, B. S., Coyle, D., and Maguire, L. P. (2011). A thalamo-cortico-
thalamic neural mass model to study α rhythms in Alzheimer's disease. *Neural
Netw.* 24, 631–645. doi: 10.1016/j.neunet.2011.02.009

Breakspear, M. (2017). Dynamic models of large-scale brain activity. *Nat. Neurosci.* 20, 340–352. doi: 10.1038/nn.4497

Buzsáki, G., and Wang, X.-J. (2012). Mechanisms of γ oscillations. *Annu. Rev. Neurosci.* 35, 203–225. doi: 10.1146/annurev-neuro-062111-150444

Casali, A. G., Gosseries, O., Rosanova, M., Boly, M., Sarasso, S., Casali, K. R., et al. (2013). A theoretically based index of consciousness independent of sensory processing and behavior. *Sci. Transl. Med.* 5:198ra105. doi: 10.1126/scitranslmed.3006294

Chen, G., Zhang, Y., Li, X., Zhao, X., Ye, Q., Lin, Y., et al. (2017). Distinct inhibitory circuits orchestrate cortical β and γ band oscillations. *Neuron* 96, 1403.e6–1418.e6. doi: 10.1016/j.neuron.2017.11.033

Cona, F., Lacanna, M., and Ursino, M. (2014). A thalamo-cortical neural mass model for the simulation of brain rhythms during sleep. *J. Comput. Neurosci.* 37, 125–148. doi: 10.1007/s10827-013-0493-1

Crone, J. S., Schurz, M., Höller, Y., Bergmann, J., Monti, M., Schmid, E., et al. (2015). Impaired consciousness is linked to changes in effective connectivity of the posterior cingulate cortex within the default mode network. *Neuroimage* 110, 101–109. doi: 10.1016/j.neuroimage.2015.01.037

Cruikshank, S. J., Lewis, T. J., and Connors, B. W. (2007). Synaptic basis for intense thalamocortical activation of feedforward inhibitory cells in neocortex. *Nat. Neurosci.* 10, 462–468. doi: 10.1038/nn1861
Crunelli, V., David, F., Lorincz, M. L., and Hughes, S. W. (2015). The thalamocortical network as a single slow wave-generating unit. *Curr. Opin. Neurobiol.* 31, 72–80. doi: 10.1016/j.conb.2014.09.001

Dale, A. M., Fischl, B., and Sereno, M. I. (1999). Cortical surface-based analysis. *I. Segmentation and surface reconstruction*. *Neuroimage* 9, 179–194. doi: 10.1006/nimg.1998.0395

David, O., and Friston, K. J. (2003). A neural mass model for MEG/EEG: coupling and neuronal dynamics. *Neuroimage* 20, 1743–1755. doi: 10.1016/j.neuroimage.2003.07.015

Dehaene, S., and Changeux, J. P. (2011). Experimental and theoretical approaches to conscious processing. *Neuron* 70, 200–227. doi: 10.1016/j.neuron.2011.03.018

Dehaene, S., Naccache, L., Cohen, L., Bihan, D. L., Mangin, J. F., Poline, J. B., et al. (2001). Cerebral mechanisms of word masking and unconscious repetition priming. *Nat. Neurosci.* 4, 752–758. doi: 10.1038/89551

Dehaene, S., Kerszberg, M., and Changeux, J. P. (1998). A neuronal model of a global workspace in effortful cognitive tasks. *Proc. Natl. Acad. Sci. U S A* 95, 14529–14534. doi: 10.1073/pnas.95.24.14529

Desikan, R. S., Segonne, F., Fischl, B., Quinn, B. T., Dickerson, B. C., Blacker, D., et al. (2006). An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral regions of interest. *Neuroimage* 31, 968–980. doi: 10.1016/j.neuroimage.2006.01.021

Di Lazzaro, V., and Ziemann, U. (2013). The contribution of transcranial magnetic stimulation in the functional evaluation of microcircuits in human motor cortex. *Front. Neural Circuits* 7:18. doi: 10.3389/fncir.2013.00018

Dipoppa, M., Ranson, A., Krumin, M., Pachitariu, M., Carandini, M., and Harris, K. D. (2018). Vision and locomotion shape the interactions between neuron types in mouse visual cortex. *Neuron* 98, 602.e8–615.e8. doi: 10.1016/j.neuron.2018.03.037

Esser, S. K., Hill, S., and Tononi, G. (2009). Breakdown of effective connectivity during slow wave sleep: investigating the mechanism underlying a cortical gate using large-scale modeling. *J. Neurophysiol.* 102, 2096–2111. doi: 10.1152/jn.00059.2009

Ferrarelli, F., Massimini, M., Sarasso, S., Casali, A., Riedner, B. A., Angelini, G., et al. (2010). Breakdown in cortical effective connectivity during midazolam-induced loss of consciousness. *Proc. Natl. Acad. Sci. U S A* 107, 2681–2686. doi: 10.1073/pnas.0913008107

Fino, E., and Yuste, R. (2011). Dense inhibitory connectivity in neocortex. *Neuron* 69, 1188–1203. doi: 10.1016/j.neuron.2011.02.025

Fischl, B., Sereno, M. I., and Dale, A. M. (1999). Cortical surface-based analysis. *II: Inflation, flattening, and a surface-based coordinate system*. *Neuroimage* 9, 195–207. doi: 10.1006/nimg.1998.0396

Freeman, W. J. (1978). Models of the dynamics of neural populations. *Electroencephalogr. Clin. Neurophysiol. Suppl.* 34, 9–18.

Fries, P. (2005). A mechanism for cognitive dynamics: neuronal communication through neuronal coherence. *Trends Cogn. Sci.* 9, 474–480. doi: 10.1016/j.tics.2005.08.011

Fries, P. (2009). Neuronal γ-band synchronization as a fundamental process in cortical computation. *Annu. Rev. Neurosci.* 32, 209–224. doi: 10.1146/annurev.neuro.051508.135603

Friston, K. J. (2011). Functional and effective connectivity: a review. *Brain Connect.* 1, 13–36. doi: 10.1089/brain.2011.0008

Fu, Y., Tucciarone, J. M., Espinosa, J. S., Sheng, N., Darcy, D. P., Nicoll, R. A., et al. (2014). A cortical circuit for gain control by behavioral state. *Cell* 156, 1139–1152. doi: 10.1016/j.cell.2014.01.050

Funk, C. M., Peelman, K., Bellesi, M., Marshall, W., Cirelli, C., and Tononi, G. (2017). Role of somatostatin-positive cortical interneurons in the generation of sleep slow waves. *J. Neurosci.* 37, 9132–9148. doi: 10.1523/jneurosci.1303–17.2017

Gentet, L. J., Kremer, Y., Taniguchi, H., Huang, Z. J., Staiger, J. F., and Petersen, C. C. H. (2012). Unique functional properties of somatostatin-expressing GABAergic neurons in mouse barrel cortex. *Nat. Neurosci.* 15, 607–612. doi: 10.1038/nn.3051

Gidon, A., and Segev, I. (2012). Principles governing the operation of synaptic inhibition in dendrites. *Neuron* 75, 330–341. doi: 10.1016/j.neuron.2012.05.015

Gómez, F., Phillips, C., Soddu, A., Boly, M., Boveroux, P., Vanhaudenhuyse, A., et al. (2013). Changes in effective connectivity by propofol sedation. *PLoS One* 8:e71370. doi: 10.1371/journal.pone.0071370

Gonçalves, S. I., de Munck, J. C., Verbunt, J. P., Bijma, F., Heethaar, R. M., and Lopes da Silva, F. (2003). *In vivo* measurement of the brain and skull resistivities using an EIT-based method and realistic models for the head. *IEEE Trans. Biomed. Eng.* 50, 754–767. doi: 10.1109/tbme.2003.812164

Gramfort, A., Papadopoulo, T., Olivi, E., and Clerc, M. (2010). OpenMEEG: opensource software for quasistatic bioelectromagnetics. *Biomed. Eng. Online* 9:45. doi: 10.1186/1475-925x-9-45

Hagmann, P., Cammoun, L., Gigandet, X., Meuli, R., Honey, C. J., Wedeen, V. J., et al. (2008). Mapping the structural core of human cerebral cortex. *PLoS Biol.* 6:e159. doi: 10.1371/journal.pbio.0060159

Hale, J. R., White, T. P., Mayhew, S. D., Wilson, R. S., Rollings, D. T., Khalsa, S., et al. (2016). Altered thalamocortical and intra-thalamic functional connectivity during light sleep compared with wake. *Neuroimage* 125, 657–667. doi: 10.1016/j.neuroimage.2015.10.041

Harris, A. Z., and Gordon, J. A. (2015). Long-range neural synchrony in behavior. *Annu. Rev. Neurosci.* 38, 171–194. doi: 10.1146/annurev-neuro-071714-034111

Harris, K. D., and Mrsic-Flogel, T. D. (2013). Cortical connectivity and sensory coding. *Nature* 503, 51–58. doi: 10.1038/nature12654

Hill, S., and Tononi, G. (2005). Modeling sleep and wakefulness in the thalamocortical system. *J. Neurophysiol.* 93, 1671–1698. doi: 10.1152/jn.00915.2004

Holmes, C. J., Hoge, R., Collins, L., Woods, R., Toga, A. W., and Evans, A. C. (1998). Enhancement of MR images using registration for signal averaging. *J. Comput. Assist. Tomogr.* 22, 324–333. doi: 10.1097/00004728-199803000-00032

Hudetz, A. G. (2012). General anesthesia and human brain connectivity. *Brain Connect.* 2, 291–302. doi: 10.1089/brain.2012.0107

Jansen, B. H., and Rit, V. G. (1995). Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. *Biol. Cybern.* 73, 357–366. doi: 10.1007/bf00199471

Ji, X.-Y., Zingg, B., Mesik, L., Xiao, Z., Zhang, L. I., and Tao, H. W. (2016). Thalamocortical innervation pattern in mouse auditory and visual cortex: laminar and cell-type specificity. *Cereb. Cortex* 26, 2612–2625. doi: 10.1093/cercor/bhv099

Jin, S.-H., and Chung, C. K. (2012). Messages from the brain connectivity regarding neural correlates of consciousness. *Exp. Neurobiol.* 21, 113–122. doi: 10.5607/en.2012.21.3.113

Kapfer, C., Glickfeld, L. L., Atallah, B. V., and Scanziani, M. (2007). Supralinear increase of recurrent inhibition during sparse activity in the somatosensory cortex. *Nat. Neurosci.* 10, 743–753. doi: 10.1038/nn1909

Karnani, M. M., Agetsuma, M., and Yuste, R. (2014). A blanket of inhibition: functional inferences from dense inhibitory connectivity. *Curr. Opin. Neurobiol.* 26, 96–102. doi: 10.1016/j.conb.2013.12.015

Kim, S.-P., Hwang, E., Kang, J.-H., Kim, S., and Choi, J. H. (2012). Changes in the thalamocortical connectivity during anesthesia-induced transitions in consciousness. *Neuroreport* 23, 294–298. doi: 10.1097/wnr.0b013e3283509ba0

Koch, C., Massimini, M., Boly, M., and Tononi, G. (2016). Neural correlates of consciousness: progress and problems. *Nat. Rev. Neurosci.* 17, 307–321. doi: 10.1038/nrn.2016.61

Laureys, S. (2004). Functional neuroimaging in the vegetative state. *NeuroRehabilitation* 19, 335–341.

Laureys, S., Boly, M., and Maquet, P. (2006). Tracking the recovery of consciousness from coma. *J. Clin. Invest.* 116, 1823–1825. doi: 10.1172/jci29172

Laureys, S., Goldman, S., Phillips, C., Van Bogaert, P., Aerts, J., Luxen, A., et al. (1999). Impaired effective cortical connectivity in vegetative state: preliminary investigation using PET. *Neuroimage* 9, 377–382. doi: 10.1006/nimg.1998.0414

Lee, S., Kruglikov, I., Huang, Z. J., Fishell, G., and Rudy, B. (2013). A disinhibitory circuit mediates motor integration in the somatosensory cortex. *Nat. Neurosci.* 16, 1662–1670. doi: 10.1038/nn.3544

Lee, J. H., and Mihalas, S. (2017). Visual processing mode switching regulated by VIP cells. *Sci. Rep.* 7:1843. doi: 10.1038/s41598-017-01830-0
Lewis, D. A., Curley, A. A., Glausier, J. R., and Volk, D. W. (2012). Cortical parvalbumin interneurons and cognitive dysfunction in schizophrenia. *Trends Neurosci.* 35, 57–67. doi: 10.1016/j.tins.2011.10.004

Lewis, L. D., Voigts, J., Flores, F. J., Schmitt, L. I., Wilson, M. A., Halassa, M. M., et al. (2015). Thalamic reticular nucleus induces fast and local modulation of arousal state. *Elife* 4:e08760. doi: 10.7554/efelife.08760

Liu, F., Wang, J., Liu, C., Li, H., Deng, B., Fietkiewicz, C., et al. (2016). A neural mass model of basal ganglia nuclei simulates pathological β rhythm in Parkinson's disease. *Chaos* 26:123113. doi: 10.1063/1.4972200

Liu, C., Zhu, Y., Liu, F., Wang, J., Li, H., Deng, B., et al. (2017). Neural mass models describing possible origin of the excessive β oscillations correlated with Parkinsonian state. *Neural Netw.* 88, 65–73. doi: 10.1016/j.neunet.2017.01.011

Llinás, R., Ribary, U., Contreras, D., and Pedroarena, C. (1998). The neuronal basis for consciousness. *Philos. Trans. R. Soc. Lond. B Biol. Sci.* 353, 1841–1849. doi: 10.1098/rstb.1998.0336

Lopes da Silva, F. H., van Rotterdam, A., Barts, P., van Heusden, E., and Burr, W. (1976). Models of neuronal populations: the basic mechanisms of rhythmicity. *Prog. Brain Res.* 45, 281–308. doi: 10.1016/s0079-6123(08)60995-4

Lovett-Barron, M., Turi, G. F., Kaifosh, P., Lee, P. H., Bolze, F., X.-Sun, H., et al. (2012). Regulation of neuronal input transformations by tunable dendritic inhibition. *Nat. Neurosci.* 15, 423–430. doi: 10.1038/nn.3024

Ma, W.-P., Liu, B.-H., Li, Y.-T., Huang, Z. J., Zhang, L. I., and Tao, H. W. (2010). Visual representations by cortical somatostatin inhibitory neurons—selective but with weak and delayed responses. *J. Neurosci.* 30, 14371–14379. doi: 10.1523/JNEUROSCI.3248-10.2010

Maex, R., and De Schutter, E. (2007). Mechanism of spontaneous and self-sustained oscillations in networks connected through axo-axonal gap junctions. *Eur. J. Neurosci.* 25, 3347–3358. doi: 10.1111/j.1460-9568.2007.05593.x

Massimini, M., Ferrarelli, F., Murphy, M. J., Huber, R., Riedner, B. A., Casarotto, S., et al. (2010). Cortical reactivity and effective connectivity during REM sleep in humans. *Cogn. Neurosci.* 1, 176–183. doi: 10.1080/17588921003731578

Mina, F., Benquet, P., Pasnicu, A., Biraben, A., and Wendling, F. (2013). Modulation of epileptic activity by deep brain stimulation: a model-based study of frequency-dependent effects. *Front. Comput. Neurosci.* 7:94. doi: 10.3389/fncom.2013.00094

Miranda, P. C., Hallett, M., and Basser, P. J. (2003). The electric field induced in the brain by magnetic stimulation: a 3-D finite-element analysis of the effect of tissue heterogeneity and anisotropy. *IEEE Trans. Biomed. Eng.* 50, 1074–1085. doi: 10.1109/tbme.2003.816079

Modolo, J., Hassan, M., Ruffini, G., and Legros, A. (2018). Probing the circuits of conscious perception with magnetophosphenes. *Preprint*

Molaee-Ardekani, B., Benquet, P., Bartolomei, F., and Wendling, F. (2010). Computational modeling of high-frequency oscillations at the onset of neocortical partial seizures: from 'altered structure' to 'dysfunction'. *Neuroimage* 52, 1109–1122. doi: 10.1016/j.neuroimage.2009.12.049

Nakajima, M., and Halassa, M. M. (2017). Thalamic control of functional cortical connectivity. *Curr. Opin. Neurobiol.* 44, 127–131. doi: 10.1016/j.conb.2017.04.001

Neske, G. T. (2015). The slow oscillation in cortical and thalamic networks: mechanisms and functions. *Front. Neural Circuits* 9:88. doi: 10.3389/fncir.2015.00088

Ni Mhuircheartaigh, R., Warnaby, C., Rogers, R., Jbabdi, S., and Tracey, I. (2013). Slow-wave activity saturation and thalamocortical isolation during propofol anesthesia in humans. *Sci. Transl. Med.* 5:208ra148. doi: 10.1126/scitranslmed.3006007

Nunez, P. L. (1974). The brain wave equation: a model for the EEG. *Math. Biosci.* 21, 279–297. doi: 10.1016/0025-5564(74)90020-0

Owen, A. M., Coleman, M. R., Boly, M., Davis, M. H., Laureys, S., and Pickard, J. D. (2006). Detecting awareness in the vegetative state. *Science* 313:1402. doi: 10.1126/science.1130197

Paulus, W., and Rothwell, J. C. (2016). Membrane resistance and shunting inhibition: where biophysics meets state-dependent human neurophysiology. *J. Physiol.* 594, 2719–2728. doi: 10.1113/jp271452

Pfeffer, C. K. (2014). Inhibitory neurons: vip cells hit the brake on inhibition. *Curr. Biol.* 24, R18–R20. doi: 10.1016/j.cub.2013.11.001

Pi, H.-J., Hangya, B., Kvitsiani, D., Sanders, J. I., Huang, Z. J., and Kepecs, A. (2013). Cortical interneurons that specialize in disinhibitory control. *Nature* 503, 521–524. doi: 10.1038/nature12676

Povysheva, N. V., Zaitsev, A. V., Rotaru, D. C., Gonzalez-Burgos, G., Lewis, D. A., and Krimer, L. S. (2008). Parvalbumin-positive basket interneurons in monkey and rat prefrontal cortex. *J. Neurophysiol.* 100, 2348–2360. doi: 10.1152/jn.90396.2008

Roberts, J. A., and Robinson, P. A. (2012). Corticothalamic dynamics: structure of parameter space, spectra, instabilities, and reduced model. *Phys. Rev. E Stat. Nonlin. Soft Matter Phys.* 85:011910. doi: 10.1103/physreve.85.011910

Rudy, B., Fishell, G., Lee, S., and Hjerling-Leffler, J. (2011). Three groups of interneurons account for nearly 100% of neocortical GABAergic neurons. *Dev. Neurobiol.* 71, 45–61. doi: 10.1002/dneu.20853

Ruffini, G. (2017). An algorithmic information theory of consciousness. *Neurosci. Conscious.* 2017:nix019. doi: 10.1093/nc/nix019

Saalmann, Y. B. (2014). Intralaminar and medial thalamic influence on cortical synchrony, information transmission and cognition. *Front. Syst. Neurosci.* 8:83. doi: 10.3389/fnsys.2014.00083

Schoenberg, P. L. A., Ruf, A., Churchill, J., Brown, D. P., and Brewer, J. A. (2018). Mapping complex mind states: EEG neural substrates of meditative unified compassionate awareness. *Conscious. Cogn.* 57, 41–53. doi: 10.1016/j.concog.2017.11.003

Sen Bhattacharya, B., Cakir, Y., Serap-Sengor, N., Maguire, L., and Coyle, D. (2013). Model-based bifurcation and power spectral analyses of thalamocortical α rhythm slowing in Alzheimer's disease. *Neurocomputing* 115, 11–22. doi: 10.1016/j.neucom.2012.10.023

Silberberg, G., and Markram, H. (2007). Disynaptic inhibition between neocortical pyramidal cells mediated by martinotti cells. *Neuron* 53, 735–746. doi: 10.1016/j.neuron.2007.02.012

Sohn, J., Okamoto, S., Kataoka, N., Kaneko, T., Nakamura, K., and Hioki, H. (2016). Differential inputs to the perisomatic and distal-dendritic compartments of VIP-positive neurons in layer 2/3 of the mouse barrel cortex. *Front. Neuroanat.* 10:124. doi: 10.3389/fnana.2016.00124

Sotero, R. C., Trujillo-Barreto, N. J., Iturria-Medina, Y., Carbonell, F., and Jimenez, J. C. (2007). Realistically coupled neural mass models can generate EEG rhythms. *Neural Comput.* 19, 478–512. doi: 10.1162/neco.2007.19.2.478

Suffczynski, P., Kalitzin, S., and Lopes Da Silva, F. H. (2004). Dynamics of non-convulsive epileptic phenomena modeled by a bistable neuronal network. *Neuroscience* 126, 467–484. doi: 10.1016/j.neuroscience.2004.03.014

Tadel, F., Baillet, S., Mosher, J. C., Pantazis, D., and Leahy, R. M. (2011). Brainstorm: a user-friendly application for MEG/EEG analysis. *Comput. Intell. Neurosci.* 2011:879716. doi: 10.1155/2011/879716

Tan, Z., Hu, H., Huang, Z. J., and Agmon, A. (2008). Robust but delayed thalamocortical activation of dendritic-targeting inhibitory interneurons. *Proc. Natl. Acad. Sci. USA* 105, 2187–2192. doi: 10.1073/pnas.0710628105

Timofeev, I., and Steriade, M. (1996). Low-frequency rhythms in the thalamus of intact-cortex and decorticated cats. *J. Neurophysiol.* 76, 4152–4168. doi: 10.1152/jn.1996.76.6.4152

Tononi, G. (2004). An information integration theory of consciousness. *BMC Neurosci.* 5:42. doi: 10.1186/1471-2202-5-42

Tononi, G. (2012). Integrated information theory of consciousness: an updated account. *Arch. Ital. Biol.* 150, 56–90. doi: 10.4449/aib.v149i5.1388

Tononi, G., and Edelman, G. M. (1998). Consciousness and complexity. *Science* 282, 1846–1851. doi: 10.1126/science.282.5395.1846

Tononi, G., and Sporns, O. (2003). Measuring information integration. *BMC Neurosci.* 4:31. doi: 10.1186/1471-2202-4-31

Traub, R. D., Contreras, D., Cunningham, M. O., Murray, H., LeBeau, F. E. N., Roopun, A., et al. (2005). Single-column thalamocortical network model exhibiting γ oscillations, sleep spindles, and epileptogenic bursts. *J. Neurophysiol.* 93, 2194–2232. doi: 10.1152/jn.01147.2004

Tremblay, R., Lee, S., and Rudy, B. (2016). GABAergic interneurons in the neocortex: from cellular properties to circuits. *Neuron* 91, 260–292. doi: 10.1016/j.neuron.2016.06.033

Urban-Ciecko, J., and Barth, A. L. (2016). Somatostatin-expressing neurons in cortical networks. *Nat. Rev. Neurosci.* 17, 401–409. doi: 10.1038/nrn.2016.53

Ursino, M., Cona, F., and Zavaglia, M. (2010). The generation of rhythms within a cortical region: analysis of a neural mass model. *Neuroimage* 52, 1080–1094. doi: 10.1016/j.neuroimage.2009.12.084
Usami, K., Matsumoto, R., Kobayashi, K., Hitomi, T., Shimotake, A., Kikuchi, T., et al. (2015). Sleep modulates cortical connectivity and excitability in humans: direct evidence from neural activity induced by single-pulse electrical stimulation. *Hum. Brain Mapp.* 36, 4714–4729. doi: 10.1002/hbm.22948

van Vugt, B., Dagnino, B., Vartak, D., Safaai, H., Panzeri, S., Dehaene, S., et al. (2018). The threshold for conscious report: signal loss and response bias in visual and frontal cortex. *Science* 360, 537–542. doi: 10.1126/science.aar7186

Varga, C., Oijala, M., Lish, J., Szabo, G. G., Bezaire, M., Marchionni, I., et al. (2014). Functional fission of parvalbumin interneuron classes during fast network events. *eLife* 3:e04006. doi: 10.7554/elife.04006

Varotto, G., Fazio, P., Rossi Sebastiano, D., Duran, D., D'Incerti, L., Parati, E., et al. (2014). Altered resting state effective connectivity in long-standing vegetative state patients: an EEG study. *Clin. Neurophysiol.* 125, 63–68. doi: 10.1016/j.clinph.2013.06.016

Walker, F., Möck, M., Feyerabend, M., Guy, J., Wagener, R. J., Schubert, D., et al. (2016). Parvalbumin- and vasoactive intestinal polypeptide-expressing neocortical interneurons impose differential inhibition on Martinotti cells. *Nat. Commun.* 7:13664. doi: 10.1038/ncomms13664

Wang, X.-J., and Buzsáki, G. (1996). γ oscillation by synaptic inhibition in a hippocampal interneuronal network model. *J. Neurosci.* 16, 6402–6413. doi: 10.1523/JNEUROSCI.16-20-06402.1996

Ward, L. M. (2013). The thalamus: gateway to the mind. *Wiley Interdiscip. Rev. Cogn. Sci.* 4, 609–622. doi: 10.1002/wcs.1256

Wendling, F., Bartolomei, F., Bellanger, J. J., and Chauvel, P. (2002). Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition. *Eur. J. Neurosci.* 15, 1499–1508. doi: 10.1046/j.1460-9568.2002.01985.x

Wendling, F., Benquet, P., Bartolomei, F., and Jirsa, V. (2016). Computational models of epileptiform activity. *J. Neurosci. Methods* 260, 233–251. doi: 10.1016/j.jneumeth.2015.03.027

Whittington, M. A., Traub, R. D., Kopell, N., Ermentrout, B., and Buhl, E. H. (2000). Inhibition-based rhythms: experimental and mathematical

observations on network dynamics. *Int. J. Psychophysiol.* 38, 315–336. doi: 10.1016/s0167-8760(00)00173-2

Williams, L. E., and Holtmaat, A. (2019). Higher-order thalamocortical inputs gate synaptic long-term potentiation *via* disinhibition. *Neuron* 101, 91.e4–102.e4. doi: 10.1016/j.neuron.2018.10.049

Wilson, H. R., and Cowan, J. D. (1973). A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. *Kybernetik* 13, 55–80. doi: 10.1007/bf00288786

Womelsdorf, T., Valiante, T. A., Sahin, N. T., Miller, K. J., and Tiesinga, P. (2014). Dynamic circuit motifs underlying rhythmic gain control, gating and integration. *Nat. Neurosci.* 17, 1031–1039. doi: 10.1038/nn.3764

Yang, W., Carrasquillo, Y., Hooks, B. M., Nerbonne, J. M., and Burkhalter, A. (2013). Distinct balance of excitation and inhibition in an interareal feedforward and feedback circuit of mouse visual cortex. *J. Neurosci.* 33, 17373–17384. doi: 10.1523/JNEUROSCI.2515-13.2013

Yang, G. R., Murray, J. D., and Wang, X.-J. (2016). A dendritic disinhibitory circuit mechanism for pathway-specific gating. *Nat. Commun.* 7:12815. doi: 10.1038/ncomms12815

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling Editor declared a past collaboration with one of the authors FW.

Copyright © 2019 Bensaid, Modolo, Merlet, Wendling and Benquet. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.