Adaptation Shapes Local Cortical Reactivity: From Bifurcation Diagram and Simulations to Human Physiological and Pathological Responses

Abstract Human studies employing intracerebral and transcranial perturbations suggest that the input-output properties of cortical circuits are dramatically affected during sleep in healthy subjects as well as in awake patients with multifocal and focal brain injury. In all these conditions, cortical circuits react to direct stimulation with an initial activation followed by suppression of activity (Off-period) that disrupts the build-up of sustained causal interactions typically observed in healthy wakefulness. The transition to this stereotypical response has important clinical implications, being associated with loss of consciousness or loss of functions. Here, we provide a mechanistic explanation of these findings by means of simulations of a cortical-like module endowed with activity-dependent adaptation and mean-field theory. First, we show that fundamental aspects of the local responses elicited in humans by direct cortical stimulation can be replicated by systematically varying the relationships between adaptation strength and excitation level in the network. Then, we reveal a region in the adaptation-excitation parameter space of crucial relevance for both physiological and pathologic conditions, where spontaneous activity and responses to perturbation diverge in their ability to reveal Off-periods. Finally, we substantiate through simulations of connected cortical-like modules the role of adaptation mechanisms in preventing cortical neurons from engaging in reciprocal causal interactions, as suggested by empirical studies. These modeling results provide a general theoretical framework and a mechanistic interpretation for a body of neurophysiological measurements that bears critical relevance for physiological states as well as for the assessment and rehabilitation of brain-injured patients.


Introduction
Studies combining transcranial magnetic stimulation (TMS) with electroencephalographic (EEG) recordings have demonstrated a dramatic impairment of the capacity of cortical networks to engage in large-scale interactions during sleep and anesthesia (Massimini et al., 2005;Ferrarelli et al., 2010;Sarasso et al., 2015): while during wakefulness the initial cortical activation triggers a chain of recurrent waves of activity and a long-range, complex pattern of interactions (Casali et al., 2013;Comolatti et al., 2019), during both non-rapid eye movement (NREM) sleep and anesthesia such distributed, rich spatiotemporal activity is lost and replaced by a response that is simple and stereotypical. Intracranial explorations employing singlepulse intracortical electrical stimulations and local field potential recordings in humans (Pigorini et al., 2015;Usami et al., 2015) showed that this stereotypical response consists of an initial activation rapidly followed by silencing of neuronal firing (Off-period), as assessed by a significant suppression of high-frequency activity (Mukovski et al., 2007;Cash et al., 2009). This tendency for cortical circuits to fall into an Off-period after initial activation during sleep and anesthesia, also known as cortical bistability (Tononi and Massimini, 2008;Nir et al., 2011), is in a key position to prevent the emergence of the complex interactions observed during wakefulness. Indeed, Off-periods not only do temporarily interrupt but also disrupt neuronal activity: as demonstrated by phase locking analysis, when neuronal activity resumes after each Off-period, it does so in a stochastic manner retaining no causal relationship with the initial input.
Crucially, this dramatic change in the input-output properties of cortical circuits is not only typical of physiological sleep and anesthesia but can also occur during wakefulness in pathologic conditions. Accumulating evidence shows that cortical circuits react to perturbations with an Off-period also in awake unresponsive wakefulness syndrome patients (UWS) (Rosanova et al., 2018), previously known as vegetative state patients, as well as in the perilesional area surrounding focal cortical lesions in awake stroke patients (Sarasso et al., 2020;Tscherpel et al., 2020) with profound clinical implications. When cortical bistability involves most of the cortex, such as in the UWS, large-scale brain interaction collapse leading to loss of consciousness (Rosanova et al., 2018). When cortical bistability is local, such as that found in the perilesional areas of stroke patients, it leads to regional circuit impairment and selective functional deficits (Sarasso et al., 2020;Tscherpel et al., 2020).
Given the clinical relevance of these findings, it is crucial to understand the nature of evoked Off-periods and their impact on cortical responsiveness within a mechanistic framework. An interesting hypothesis is that adaptation mechanisms play an important role in generating the Offperiods observed after cortical stimulation. Activity-dependent adaptation accounts for local fatigue mechanisms, i.e., self-inhibition, which lowers the firing rate of neuronal populations. At the level of cortical circuits, adaptation can be induced by several microscopic biophysical mechanisms involving calcium-dependent potassium channels, short-term synaptic depression (Gerstner et al., 2014) and/or GABAergic synaptic transmission . As suggested by mean-field theories, activity-dependent adaptation is involved in the generation of the rhythmic alternation between On-periods (highactivity up state) and Off-periods (low-activity down state) spontaneously occurring during slow-wave sleep (Latham et al., 2000;Gigante et al., 2007).
Here, we ask whether activity-dependent adaptation can reproduce and explain the fundamental features of the alteration of cortical responsiveness observed empirically during physiological sleep and during wakefulness in pathologic conditions. We address this question within a general formal framework resting on bifurcation analysis and in-silico simulations of cortical modules endowed with activity-dependent adaptation (Mattia and Sanchez-Vives, 2012). We show how such a model can provide a parsimonious explanation for the three key features characterizing the alteration of cortical reactivity observed in humans. First, using bifurcation analysis of the model we classify its dynamical regimes and reproduce parametrically both spontaneous and stimulus-evoked dynamics empirically observed during sleep and pathologic wakefulness. Second, using the same framework, we explain the empirical finding that direct cortical stimulation is more effective than the observation of ongoing dynamics in revealing Off-periods in both physiological and pathologic conditions. Third, by considering two linked modules we show how changing the adaptation level alone results in a break-off of reciprocal corticocortical interactions as observed in vivo.

Spiking neuron network
The neuronal network model is adapted from (Torao-Angosto et al., 2021), which has been proven to quantitatively reproduce the statistical features of the spontaneous up-down slow oscillations recorded in Layer 5 of the visual cortex of sleeping and anesthetized rats. Briefly, the network is composed of 6300 excitatory (E) and 2580 inhibitory leaky integrate-and-fire (LIF) neurons. The membrane potential V i of the i-th neuron evolves as emitting its k-th spike at time t ik if V i t ik ð Þ ! V thr for the first time from the emission of its previous action potential. Following the spike emission, V i ¼ V res for a refractory absolute period t 0 [2 ms (1 ms) for E (I) neurons] before restarting its time evolution (Eq. 1). Reset potential V res ¼ 15 mV and the emission threshold V thr ¼ 20 mV are the same for both neuron types, while decay constant t is 20 and 10 ms for E and I neurons, respectively. The synaptic cur- k Þ results from the spiking activity of the presynaptic neurons mediated by recurrent synaptic efficacies J ij randomly selected to be different from zero with connection probability c EE ; c EI ; c IE ; c II ¼ f0:6; 5; 0:2; 1:7g%. J ij takes positive or negative values according to the type of the j-th presynaptic neuron, and it is randomly extracted from a truncated Gaussian distribution with mean J EE ; J EI ; J IE ; J II ¼ 1:9; À1:1; 2:2; À1:1 f g mV and a relative standard deviation of 25%. Presynaptic spikes are delivered with an axonal delay d ij randomly sampled from an exponential distribution with an average of 22.6 and 5.7 ms for E and I presynaptic neurons, respectively, modeling noninstantaneous synaptic transmission . External neurons contribute to I i t ð Þ at baseline as a Poissonian spike train made of C ext independent sources each with an average firing rate of ext ¼ 0:25 Hz. Here, for E neurons, C ext range from 3200 and 3350 in the bifurcation diagram of Figure 1, while C ext ¼ 733 for I neurons. The external stimulations consist of increasing the frequency ext by a factor 4 (i.e., the stimulation intensity) for a short period of 2 ms. The external k-th spike occurring at t ext;k affects I i t ð Þ with efficacy J ext ¼ 0.48 and 2.2 mV for E and I neurons, respectively.
The second equation in Equation 1 provides the description of the activity-dependent fatigue mechanism modeling the spike-frequency adaptation (Benda and Herz, 2003;Gigante et al., 2007;Mattia and Sanchez-Vives, 2012) that affects the excitatory neurons. Specifically, the dynamics of the variable a i t ð Þ is defined by the spiking activity of the excitatory neurons: at each spike, the adaptation a i t ð Þ undergoes a unitary increase and relaxes exponentially until the next spike is emitted. Thus, the higher the excitatory population firing rate, the higher the variable a i t ð Þ. However, how much the a i t ð Þ dynamic modulates the excitatory population firing rate is also determined by the adaptation level g a ranging from 20 to 100 mV/s and the relaxation time constant for the adaptation level set at t a ¼ 150 ms. In general, the higher the a i t ð Þ because of the spikes emitted by the i-th neuron and adaptation level g a , the higher the excitatory population self-inhibition (eventually reducing its firing rate). This fatigue mechanism mimics the extracellular ionic concentrations (calcium and/or sodium) that drive a hyperpolarizing potassium current. The choice to model only the activity-dependent adaptation associated with the calciumdependent potassium currents described above does not affect the generality of the result we presented. Indeed, under the mean-field approximation detailed in the next section, the nonlinear dynamics of the network firing rate governed by other forms of activity adaptation are qualitatively similar. More specifically, the phase-plane analysis of spiking neuron networks incorporating adaptation mechanisms like short-term synaptic depression (Holcman and Tsodyks, 2006) and slow GABAergic inhibition (Parga and Abbott, 2007) are known to display similar nullclines and equilibrium states.
We also performed simulations of a network made of two identical modules. Each module is made by E-neurons and I-neurons described by Equation 1 with all the parameters set as in the single module, with the exception of the relative standard deviation for synaptic efficacies. Since the relative standard deviation for synaptic efficacies is set to zero in the two-module case, synaptic efficacies are J EE ; J EI ; J IE ; J II ¼ 1:9; À1:1; 2:2; À1:1 f g mV. Furthermore, the two modules interact with each other through connections established between excitatory neurons only, with synaptic efficacy J 12 ¼ J 21 set to 1:18 mV, where 1 and 2 identify the two modules. The probability of connection c 12 ; c 21 for each excitatory neuron is 0.1%, and average axonal delays are set to d 12 ¼ 55 ms for connections from the first to the second module, and d 21 ¼ 50 ms from the second to the first module. Furthermore, we considered the modules in a low adaptation regime, where both the modules are characterized by g a ¼ 48 mV/s, and high adaptation regime, where g a ¼ 78 mV/s.
Simulations of the spiking neuron network were performed with an open-source program relying on an eventdriven numerical integration described in (Mattia and Del Giudice, 2000). To run this software with the appropriate parameters, we used a MATLAB code available at https:// github.com/annacatt/Adaptation_and_cortical_responses and available as Extended Data 1.

Population rate model
Under mean-field approximations requiring many presynaptic contacts c ab N and limited firing rates a (Amit and Tsodyks, 1991;Amit and Brunel, 1997), where a and b identify different interacting populations, the evolution of the instantaneous firing rate a ðtÞ of an infinite-size network is described by t _ a ¼ À a 1 Uð a Þ (Treves, 1993;Brunel and Hakim, 1999;Mattia and Del Giudice, 2002). Here, the current-to-rate gain function U is the Siegert-Ricciardi one (Siegert, 1951;Capocelli and Ricciardi, 1971): ; where m and s 2 are, respectively, the infinitesimal mean and variance of the input current (Amit and Tsodyks, 1991;Amit and Brunel, 1997) whose form is shown below,   Figure 1. Dynamical regimes of the spiking neuron network accompanied by the bifurcation analysis of the population rate model. A, Dynamical regimes of the spiking neuron network in the excitation-adaptation parameter space, namely (C ext ; g a ), classified according to its spontaneous activity (for details, see Materials and Methods, Analysis of the dynamical regimes of the spiking neuron network) superimposed to the codimension-two bifurcation diagram of the population rate model. The codimension-one bifurcations encompass the subcritical Andronov-Hopf bifurcation curves (dashed magenta lines), the supercritical Andronov-Hopf bifurcation curves (solid magenta line), the and erf is the error function. The first-order dynamical equation t _ a ¼ À a 1 Uð a Þ can be seen as the simplest relaxation dynamics leading to the asymptotic firing rates given by Uð a Þ. Here, 1=Uð a Þ is the mean of the interspike intervals (ISI) expected for the leaky integrate-and-fire neuron driven by a stochastic current with infinitesimal mean m and variance s 2 . The ISI in this model is the time when the membrane potential VðtÞ crosses the emission threshold V thr for the first time starting from the reset potential V res at time t ¼ 0. This interval also includes the absolute refractory period (t 0 ). The above dynamics for a ðtÞ is obtained by resorting to a spectral expansion of the Fokker-Planck equation and determining the evolution of the probability density of the membrane potential. The first-order ordinary differential equation results from taking into account only the slowest eigenmode, i.e., the eigenfunctions of the Fokker-Planck operator with the highest eigenvalue (Treves, 1993;Mattia and Del Giudice, 2002). This first-order ODE recovers the phenomenological dynamics for the firing rate introduced by Wilson and Cowan (1972). Since the cortical module we are considering includes two interacting populations of excitatory and inhibitory neurons, the mean-field dynamics is: where a 2 fE; Ig; U a ðm a ; s a Þ is the above gain function computed with the single-neuron parameters of type a and t a is the decay constant of the LIF membrane potential. The infinitesimal mean and variance of the input current are defined as: with a; b 2 E; I f g and D ab are the relative standard deviations of the synaptic efficacies J ij for the different neuronal types. This mean-field approximation effectively describes the collective dynamics of networks of LIF neurons with spike-frequency adaptation even far from the equilibrium (Gigante et al., 2007;Mattia and Sanchez-Vives, 2012). Here, only excitatory neurons have adaptation: g E g a and g I ¼ 0.
Finally, because of the breaking of the diffusion approximation valid only in the limit J ab ! 0, we shifted horizontally the critical point of this bifurcation diagram by an appropriate amount (DC E;ext ¼ À310) of excitation level. Indeed, in the diffusion limit the currents received by the neurons are continuous stochastic processes determined by the continuous barrage of the synaptic input arriving at high rates and inducing small jumps in the membrane potential (Tuckwell, 1988). However, in the model network we simulated, synaptic efficacies are not negligible and the membrane potentials have a jump-like evolution in time. In this shot-noise regime, firing rates are lower compared with the ones predicted under the diffusion approximation (Richardson and Swarbrick, 2010), and this explains why we need to incorporate this effect by reducing the excitation level in our mean-field theory, eventually leading to a rightward shift of the theoretical bifurcation diagram. Furthermore, we also shifted vertically the critical points of the bifurcation diagram by an appropriate amount (Dg a ¼ À7 mV/s) of adaptation. Besides the breaking of the diffusion approximation, the shifts DC E;ext and Dg a are expected to be theoretically justified by the fact that the simulated networks are composed of a finite number of spiking neurons challenging the infinite-size limit on which extended mean-field approximation relies on (Mattia and Del Giudice, 2002).

Simulated data and data preprocessing
Simulations consist of N trials (N ¼ 50 for the single cortical module and N ¼ 250 for the case of two interacting modules). The external stimulus in each trial was delivered at time t ¼ 0, with a prestimulus interval lasting randomly between 5000 and 5300 ms and poststimulus interval lasting 5000 ms. Stimulation artifact was reduced continued saddle-node bifurcation curve (blue curves), and the limit point of cycles (green lines, also known as saddle nodes of periodic orbits) that turn the periodic orbits originating from Andronov-Hopf stable. The codimension-two bifurcation Bogdanov-Takens (BT) is the contact point between the saddle-node bifurcation curve and Andronov-Hopf bifurcation curves. The saddle-node curve shows a cusp bifurcation (CP). The areas among the curves have been color-coded according to the regime shown by the spontaneous dynamics of the spiking neuron network. The five different dynamical regimes are characterized by HAS (high-asynchronous state; blue area), HAS/SO (high-asynchronous state with incursions of slow oscillations, namely with incursions of Off-periods because of finite-size effects; orange area), SO (slow oscillations; purple area), LAS/SO (low-asynchronous state with incursions of slow oscillations; green area), and, finally, LAS (low-asynchronous state; green area). A sixth region (gray area), not of interest in the present work, is characterized by HAS and LAS alternating at an irregular pace. B, Spontaneous and stimulus-evoked signals with the same color-coding used for the dynamical regimes of the spiking neuron network (panel A) for a fixed excitation (C ext ¼ 3297:5Þ and several adaptation levels [Case 1: 30; Case 2: 40; Case 3: 45; Case 4: 57.5; Case 5: 90 (mV/s)]. Left column, Time series encompassing both spontaneous activity (up to 6000 time steps) and stimulus-evoked activity (remaining interval). Black triangles indicate the occurrence of the stimulation. Central column, Spontaneous activity as a function of fatigue (g a Á a mV/s). Right column, Three superimposed orbits because of perturbations as a function of the fatigue. In both the central and right columns, the spontaneous and stimulus-evoked dynamics are superimposed to the nullclines of the population rate model. The intersection of the nullclines represents the equilibrium point of the corresponding population rate model. The equilibria of the spiking neuron network and the population rate model show a close albeit not complete overlap. The dynamics of the spiking neuron network for a fixed excitation and different adaptation levels are provided in Extended Data Figure 1-1. C, Probability to evoke Off-periods in wakefulness (Case 1) and deep sleep (Case 3) as a function of several stimulation intensities spanning from 2 to 6 ([a.u.]). The probability to evoke Off-periods is computed as the number of trials, lasting 2 s each, showing at least one Off-period divided by the total number of trials (50). by applying a Tukey-windowed median filtering, as in (Chang et al., 2012), between À5 and 5 ms (Pigorini et al., 2015).

Bifurcation analysis of the population rate model
A detailed analysis of local bifurcations of Equation 2 was performed through MatCont (Dhooge et al., 2003), a MATLAB Toolbox that allows computing curves of equilibria and bifurcation points through a prediction-correction continuation algorithm, as described previously (Kuznetsov, 2004).

Analysis of the dynamical regimes of the spiking neuron network
The dynamical regimes of the spiking neuron network (Eq. 1) were determined according to the frequency of On-periods and Off-periods (up and down states, respectively) detected across 50 trials of spontaneous activity lasting 2 s each. The detection of the On-periods and Offperiods was conducted by looking at the firing rate in each trial. If the firing rate exceeded 20 Hz, the trial contained at least one On-period. On the contrary, if the firing rate fell under 5 Hz, the trial contained at least one Off-period. Once we collected the probability of On-periods and Off-periods across 50 trials for each value of excitation and adaptation level, we classified five dynamical regimes of interest. Following nomenclature previously published (Gigante et al., 2007;Mattia and Sanchez-Vives, 2012), high-asynchronous state (HAS) is characterized by a probability of On-periods equal to 1 (i.e., sustained high firing rate in every trial) and absence of Off-periods. Low-asynchronous state (LAS) presents a probability of Off-periods equal to 1 (i.e., persistent low firing rates in every trial) and absent On-periods. Slow oscillations state (SO) is characterized by both probabilities of On-periods and Off-periods equal to 1 since each trial sees On-periods and Off-periods alternating at a regular pace. Furthermore, on the border between LAS and SO, there is a region where the probability of Off-periods is equal to 1 with incursions of On-periods states (nonzero probability of On-periods , 1). Vice versa, the region between HAS and SO is characterized by On-periods in all the trials with incursions of Off-periods (non-zero probability of Off-periods ,1). Finally, a sixth region, which is not the focus of this work, encompasses both HAS and LAS with non-zero probability of On-periods and Off-periods lower than 1. This is the case of heterogeneous trials, some containing only On-periods, some only Off-periods, and some exhibiting both On-periods and Off-periods alternating at an irregular pace.

Simulated data analysis in the power and phase domain
Data analysis was performed using MATLAB R2020b (The MathWorks Inc.). We used the timef function implemented in EEGLAB (Delorme and Makeig, 2004) to detect transient event-related spectral perturbation (ERSP) and event-related phase-locking, i.e., inter-trial coherence (ITC) events, in the simulated data. More in detail, single trials were time-frequency decomposed between 5 Hz and 45 Hz using Wavelet transform (Morlet, window span: 3.5 cycles). The resulting ERSPs and ITCs were averaged across trials and normalized by subtracting the mean spectral activity of the prestimulus baseline from À1000 to À400 ms (from À5000 to À400 ms in the two-module network). To detect statistically significant activation with respect to the prestimulus baseline, we applied bootstrap statistics to both ERSP and ITC with significance level a ¼ 0.005 and 1000 permutations. Nonsignificant bins were zeroed out. Furthermore, to discard the spurious intertrial coherent events associated with the absence of firing rate activity underlying the ERSP suppression, we retained only the ITC points associated with a significant increase in the ERSP power with respect to the prestimulus baseline. Finally, we averaged the ITC values across frequencies for each time point. The resulting average ITC provides an indication of the duration of the deterministic effect in a wide frequency window of a given external input. The code for the analysis is available at https://github.com/annacatt/Adaptation_ and_cortical_responses and as Extended Data 1.
Methods and results related to human data are presented in published articles (Massimini et al., 2007;Pigorini et al., 2015;Rosanova et al., 2018;Sarasso et al., 2020) cited throughout this work and derive from analysis in the power and phase domains similar to the ones used for the simulated data presented here.

Code accessibility
The code created for this paper is freely available online at https://github.com/annacatt/Adaptation_and_cortical_ responses and is available as Extended Data 1.

Data availability
All the simulated data that support the findings of this study are available online at https://figshare.com/articles/dataset/ Simulated_data/21112603.

Results
We employ a spiking network model to investigate to what extent varying the strength of adaptation and excitation can explain fundamental dynamics of cortical responsiveness empirically observed in healthy humans and brain-injured patients. Specifically, we aim to replicate and provide a mechanistic explanation of experimental data encompassing TMS-EEG recordings in healthy subjects during NREM sleep stage 2 (N2) (Massimini et al., 2007), in UWS patients with severe brain injury (Rosanova et al., 2018), in stroke patients (Sarasso et al., 2020;Tscherpel et al., 2020), as well as intracranially-evoked potentials in wakefulness/sleep (Pigorini et al., 2015).

Bifurcation analysis and characterization of the dynamical regimes
Activity-dependent adaptation and excitation level shape the ongoing and stimulus-evoked activity of the simulated cortical module (Eq. 1). To fully characterize the dynamical regimes the model accounts for, we first performed a bifurcation analysis of its continuous counterpart (Eq. 2) exploiting a numerical continuation technique (Dhooge et al., 2003;Kuznetsov, 2004). This analysis expands on previous studies (Gigante et al., 2007;Mattia and Sanchez-Vives, 2012), providing a complete characterization of the bifurcations of the mean-field dynamics of the network. Taking into account only two key parameters, i.e., the level of excitation (C ext ) and adaptation (g a ; see Materials and Methods), the bifurcation diagram in Figure 1A shows critical points such as saddle nodes (blue traces) and subcritical Hopf bifurcations (dashed magenta traces). Subcritical Hopf bifurcations lie near the saddle nodes of periodic orbits (green traces), also known as limit points of cycles, where stable limit cycles arise. Thus, in the region limited by these two curves, stable slow oscillations are generated, corresponding to the cortical bistability observed in neurophysiological recordings performed in anesthetized rats (Tort-Colet et al., 2021). On the contrary, the region enclosed by the limit point curve (blue trace) is characterized by two stable and one unstable equilibrium solutions. The coexisting stable equilibria are characterized by high and low asynchronous firing states. In order to disambiguate, it is worth noting that this regime, which is not the focus of the present work, is also indicated with the term "bistable" as classically understood in the field of dynamical systems (Gerstner et al., 2014).
To check the effectiveness of this bifurcation diagram, we characterized the spontaneous dynamical regimes of equivalent networks of spiking (adaptive leaky integrateand-fire) neurons in Equation 1. The leftmost part of the traces in Figure 1B, left column, shows 6 s of spontaneous activity for a given excitation value and five different adaptation levels. By taking into account the probability of Off-periods and On-periods during spontaneous activity (for details, see Materials and Methods), we identified different regions. The light blue region displays HAS (high-asynchronous state), i.e., relatively high firing rate and absence of Off-periods, which is typical of wakefulness. The light red region is characterized by LAS (lowasynchronous state), i.e., low asynchronous firing rate distinctive of Off-periods and absence of On-periods, a pattern resembling burst-suppression. The purple region features a regular alternation between HAS and LAS at a regular pace generating slow oscillations (SOs), as usually observed under anesthesia. The transition between these regimes is not as sudden as in the population rate model while varying one or both parameters. Indeed, finite-size effects make a state dominant, with incursions of the state of the nearby region (Gigante et al., 2007;Mattia and Sanchez-Vives, 2012). This is the case of the orange and green regions characterized by HAS and LAS, respectively, with incursions of SO. Specifically, HAS sparsely interspacing SO well describes deep sleep recordings.

Perturbations reveal cortical bistability beyond spontaneous activity
Crucially, we further characterized the dynamical regimes of the spiking neuron model (Eq. 1) by exploring its responses to perturbations, which consist of an external current lasting 2 ms with intensity 4 ([a.u.]) injected into the excitatory neurons resembling the effects of direct magnetic and electrical stimulation. The role of the perturbation is to transiently alter the spontaneous population firing rate, potentially unveiling interesting activity-dependent properties of the system. The rightmost part of the traces in Figure 1B, left column (presenting Cases 1-5) shows 6 s of the stimulus-evoked activity for a given excitation value and five different adaptation levels. In the middle and right columns, the firing rate traces have been represented as a function of the average fatigue determining the adaptation of spike rates for both the spontaneous and the stimulus-evoked activity. The external stimulus applied to the network (defined by duration and intensity) increases the magnitude of the external Poisson spike trains. We fixed the duration to 2 ms and checked for the presence of Off-periods in the evoked activity for several stimulus intensities. As shown in Figure  1C, the stimulus intensity was finally set to never (always) evoke Off-periods in Case 1 (Case 3), which describes wakefulness (deep sleep).
In Case 1, the stimulus-evoked activity and spontaneous activity do not show relevant dissociations, as both are characterized by the absence of Off-periods, as empirically observed during wakefulness in human data. Higher adaptation levels, corresponding to Cases 3-5, are all characterized by the presence of Off-periods in both the spontaneous and evoked activity. The overall time the system spends in the Off-periods increases with the adaptation level: (1) Case 3 shows Off-periods interspersed with up states, thus resembling NREM sleep stage 3 (N3); (2) Case 4 displays regular alternation between On-periods and Off-periods, thus mirroring the effect of anesthesia on brain activity; (3) Case 5 reminds of burst suppression characterized by prolonged Off-periods. Crucially, Case 2 in Figure 1B reveals a clear-cut dissociation between the dynamical features observed in the spontaneous activity and the stimulus-evoked response. Here, while Off-periods are extremely rare (one every ;16 s on average) in spontaneous activity, they can be reliably triggered by 76% of the stimulations. Hence, the model shows the existence of a particular regime whereby an interventional approach manifests a substantial sensitivity in revealing the underlying state of the system and its position in the bifurcation graph compared with an observational approach. In other words, it provides a crucial tool to detect cortical bistability with a few stimulations instead of relying on long spontaneous recordings. Such dissociation between observable dynamics and responses to perturbations is important because it reproduces a key feature reported by empirical works. This is illustrated in Figure 2 where the results of simulations are directly compared with the results of TMS-EEG experiments both in the time domain and frequency domain to detect the presence of cortical bistability (slow waves and Off-periods). This analysis shows a fundamental qualitative correspondence between Case 2 of the model ( Fig. 2A) and spontaneous as well as evoked patterns of electrophysiological activity found in humans across different conditions encompassing N2 sleep (Massimini et al., 2007; Fig. 2B), UWS patients (Rosanova et al., 2018;Fig. 2C), and perilesional area of stroke patients (Sarasso et al., 2020;Fig. 2D). In all these cases, TMS-evoked slow waves and Off-periods, as detected by a significant high-frequency (.20 Hz) power suppression, are reliably revealed by cortical perturbations, whereas they are rare in the spontaneous prestimulus activity.
As a note, Figure 2 shows a qualitative agreement between Case 2 in the model and the experimental data. By choosing the model in Equation 1, we aimed to unveil the mechanisms behind cortical bistability through an extensive exploration of the model phenomenology. Given that the model describes the firing rates of a local neuronal population, which differ from the experimental data in both the kind of signal and the spatial scale involved, a quantitative comparison between simulated and experimental data were beyond the aim of the present work.
Considering its prevalence in real data across different physiological and pathologic conditions, we asked whether such dissociation between spontaneous and stimulus-evoked activity is a peculiar feature of the specific level of excitation, or whether it can be generalized to a large interval of excitation levels. Figure 3 shows the difference in the frequency of Off-periods detected by analyzing the spontaneous and evoked activity for each excitation-adaptation pair. We found an area surrounding the low Hopf bifurcation curve, extending to ;10% of the analyzed excitation-adaptation plane, characterized by a pronounced difference between evoked and spontaneous Off-periods. More specifically, despite a marked difference being present for all the excitation  Figure 1, and in TMS-EEG data (during physiological N2 sleep, UWS patients, and in the perilesional area of stroke patients). First column, Time series with a few stimulations. Second column, Trials without Off-periods in the baseline aligned to the stimulation onset. Third column, Event-related spectral perturbation (ERSP). Blue color indicates a significant reduction compared with the baseline, while red indicates a significant increase. Black triangles and red lines indicate the occurrence of the stimulation. A, Simulated firing rate data. B, Re-edited from Massimini et al. (2007). EEG signal recorded from a channel (Cz) located under the stimulator during one TMS-ON block over a background of spontaneous NREM sleep (single-subject data). The TMS-ON block consisted of 40 stimuli at 0.8 Hz. The red dashed section shows the slow waves triggered at the beginning of the block. C, Re-edited from (Rosanova et al., 2018). EEG recordings during TMS stimulations in a UWS patient. EEG activity of one representative electrode, Cz, while TMS was delivered with an interstimulus interval randomly jittering between 5000 and 5300 ms. Middle panel, Trials aligned to the stimulation displaying baseline activity without spontaneous slow waves. Panel D is derived from published data presented in Table 1 and Figure 2 in Sarasso et al. (2020). EEG recordings during perilesional TMS stimulations in a middle cerebral artery ischemia patient (patient 4). EEG activity recorded from a channel (Fc3) located over the perilesional area while TMS was delivered with an interstimulus interval randomly jittering between 2000 and 2300 ms. ERSPs in B-D were computed as by Rosanova et al. (2018). levels higher than C E;ext ¼ 3250, a near-complete dissociation (;84% of evoked trials detect Off-periods where the spontaneous trials do not show any) can be shown for very high excitation levels.
Given the shallow slope of the lower branch of the Hopf bifurcation, changes in adaptation are expected to have more profound effects on cortical bistability than shifts in the excitation level. However, substantial changes in excitation levels for a given adaptation (i.e., horizontal shifts toward the left in the bifurcation diagram) still result in changes in spontaneous and stimulus-evoked dynamics similar to those reported in Figure 1. Extended Data Figure 1-1 illustrates such cases, including the interesting instance (Case 2) characterized by an intermediate excitation level, showing the dissociation between spontaneous and stimulus-evoked dynamics.

High adaptation level is associated with early breakdown of the causal interactions compared with low adaptation
Empirical studies employing perturbations in humans (Pigorini et al., 2015;Usami et al., 2015;Rosanova et al., 2018), rodents (Arena et al., 2021), and cortical slices (D'Andola et al., 2018) suggest that cortical bistability and the associated Off-periods can disrupt the emergence of sustained patterns of causal interactions among cortical neurons. We thus used the model to test whether changes in adaptation level can also affect casual interaction among groups of cortical neurons. To do that, we considered two cortical modules, instead of one as in previous sections, connected through reciprocal bidirectional connections. As in this new model composed of two interacting modules, a complete characterization of the bifurcations would require exceedingly high computational cost, we here focused on a representative scenario characterized by two extreme adaptation levels (Fig. 4A). Specifically, the low adaptation level parallels the case of wakefulness as seen in the previous section (Case 1; Fig. 1B), as opposed to the high adaptation level characterizing N3 sleep or deep anesthesia (Case 4; Fig. 1B). Under conditions of low adaptation, cortical perturbations elicited reverberant interactions between the two modules as reflected by multiple, recurrent and coherent oscillations of evoked activity. Conversely, under conditions of high adaptation, the interplay between the two simulated groups of cortical neurons was short-lasting resulting in only a few oscillations followed by a prominent Off-period, as marked by the suppression of high-frequency power, similar to that observed in sleeping and brain-injured humans.
In real-life experiments, the impact of the Off-period on the ability of cortical circuits to sustain causal interactions is assessed by quantifying the duration of the deterministic effects of the perturbation in the phase domain. A key empirical finding is that Off-periods not only disrupt deterministic interaction because of the associated suppression of power but that they also scramble the subsequent phase of the signal when neurons resume firing. Such stochastic reboot is characterized by oscillations with high power but low phase-locking or intertrial coherence following the Off-period. We thus performed the same analysis in the model by computing the intertrial coherence (ITC) averaged across frequencies between 5 and 30 Hz.  . High adaptation level results in early breakdown of causal interactions compared with low adaptation level. A, Simulated data of two coupled modules for low (g a ¼ 48 mV/s) and high (g a ¼ 78 mV/s) adaptation levels. First row, Firing rate activity averaged over 250 trials. A dashed vertical line (at t ¼ 0) marks the occurrence of the external stimulus injected into module 1. Second row, Event-related spectral perturbation (ERSP) for module 2. Significance for bootstrap statistics is set at a , 0.005: absence of significant activation is colored in green, significant increases of power compared with baseline are represented in red, while As shown in Figure 4A, the Off-period occurring under conditions of high adaptation was followed by a broadband resumption of power which was however nondeterministic as indicated by a concomitant absence of significant ITC. The primary sources of such stochasticity are (1) the variable level of fatigue accumulated by the network at the stimulation time and (2) the activity fluctuations because of the finite-size of the neuronal population determining the escape time from the metastable inactive state (Off-period).
Notably, in the case of low adaptation, significant deterministic interactions were not curtailed by the Off-period and were longer-lasting. As illustrated in Figure 4B, these model predictions resulting from the manipulation of the adaptation levels match the power and phase modulations observed in experiments in which TMS or intracranial electrical stimulations are applied in different global brain states, such as during wakefulness, deep NREM sleep and in UWS patients.

Discussion
By exploiting a spiking network model endowed with activity-dependent adaptation, the present work provides a principled mechanistic explanation for empirical results revealing common alterations of cortical responsiveness in sleeping healthy subjects, and in awake patients with multifocal and focal brain injuries. Specifically, (1) we link changes in cortical reactivity to the presence of underlying activity-dependent adaptation mechanisms, above and beyond spontaneous dynamics, (2) we explain the role of perturbations in revealing adaptation mechanisms, and (3) we expound on the impact of adaptation in disrupting causal interactions among cortical modules.
The first implication of this work is that systematically varying the relationships between excitation and adaptation levels within a formal bifurcation analysis reproduces basic patterns of cortical reactivity to direct perturbations observed in humans. As shown in Figure 1A,B, increasing levels of adaptation and/or decreasing the excitation engenders a condition whereby cortical circuits react to a direct perturbation with an initial activation followed by an Off-period, yielding responses similar to those found during sleep and in awake brain-injured patients. This finding offers important elements for interpreting macroscale empirical data in terms of neuronal mechanisms.
In this framework, the Off-periods observed after cortical stimulation in sleeping humans can be ascribed to adaptation mechanisms because of slow negative (i.e., inhibitory) feedback produced by calcium-gated and sodium-gated potassium currents (Sanchez-Vives and McCormick, 2000;Sanchez-Vives et al., 2010), which are enhanced by decreased levels of neuromodulation from brainstem activating systems (McCormick, 1992). Notably, the present simulation results closely match the changes in reactivity to electrical stimulation observed in cortical slices when modulating potassium currents by applying carbachol and norepinephrine (D'Andola et al., 2018). Recent studies have also pointed to a putative role of active inhibition in conditioning the onset and duration of the Off-period (Funk et al., 2017;Zucca et al., 2017). As adaptation in the present model encompasses self-inhibition, the role of local inhibitory neurons in shaping the changes in cortical responses observed during sleep is not incompatible with the present results.
Interpreting the empirical results obtained in brain-injured and stroke patients within the theoretical framework of the bifurcation diagram discloses a more complex and interesting landscape. Indeed, brain lesions can cause global or local alterations of adaptation and excitation through different mechanisms, even in the context of preserved arousal, as assessed by eyes opening. For example, in some UWS patients, lesions, compressions, or displacements of brainstem activating systems as well as a critical load of damage to ascending fibers in subcortical white matter may enhance potassium currents (Steriade et al., 1993;Edlow et al., 2012), which corresponds to an upward shift in the bifurcation diagram. In other cases, cortical injuries and white matter lesions can engender a state of corticocortical disfacilitation by affecting the excitation term (Takahashi et al., 1981;Lemieux et al., 2014;Rocchi et al., 2022). In this case, the resulting loss of lateral and long-range excitatory input would act by producing shifts on the horizontal axis toward the left, to a point where the probability of evoking an Off-period in an awake patient becomes higher. An interesting implication of the diagram, especially considering the slope of the lower branch of the bifurcation (Figs. 1A,  3), is that changes in adaptation are expected to have more dramatic effects on cortical bistability as compared with changes in the excitation level. This prediction finds empirical confirmation in cortical slices, an extreme model of cortical injury, whereby reducing adaptation by application of carbachol and norepinephrine is more effective in recovering wakelike responses than increasing excitation by application of kainate (D'Andola et al., 2018). In real-life conditions, however, the two mechanisms (adaptation and excitation) are not mutually exclusive and may have different relative weights depending on the type/combination of injuries. In fact, they may both concur in bringing residual cortical circuits into a state in which cortical sleep-like Off-periods are generated during wakefulness, as observed in stroke (Sarasso et al., 2020) and UWS (Rosanova et al., 2018) patients. continued significant power decreases are colored in blue. Third row, Eight firing rate traces of module 2. Fourth row, Averaged intertrial coherence (ITC) for frequencies between 5 and 30 Hz. B, Data recorded in human. The first column shows results related to stereo-EEG with SPES during wakefulness (W-red) and NREM sleep (NREM-blue) re-edited from Pigorini et al. (2015). From top to bottom, average responses of a representative contact, event-related spectral perturbation (ERSP), and phase-locking factor (PLF) for frequencies higher than 8 Hz (details can be found in Pigorini et al., 2015). Second and third columns show results related to TMS-EEG re-edited from Rosanova et al. (2018), encompassing healthy wakefulness (W-red), healthy NREM sleep (NREM-blue), and UWS patients during wakefulness (gray W). As for stereo-EEG-SPES data, from top to bottom, average responses of the channel under the stimulator, event-related spectral perturbation (ERSP), and phase-locking factor (PLF) for frequencies higher than 8 Hz (details can be found in Rosanova et al., 2018).
The second implication of this work lies in the dissociation between observational and perturbational approaches in their ability to unveil the presence of cortical bistability and Off-periods, and bears relevance for the assessment of the aftermath of brain injury. Such dissociation is represented by the area in the excitation-adaptation bifurcation diagram identified in Figure 3, in which the responses to perturbations more reliably show Off-periods as compared with spontaneous activity. This portion of the diagram encompasses a relatively small space of the overall dynamics (ranging from wake-like activity to patterns resembling burst suppression) but is extremely relevant in real-life conditions. Indeed, cortical perturbations, either electrical or magnetic, are often capable of evoking clear-cut Off-periods, which are not otherwise present in spontaneous activity not only during N2 sleep but also in many stroke and traumatic brain-injury patients (Rosanova et al., 2018;Sarasso et al., 2020). This is a critical working point along with the sleep-wake transition, in which cortical dynamics are intrinsically unstable as two global activity modes compete (Parrino et al., 2012;Tort-Colet et al., 2021). From a therapeutic perspective, this offers a window of opportunity as small physical or pharmacological perturbation of the network parameters can induce dramatic changes in the brain global behavior (Deco et al., 2019). Thus, relating empirical findings to this region of the bifurcation diagram is important for two reasons. First, it suggests that a significant number of patients may lay in a state where the inputoutput properties of cortical circuits are altered because of a critical, albeit potentially reversible, shift in adaptation (and/or excitation). Second, it shows that, because of its inherent activity-dependent nature, this state of affairs can be better revealed, above and beyond the observation of spontaneous dynamics, by challenging cortical circuits with direct stimulation.
The third result of this theoretical and computational work is that altered levels of adaptation can have profound effects on the capability of cortical neurons to engage in reciprocal interactions, as indicated by the complex set of effects observed when changing adaptation in two reciprocally connected cortical modules without modulating their connectivity strength. Under conditions of low adaptation, the two simulated modules engage in a long-lasting series of feed-forward and feed-back interactions leading to multiple waves of activity time-locked to the stimulus, resembling the general pattern of responsiveness found in cortical slices under carbachol and norepinephrine (D'Andola et al., 2018) and in healthy awake subjects (Massimini et al., 2005). For high adaptation, this deterministic pattern of interaction is drastically curtailed; not only Offperiods temporarily obliterate activity, but they also disrupt phase-locking to the stimulus once activity resumes. This peculiar condition, whereby high power is associated with minimal levels of phase-locking to the stimulus, is strikingly similar to the pattern found in multiscale empirical measurements ranging from cortical slices (D'Andola et al., 2018) and anesthetized rodents (Arena et al., 2021) to human intracranial and extracranial measurements during NREM sleep and after severe brain injury (Pigorini et al., 2015;Rosanova et al., 2018). Notably, in the brain of patients, Off-periods and the ensuing disruption of causal interactions are empirically associated with low values of whole-brain complexity and with loss of consciousness (Rosanova et al., 2018). Perhaps more importantly, the progressive disappearance of evoked Off-periods is associated with recovery from disorders of consciousness (Rosanova et al., 2018) and stroke (Tscherpel et al., 2020). In light of the present theoretical framework, this clinical evolution would correspond to the descending trajectory represented in the bifurcation diagram of Figure 1A, with potential implications for stratification, follow-up, and rehabilitation in the aftermath of brain injury. For example, detecting cortical bistability by perturbations in a stroke patient points to the presence of functional disruption, adding to the structural damage, and suggests that neuromodulation or pharmacological treatment should aim at reducing adaptation mechanisms and/ or strengthening local excitation until the occurrence of evoked Off-period is minimized.
Clearly, the present theoretical framework only represents a first stepping stone on which more realistic and complex models can be built, such as those incorporating the topological organization of cortical networks (Capone et al., 2019;Barbero-Castillo et al., 2021;Pazienti et al., 2022). To be comprehensive such models should also include other subcortical structures like the thalamus, whose input is known to play an important role in shaping Off-periods (van Wijngaarden et al., 2016;Zucca et al., 2019). Besides, the present work only provides a minimal account, limited to the proof of principle of two connected modules of the effects of adaptation on corticocortical interactions. Hence, a fundamental development will consist in embedding the present framework within a large-scale, connectome-based simulation, such as the "The Virtual Brain," encompassing a multitude of interacting modules (Sanz Leon et al., 2013;Kringelbach et al., 2020;Goldman et al., 2022;Schirner et al., 2022). This would offer a tool to better understand the effects of local intrusions of cortical bistability within the awake brain, such as those occurring after sleep deprivation (local sleep; Hung et al., 2013;Sarasso et al., 2014;Bernardi et al., 2015;Nir et al., 2017) as well as those occurring when adaptation and excitation are altered in a regional-specific manner by focal and multifocal structural lesions.
Several computational models have been proposed in the last years to describe in detail the effect of electrical stimulations on brain circuits (Bai et al., 2013;McIntyre and Foutz, 2013;Seo and Jun, 2017;Farokhniaee and McIntyre, 2019;Komarov et al., 2019). Despite some lack of biophysical details in our model network, we implemented a parsimonious description of the general effect of physical (magnetic or electric) cortical perturbations. Indeed, we considered both electrical and magnetic stimulations as capable of inducing intraparenchymal electric fields (Nieminen et al., 2015;Rahman et al., 2015;Laakso et al., 2018). Such local and exogenous fields are known to modulate neuronal excitability inducing a polarization of the membrane potential of pyramidal cells with orientation parallel to the field vector (Bikson et al., 2004;Radman et al., 2009). Based on this experimental and theoretical evidence, we eventually modeled TMS and electrical stimulations as a transient change in the input currents received by excitatory (i.e., pyramidal) cells in our spiking neuron network, leading in turn to a depolarization of their membrane potential. Remarkably, such a simple model allowed us to qualitatively reproduce the evoked responses shown in experimental works and to reveal that the adaptation level in the network strongly shapes those evoked responses. In the future, our network will benefit from increased biophysical details to differentiate between electrical and magnetic stimulations, and to explore the fine grained, intracolumnar events underlying activity-dependent adaptation.