The authors have declared that no competing interests exist.
Segregation and integration are two fundamental principles of brain structural and functional organization. Neuroimaging studies have shown that the brain transits between different functionally segregated and integrated states, and neuromodulatory systems have been proposed as key to facilitate these transitions. Although whole-brain computational models have reproduced this neuromodulatory effect, the role of local inhibitory circuits and their cholinergic modulation has not been studied. In this article, we consider a Jansen & Rit whole-brain model in a network interconnected using a human connectome, and study the influence of the cholinergic and noradrenergic neuromodulatory systems on the segregation/integration balance. In our model, we introduce a local inhibitory feedback as a plausible biophysical mechanism that enables the integration of whole-brain activity, and that interacts with the other neuromodulatory influences to facilitate the transition between different functional segregation/integration regimes in the brain.
Segregation of brain activity refers to the fact that some brain regions are specialized to handle particular features of external and internal stimuli. However, to produce a coherent behavioral outcome, the brain must coordinate the activity of these specialized brain areas, and this is called integration of brain activity. Based on a fixed connectome (the brain anatomical structure), the neuromodulatory systems are one of the plausible candidates to manage the transitions of brain states in short timescales. Understanding the role of neuromodulators in brain dynamics and the segregation/integration balance is relevant, in particular, as it is known that in several neuropsychiatric disorders the segregation/integration balance its impaired. Here, we used a computational model of the whole brain to study the dual effect of the cholinergic and noradrenergic neuromodulatory systems in the switching from segregated to integrated brain states. The novelty of our work is the inclusion of a homeostatic local inhibitory loop. This specific inhibition, modulated by the cholinergic system, maintains the excitation/inhibition balance while promoting integration. Our work links the local effects of cholinergic neuromodulation, with the more global influences of the structural connectivity and neuromodulatory systems. This constitutes a step forward in the understanding of the neural mechanisms behind the segregation/integration balance of brain activity.
Integration and segregation of brain activity are nowadays two well-established brain organization principles [1–4]. Functional segregation refers to the existence of specialized brain regions, allowing the local processing of information. Integration coordinates these local activities in order to produce a coherent response to complex tasks or environmental contexts [1, 2]. Both segregation and integration are required for the coherent global functioning of the brain; the balance between them constitutes a key element for cognitive flexibility, as highlighted by the theory of coordination dynamics [5, 6].
From a structural point of view, the complex functional organization of the brain is possible thanks to an anatomical connectivity that combines both integrated and segregated network characteristics, having small-world and modular properties [7]. In spite of this structural connectivity (SC) remaining fixed over short timescales, different patterns of functional connectivity (FC) can be observed during the execution of particular behavioral tasks [2]. Moreover, functional Magnetic Resonance Imaging (fMRI) neuroimaging studies show that during a resting state the FC is not static, but rather evolves over the recording time [8–10], highlighting the non-linear and non-stationary properties of the FC [11]. In a similar way, the integration and segregation of brain activity are not static over time [3, 12]. In this context, an interesting question emerges: How does the brain manage to produce dynamical transitions between different functional states from a rigid anatomical structure?.
Neuromodulatory systems tune the firing properties of neurons, providing a mechanism to change the flow of information within the brain, and allowing the transitions between different FC patterns. A recent hypothesis proposed by Shine [13] argues that neuromodulation allows the transition between integrated and segregated states, manipulating the neural gain function [14]. In that line, the cholinergic and noradrenergic systems have been proposed as candidates to influence the cognitive processing within the brain [15, 16], in spite of not being the unique neuromodulatory systems in the central nervous system which can tune the firing properties of neurons [14, 17].
The cholinergic system is involved in cognitive and attentional selectivity [16], and in the cerebral cortex the main source of acetylcholine are projections from the basal forebrain [18]. Acetylcholine increments the overall excitability [19, 20], and consequently rises population activity above noise, a mechanism referred as response gain [14]. The increase in signal-to-noise ratio, especially in brain areas that are close to each other, promotes segregation when considering the response gain by itself [13]. On the other hand, the noradrenergic system is related to the exploratory behavior [15], and the principal source of noradrenergic projections to the cerebral cortex comes from the locus coeruleus [21]. Noradrenaline increases the reponsivity (or selectivity) of neuronal populations to input-driven activity (e.g., sensory stimuli, inputs for distant brain areas relevant to a task) with respect to spontaneous activity (or the internal state of the brain) [22–24], filtering out noise [25] in a mechanism called filter gain [14]. The effect of the noradrenergic system in increasing the signal-to-noise ratio facilitates the detection of signals embedded in a noisy environment [25], boosting the signal detection and promoting integration [13]. Therefore, a complex interaction between the cholinergic and noradrenergic system seems to manage the balance between integration and segregation.
Using a whole-brain model, Shine et al. [26] showed that neuromodulation and integration follow an inverted-U relationship. If one considers that neuromodulation also shows an inverted-U relationship with in-task performance, [15, 27], it is possible to hypothesize that neuromodulatory systems boost cognitive and attentional performance by increasing the functional integration in the brain, as proposed by Shine et al‥ [13, 26].
There are still unanswered questions about the specific effects of neuromodulation on integration and segregation. Experimental research points out that the cholinergic system, through both nicotinic and muscarinic receptors, boosts the signal-to-noise ratio in two principal ways [14, 28]: first, increasing the excitability of pyramidal neurons [29–31], and second, enhancing the firing rates of dendritic-targeting GABAergic interneurons –an effect that promotes a focused intra-columnar inhibition, reducing the local excitatory feedback to pyramidal neurons [31–33]. Consequently, pyramidal neurons become more responsive to stimulus from other distant regions respect to the stimulus of its own cortical column [28, 29, 34]. The particular effect of the cholinergic system on excitatory neurons was one of the focus of the whole-brain simulation work by Shine et al‥ [26]. However, the cholinergic modulation of inhibitory interneurons and its effect on the segregation/integration balance has not been analyzed at the whole-brain level. This is the main focus of the present work.
Here, we use an in silico approach to analyze the effect of neuromodulatory systems on functional integration in the brain, focusing on the cholinergic action in inhibitory interneurons. We combined a real human structural connectivity with the Jansen & Rit neural mass model of cortical columns [35, 36]. The mesoscopic properties of the model enable us to study more specifically the effects of neuromodulators in whole-brain dynamics. To make our simulations more comparable to experimental findings [3, 12, 37, 38], and also following Shine et al. [26], fMRI blood-oxygen-level-dependent (BOLD) signals were generated from the firing rates of pyramidal neurons. Integration and segregation were then assessed in the functional connectivity matrices derived from the BOLD-like signals, using a graph theoretical approach.
The neuromodulation was discerned in three components. First, we included an “excitatory gain”, which increases the inter-columnar coupling. In our model, this gain mechanism is mediated by the action of the cholinergic system on pyramidal neurons, having an indirect effect on their excitability [13, 14, 28]. Second, we added an “inhibitory gain”, also mediated by the cholinergic system, that controls the inputs from inhibitory to excitatory interneurons and reduces the local feedback excitation. This additional connection, well described in cortical columns [39, 40], represents a modification of the original neural mass model proposed by Jansen & Rit [35, 36]. Finally, we incorporated a “filter gain”, that increments the pyramidal neurons’ sigmoid function slope [14]. This gain mechanism is mediated by the noradrenergic system; it acts as a filter, decreasing (increasing) the responsivity to weak (strong) stimuli [23, 25], boosting signal-to-noise ratio. We show that the increase of the signal-to noise ratio, mediated both by the excitatory and filter gains, and the decrease in the feedback excitation, related to the inhibitory gain, promote functional integration.
We assessed the effect of the neuromodulatory systems using a whole-brain neural mass model of brain activity. In the model, each node corresponds to a brain area and is represented by a neural mass consisting of three populations [35, 36]: pyramidal neurons (that reside in cortical column layer V), excitatory interneurons (nearby pyramidal cells which reside in the same layer than the principal pyramidal population), and inhibitory interneurons (Fig 1A). Based on Silberberg & Markram [39] and Fino et al. [40], we have added a connection from inhibitory interneurons from excitatory interneurons (thick line in Fig 1A), allowing us to study the effect of its modulation by cholinergic influences (see below). The nodes are connected through a weighted and undirected structural connectivity matrix derived from human data [41], parcellated in 90 cortical and sub-cortical regions with the automated anatomical labeling (AAL) atlas [42] (Fig 1B). Connections between nodes are made by pyramidal neurons, considering that long-range projections are mainly excitatory [43, 44]. Using the firing rates of each node as inputs to a generalized hemodynamic model [45], we obtained fMRI BOLD signals from which we calculated integration and segregation of the resulting FC matrices.
Whole-brain neural mass model.
A) The Jansen & Rit model is constituted by a population of pyramidal neurons with excitatory and inhibitory feedback mediated by interneurons (INs). Each population is connected by a series of constants Ci. The outputs are transformed from average pulse density to average postsynaptic membrane potential by an excitatory (inhibitory) impulse response function hE(t) (hI(t)). Then, a sigmoid function S performs the inverse operation. Pyramidal neurons project to distant cortical columns, and receive both uncorrelated Gaussian-distributed inputs p(t) and inputs from other cortical columns z(t). Neuromodulation is constituted by three parameters, colored in red: excitatory gain α, which scales z(t), inhibitory gain β, which increases the inhibitory input to excitatory INs (thick line), and filter gain, r0, which modifies the slope of the sigmoid function in pyramidal neurons. B) Each node represents a cortical column, whose dynamics is ruled by the Jansen & Rit equations. Nodes are connected through a structural connectivity matrix, M C) Neuromodulation modifies the coupling between neurons and the properties of the input (average postsynaptic membrane potential) to output (average pulse density) sigmoid function. The cholinergic system modifies the global coupling and local inhibition. α amplifies the response of pyramidal neurons to other columns’ input; it also increases pyramidal neurons excitability. β amplifies the effect of inhibitory INs to excitatory INs, damping pyramidal cells excitability. The noradrenergic system increments the responsivity of pyramidal neurons to relevant stimuli respect to noise, as a filter, by increasing the slope r0 of their sigmoid function.
Following Shine et al. [26], we modeled the influence of the cholinergic and noradrenergic systems through the manipulation of the response and filter gains, respectively (Fig 1C). The principal difference in our approach is that we split the response gain in excitatory gain (long-range pyramidal to pyramidal coupling), α, and inhibitory gain (local inhibitory to excitatory interneurons coupling), β. While the excitatory gain boosts the pyramidal neurons responsivity to long-range inputs, and indirectly increases the pyramidal cells excitability, the inhibitory gain reduces the local excitatory feedback from interneurons. Finally, the filter gain r0 modifies the sigmoid function slope of pyramidal neurons, increasing its responsivity to relevant stimuli and boosting signal-to-noise ratio. Here, we studied the combined effect of the three gain mechanisms to understand how neuromodulatory systems shape the global neuronal dynamics in two different timescales: EEG-like and BOLD-like signals. Our hypothesis is that the inhibitory gain will play a significant role in increasing the likelihood of integration.
We first studied the combined influence of the excitatory and inhibitory response gains, by fixing r0 = 0.56 mV−1 (its default value) and then simulating neuronal activity at different combinations of α ∈ [0, 1] and β ∈ [0, 0.5]. Then, we analyzed the graph properties of the static (time-averaged) functional connectivity (sFC) matrices, obtained from the pairwise Pearson’s correlations of BOLD-like signals. Namely, we calculated the global efficiency Ew, a measure of integration defined as the inverse of the characteristic path length [46], and modularity Qw, a measure of segregation based on the detection of network communities or modules [46]. High values of Ew represent an efficient coordination between all pairs of nodes in the network, a signature of integration. In contrast, a high modularity Qw is associated to segregation [46].
Fig 2A shows that functional integration (Ew) is maximized in an intermediate region of the (α, β) parameter space; and that integration is accompanied by a decrease in the segregation (Qw). Also, the system undergoes a sharp transition crossing a critical boundary. The transitions between different regimes are better appreciated in Fig 2B, where we show a 1-D sweep of α at β = 0.25. Dashed lines at α = 0.3 and α = 0.8 correspond to points in the parameter space where drastic changes in dynamic properties of the network occur. Further, the global efficiency Ew follows an inverted-U relationship with the excitatory neuromodulation, suggesting (in our whole-brain model) that optimal levels of cholinergic neuromodulation maximize functional integration (see Discussion). Also, Qw peaks higher at the right critical boundary (dashed lines). A 1-D sweep of β at α = 0.5 (Fig 2C), shows an increase in integration crossing the critical transition at β = 0.1. These results are replicated using other measures of integration and segregation: the mean participation coefficient PCw, an integration metric that quantifies the between-modules connectivity, and the transitivity Tw, which accounts for segregation counting triangular motifs [46] (S1 Fig).
Network features in the (α, β) parameter space.
A) Global efficiency Ew (integration) and modularity Qw (segregation) of the graphs derived from the sFCs of the BOLD-like signals. B) Transitions through critical boundaries in the α axis, for a fixed β = 0.25. Transition points are represented by black dashed lines at α = 0.3 and α = 0.8. C) Transitions in the β axis, for a fixed α = 0.5, with a critical transition at β = 0.1.
The modulation of the inhibitory gain (β) shows an important effect on the integration and segregation properties of the whole network measured by the global efficiency and modularity, respectively. This could be due to the reduction of excitatory feedback only, or to a more specific effect of the newly introduced connection from inhibitory to excitatory interneurons. In the first case, we expect a similar effect by reducing the C1 parameter (see Fig 1A) because this also reduces the excitatory feedback loop of the cortical columns. As shown in S2 Fig, this is only partially the case. The reduction of the C1 connection weight –in the absence of the inhibitory to excitatory interneuron connection– enables the network to reach integration but in a smaller region of the parameter space and to a lower extent than the inhibitory modulation that we introduced in our model.
To show in more detail how each gain mechanism produces integrated or segregated functional network states, we present in Fig 3 some BOLD-like signals and their respective sFC matrices. We chose five tuples of (α, β) parameters, marked with the red circles in the Fig 3A. Functional integration measured by global efficiency is maximal in the middle (α = 0.5, β = 0.25), and segregation measured by modularity is maximized far away from this point (α = 0.25, β = 0.125, and α = 0.75, β = 0.375). In the extreme cases (α = 0, β = 0, and α = 1, β = 0.5) there is neither integration nor segregation; in the first case the network is disconnected, and in the second one the system crossed the second bifurcation point and pyramidal neurons are not oscillating (neurons are over-excited).
fMRI-like sFCs at different values of α and β.
A) The red circles represent pairs of (α, β) values in which different integration/segregation profiles can be observed. B-F) BOLD-like signals, and their respective sFC matrices, for the (α, β) values shown in A. The sFC networks evolve from neither integration nor segregation (B, the nodes are disconnected), to a more integrated sFC (C). In D the integration is maximal, and a further increase of both parameters produces a more segregated sFC matrix (E). Finally, in F there is neither integration nor segregation (the pyramidal neurons are over-excited). We shown only 120 s of BOLD-like signals, while sFC matrices were built with the full-length time series (600 s).
To further validate our model, we sought to reproduce the results of the neuromodulatory paradigm proposed by Shine et al‥ [13, 26]. We characterized the relationship between neuromodulation and integration in the (α, r0) parameter space, with α ∈ [0, 1] and r0 ∈ [0, 1] while leaving β fixed at 0 or 0.4 (without and with inhibitory gain, respectively). The results for β = 0 (Fig 4A) show no integration in the entire parameter space. On the other hand, the observations of Shine et al. [26] are fully reproduced with β = 0.4 (Fig 4B). Similar results hold for the mean PCw and Tw, as shown in S3 Fig.
Network features in the (α, r0) parameter space.
A-B) Global efficiency Ew (integration) and modularity Qw (segregation) of the graphs derived from the sFCs of the BOLD-like signals, for A) β = 0 (no action of the inhibitory gain) and B) β = 0.4. C) Transitions through the critical boundary in the α axis, with a fixed r0 = 1 mV−1 and β = 0.4. Critical transition points represented by black dashed lines at α = 0.3 and α = 0.8. D) Transitions in the r0 axis, for a fixed α = 0.6 and β = 0.4, with a critical transition at r0 = 0.33 mV−1.
As observed previously in Fig 2, critical boundaries delimit asynchronous and synchronous states in the (α, r0) parameter space. A 1-D sweep of α at r0 = 1 mV−1 shows a sharp transition (Fig 4C) (Fig 4B shows that this is also true for lower values of r0). Global efficiency Ew increases alongside the decrease of modularity Qw, and further increments of α produce network desynchronization. On the other hand, a 1-D sweep of r0 at α = 0.6 (Fig 4D) produces similar observations, but just one boundary is visible. As in the (α, β) parameter space, the excitatory gain α follows an inverted-U relationship with integration. This relationship was not observed between the filter gain r0 and integration.
In the whole-brain model, the cholinergic system exerts its effect by changing both α and β parameters. Under this assumption, a logical consequence of the cholinergic neuromodulation is the possibility of α and β increasing/decreasing simultaneously. For that reason, we repeated the analysis previously performed in the (α, r0) parameter space, but this time we changed β alongside α following the relationship β = 0.5α. The results are shown in S4 Fig. They are similar to those in Fig 4, but this time the relationship between r0 and Ew is no longer a sigmoid-like function, and instead it follows an inverted-U relationship.
Previous experimental and theoretical works [13, 14, 28] suggest that neuromodulatory systems increase the signal-to-noise ratio, allowing neuronal populations to be sensitive to local or distant populations to a greater extent than noise. To test that, we measured the signal-to noise-ratio (SNR) using the power spectral density (PSD) function of each EEG-like signal (see Methods) and report the average value over all nodes. Additionally, we computed the average Kuramoto order parameter [47], as a measure of global synchronization in the fast timescale of EEG. Values of closer to 1 indicate a perfect in-phase synchronization, and values closer to 0 indicate complete asynchrony.
Both SNR and match the region of integration, measured as the global efficiency Ew in the slowest BOLD timescale (Fig 5), supporting the idea that neuromodulatory systems promote integration by increasing SNR. In consequence, our results link in two different timescales the effect of neuromodulation in the coordination of brain activity. These results are not possible without the action of the inhibitory gain (β = 0, Fig 5A).
EEG features in the (α, r0) parameter space.
A-B) Average phase synchrony and signal-to-noise ratio (SNR) measured from the EEG-like signals. A) No action of inhibitory gain (β = 0). B) The increase of and the SNR matches with functional integration for β = 0.4. C) Transitions through the critical boundary in the α axis, with a fixed r0 = 1 mV−1. Critical transitions represented by black dashed lines at α = 0.3 and α = 0.8. D) Transitions in the r0 axis, for a fixed α = 0.6, with a critical transition at r0 = 0.33 mV−1.
As suggested by experimental [3] and computational studies [26], a shift to more segregated or integrated functional states decreases the topological variability of the network. Also, near the critical transitions for segregation to integration, network variability and communicability peak [26]. We tested this hypothesis performing a Functional Connectivity Dynamics (FCD) analysis [9, 10] on the BOLD-like signals, using the sliding windows approach depicted in Fig 6A–6C [48]. The resulting time vs time FCD matrix captures the concurrence of FC patterns, visualized as square blocks. We computed the variance of the FCD, var(FCD), as a multistability index [48], where values greater than 0 indicate the switching between different FC patterns. Additionally, we calculated the FCD speed dtyp as described by Battaglia et al. [49], which captures how fast the FC patterns fluctuate over time. Values closer to 1 indicate a continuous change of diverse FC patterns, and closer to 0 the concurrence of stable and similar states over time.
Analysis of functional connectivity dynamics.
A) Sample fMRI BOLD time series showing the fixed length and overlapping time windows at the begginging. In color, the time windows corresponding to the FCs shown in B. B) FC matrices obtained in the colored time windows. C) Functional Connectivity Dynamics (FCD) matrix, where all the FCs obtained were vectorized and then compared against each other using a vector-based distance (Clarkson distance). D) FCD matrices through the critical boundary, in both α and r0 direction. Below each FCD, a histogram of its upper triangular values is shown. The variance of these values constitutes a measure of multistability.
In Fig 6D we show a set of FCD matrices obtained at different values of α and r0, together with histograms of their off-diagonal values. Red FCD matrices (with high values) correspond to incoherent states, as the FC continuously evolves in time. On the other hand, a blue FCD matrix (with low values) indicates a fixed FC throughout the simulation. Multistability is higher for green/yellow patchy matrices, because this indicates FC patterns that change and also repeat over time. As can be inferred observing the FCD distributions, the variance of the values in the histograms –var(FCD)– can be used as a measure of multistability [48].
Fig 7 shows how multistability (var(FCD)) and FCD speed change in the whole (α, r0) space, for β = 0.4. At low levels of both α and r0, the neuronal activity is constituted mainly by noisy asynchronous signals, conditions associated to low (near 0) values of var(FCD), and with a high dtyp (all FC patterns differ from each other, as expected for noise-driven signals) (Fig 7A). In the other extreme, for r0 > 0.5 mV−1 and α ∈ [0.5, 0.6], values that correspond to the integrated states, var(FCD) is also small and dtyp falls close to 0. In consequence, integrated states are more stable and less susceptible to network reconfiguration over time. In contrast, var(FCD) peaks near the critical boundary, through the α and r0 axes (Fig 7B and 7C). Moreover, crossing the boundaries is associated with a continuous decrease of dtyp: the emerging integration mediated by gain mechanisms is associated with more stable FC patterns over time.
Dynamical features of the system in the (α, r0) parameter space for β = 0.4.
A) Multistability var(FCD) and typical FCD speed dtyp measured from the Functional Connectivity Dynamics analysis (BOLD-like signals). B) Transitions through the critical boundary in the α axis, with a fixed r0 = 1 mV−1. Critical transitions represented by black dashed lines at α = 0.3 and α = 0.8. C) Transitions in the r0 axis, for a fixed α = 0.6, with a critical transition at r0 = 0.33 mV−1.
In this work, we used a whole-brain neural mass model to investigate how local (meso-scale) neuromodulatory effects can impact the global functional network properties. Importantly, we studied the effect that the cholinergic system has in both, excitatory and inhibitory neurons, along with the noradrenergic modulatory influence. Our model shows an increase in functional integration at intermediate values of the parameters that resemble the cholinergic and noradrenergic systems, following an inverted-U relation with the neuromodulation. In addition, the modulation of an intra-columnar inhibitory gain can promote functional integration and facilitates the effect of the other neuromodulatory systems. Finally, we show that our results hold for both the EEG and fMRI timescales, and that integration is accompanied by an increment in the signal-to-noise ratio, as well as a reduction of dynamical variability captured by the FCD analysis.
Our main motivation was to study the large-scale effects of the cholinergic neuromodulation of local inhibitory circuits. Although the cholinergic system increases the response of pyramidal neurons to external afferences through nicotinic and muscarinic receptors [34], the same system can promote intra-columnar inhibition, an effect mediated by nicotinic receptors expressed by inhibitory interneurons [28, 31, 50]. A possible consequence is the increase of the influence of external inputs, in comparison with the local intra-columnar inputs, shifting the flow of information from local to global processing. Based on these experimental findings [28, 31, 34, 50], we hypothesized that the cholinergic neuromodulation of the inhibitory interneurons facilitates functional integration. Our model shows that the action of the cholinergic system, on both the excitation of pyramidal neurons and the intra-columnar inhibitory feedback, can shift the system towards a functionally integrated regime. In this way, we propose a plausible biophysical mechanism of inhibitory-to-excitatory interneuron connection that facilitates functional integration of brain activity (see Fig 4).
The new intra-columnar connection that we introduced in the Jansen & Rit model –represented by the inhibitory gain β–, produces a higher dampening of the excitatory input to pyramidal neurons when their excitability is high. Conversely, when the pyramidal excitability decreases, the effect of the inhibitory loop between interneurons is low, and the excitatory loop can rise the excitability of pyramidal cells. In this way, the inhibitory gain provides a simple dynamical mechanism to homeostatically preserve the excitation/inhibition (E/I) balance at the node level. In contrast, a simple reduction of the excitatory feedback (e.g., decreasing C1, S2 Fig), fails to compensate the E/I balance when pyramidal excitability is low, and thus has a limited ability to promote functional integration. This highlights the role of specific intra-columnar inhibitory feedback connections in shaping the network behavior, and justifies our modification of the model with a homeostatic mechanism. Similar types of inhibition-mediated control of the E/I balance have been implemented in a dynamic mean-field model [51] as well as in the Wilson-Cowan model [52]. Remarkably, the E/I balance has been considered a determinant element in the interplay between integration and segregation [53]
Our results also have an interpretation from the dynamical systems perspective. It has been proposed that, at rest, brain activity operates near a bifurcation point, where segregated (uncoordinated) and integrated (coordinated) regimes alternate in time. Then, a shift to more segregated or integrated states takes places with a change in behavioral context [5, 9, 10]. At the node level, the Jansen & Rit model has two Hopf supercritical bifurcations [54]. When α and r0 are low, the node dynamics is defined by a stable focus (a fixed point with non-monotonic convergence), and thus pyramidal outputs consist of low amplitude noisy signals. Increasing both parameters causes the bifurcation into an unstable focus within a limit cycle, with high amplitude oscillations. Increasing α further produces a new bifurcation (into a stable focus) and the limit cycle disappears. The inhibitory gain β constitutes a mechanism to keep the model working between the two bifurcations points, allowing the transitions between different functional states (more segregated or integrated). This again highlights the role of β as an inhibitory control loop which preserves the E/I balance and sustains a richer brain dynamics.
Several clues suggest that the model we propose is in the right track. First, our model reproduces previous results of Shine et al., also using a whole-brain model [26]: they reported and inverted-U relationship between cholinergic neuromodulation and integration, an increase in phase synchronization, and that integration is accompanied by a reduction in the time-resolved topological variability (captured, in our analysis, by the variance of the FCD). Second, functional integration matches with an increase of SNR, a common effect attributed to neuromodulatory systems [13, 14, 28, 55]. Third, integration is accompanied in our model by a reduction in oscillatory frequency –which falls within the Theta range of EEG spectrum– (S8 Fig); an effect that is also perceived in several cognitive tasks [56, 57], and reproduced by other computational studies [58, 59].
From an experimental point of view, an increase in the local and global efficiency (integration) in fMRI has been reported after the administration of nicotine both in resting state conditions [37] and during an attentional task [38]. Interestingly, the performance was positively correlated with the global efficiency, and negatively correlated with the average clustering coefficient. Some nicotinic agonists have pro-cognitive effects as well, in health and disease [60]. Considering the relationship between functional integration and cognition [2, 3, 13], our model suggests that the possible pro-cognitive effects associated with the cholinergic system are due to a selective increase in the excitability of excitatory and inhibitory neural populations within brain areas. Thus, our computational approach –in the same spirit as Wylie et al. [37] and Gießing et al. [38]– links the meso-scale consequences of inhibitory interneurons neuromodulation with the functional network topology features at the whole-brain level.
On the other hand, the inverted-U relationship between neuromodulation and integration that we are reporting in our whole-brain model, has not been experimentally observed. However, an inverted-U function between cholinergic and noradrenergic neuromodulation and in-task performance [15, 27]; as well as between in-task performance and functional integration [12] has been reported. Taking these together we propose, in accordance to Shine et al‥ [13, 26], that neuromodulation improves cognitive performance by boosting integration. The results from computational modeling should nevertheless be verified by measuring both integration and cognitive performance as functions of neuromodulatory activity, e.g., as functions of the dose of a cholinergic/noradrenergic drug.
There is a lot of room for further progress starting from this work. Future research may consider the addition of neuromodulatory maps [41, 61] in order to take into account the heterogeneous expression of the receptors, or explore models that can reproduce the effect of other neuromodulatory systems [17]. Furthermore, it is known that cholinergic and noradrenergic projections have a specific spatial organization [62, 63]. Our model considers neuromodulation to be static, that is, the parameters α, β and r0 do not change over time, as in tonic neuromodulation. An improvement to our model may be the addition of the release and reuptake dynamics of neuromodulators, as in Kringelbach et al. [64] or the characterization of the dynamics under acute neuromodulatory ‘pulses’. Other interneurons subtypes and their modulation could be included –such as fast-spiking inhibitory interneurons– to account for the faster EEG features of brain activity [65]. Finally, the graph theoretical analysis used here only considers pairwise interactions, neglecting high-order effects that may contain important information about high dimensional functional brain interactions. Information-theoretical [66, 67] and algebraic topological approaches [68–70] may provide complementary insights of high-order interdependencies in the brain.
Our findings shed light on a better understanding of neurophysiological mechanisms involved in the functional integration and segregation of the human brain activity and constitutes a step forward from the neuromodulatory framework proposed by Shine [13], including the role of a second cholinergic target and also highlighting the role of a homeostatic inhibitory feedback. This line of research may have plentiful of scientific and clinical implications, as a vast body of evidence suggest that functional integration and segregation may be altered in neuropsychiatric disorders [53, 71, 72], e.g, in Alzheimer disease and Attention-Deficit/Hyperactivity Disorder (ADHD) [73, 74]. These results point out the usefulness of graph theoretical analysis to exhibit functional markers for characterizing and understanding neuropsychiatric disorders. Understanding the neuromodulatory mechanisms that underlie the imbalances of integration and segregation will lead to a more profound understanding of how the brain works in health and disease and to future progress in pharmacological treatments.
To simulate neuronal activity we used a modified version of the Jansen & Rit neural mass model [35, 36]. In this model, a cortical column consists of a population of pyramidal neurons (that reside in cortical column layer V) with projections to other two populations: excitatory interneurons (nearby pyramidal cells which reside in the same layer than the principal pyramidal population) and inhibitory interneurons; both project back to the pyramidal population. The dynamical evolution of the three populations within the cortical column is modeled by two blocks each. The first transforms the average pulse density in average postsynaptic membrane potential (which can be either excitatory or inhibitory) (Fig 1A). This block, denominated post synaptic potential (PSP) block, is represented by an impulse response function
Additionally, the pyramidal neurons receive an external stimulus p(t), whose values were taken from a Gaussian distribution with mean μ = 2 impulses/s and standard deviation σ = 2. Different values of σ were explored; qualitatively the results are similar for different σ values, but the magnitude of integration decreases with σ. This exploration is shown in the S5 Fig. In the same manner, the mean of the Gaussian distribution μ has an effect in decreasing the synchronization and integration, as shown in the S5 Fig.
To study the effect of the neuromodulatory systems at the macro-scale level, we included long-range pyramidal-to-pyramidal neurons and short-range inhibitory-to-excitatory interneurons couplings, to mimic the effects of neuromodulation through the excitatory and inhibitory gain parameters, respectively. This short-range coupling between interneurons, well described at the meso-scale level [39, 40], constitutes a modification of the original equations. A bifurcation analysis of the model at different values of β reveals that this parameter shifts the oscillatory regime of the system towards larger values of the external input p(t) (S7 Fig). Nevertheless, the main oscillatory rhythm of the model is well maintained within the α band frequency (around 10 Hz).
In the model presented in the Fig 1A, each node i ∈ [1…N] represents a single brain area. The nodes are connected by a normalized structural connectivity matrix (Fig 1B). This matrix is derived from a human connectome [41] parcellated in n = 90 cortical and subcortical regions with the automated anatomical labelling (AAL) atlas [42]; the matrix is undirected and takes values between 0 and 1. Because long-range connections are mainly excitatory [43, 44], only links between the pyramidal neurons of a node i with pyramidal neurons of a node j are considered. We applied a local normalization procedure to the structural connectivity matrix M. The normalization consisted in dividing all the columns belonging to a node i by the in-strength of the node. The entries of the resulting normalized matrix are defined as
The local normalization procedure constitutes a form of homeostatic plasticity, which equalizes the excitatory inputs that the nodes receive, while preserving the structural topology. It has been reported that this mechanism improves the fit of a whole-brain mesoscopic model to empirical fMRI data, and leads to a better estimation of the functional connectivity [75].
The overall set of equations, for a node i, includes the within and between nodes activity
The overall input from other cortical columns j ≠ i to the column i is given by
The average PSP of pyramidal neurons in column i characterizes the EEG-like signal in the source space; it is computed as [35, 36]
The effects of the cholinergic system were modeled by the parameters α and β. The parameter α increases the long-rage pyramidal to pyramidal neuron coupling through the matrix. Although α does not control directly the excitability, increasing α amplifies the input to pyramidal neurons [13, 14]. The parameter β scales the short-range inhibitory-to-excitatory interneurons coupling, decreasing the recurrent excitation to pyramidal neurons [28]. We refer to α as the excitatory gain, and β as the inhibitory gain. In comparison with the current framework proposed by Shine [13], the novelty of our neuromodulatory approach is the inclusion of the inhibitory gain to the model. The effect of the noradrenergic system, designated as filter gain, was simulated controlling the parameter r0, which represents the sigmoid function slope of the pyramidal population, and increases the signal-to-noise ratio of pyramidal cells [14, 25].
Following Birn et al. [77], we ran simulations to generate the equivalent of 11 min real-time recordings, discarding the first 60 s. The system of differential equations (Eq (5)) was solved with the Euler–Maruyama method, using an integration step of 1 ms. We used six random seeds which controlled the initial conditions and the stochasticity of the simulations. We simulated neuronal activity sweeping the parameters α ∈ [0, 1], β ∈ [0, 0.5] and r0 ∈ [0, 1]. All the simulations were implemented in Python and the codes are freely available at https://github.com/vandal-uv/Neuromod2020.
We used the firing rates ζi(t) to simulate BOLD-like signals from a generalized hemodynamic model [45]. An increment in the firing rate ζi(t) triggers a vasodilatory response si, producing blood inflow fi, changes in the blood volume vi and deoxyhemoglobin content qi. The corresponding system of differential equations is
The system of differential equations (Eq (8)) was solved with the Euler method, using an integration step of 1 ms. The signals were band-pass filtered between 0.01 and 0.1 Hz with a 3rd order Bessel filter. These BOLD-like signals were used to build the functional connectivity (FC) matrices from which the subsequent analysis of functional network properties is performed using tools from graph theory.
Although the nodes consist of three neural masses, there is some evidence that the hemodynamic response is related mainly to excitatory activity [78]. In fact, some reports suggest that inhibitory activity does not trigger a measurable BOLD response, because the inhibitory connections are relatively few and their energy expenditure is lower [79]. Nevertheless, we reproduced the Fig 2 using a combined BOLD response, and we found no noticeable differences (see S6 Fig).
As a measure of global synchronization in the EEG timescale, we calculated the Kuramoto order parameter R(t) [47] of the EEG-like signals ν(t) derived from the Jansen & Rit model. First, the raw signals were filtered with a 3rd order Bessel band-pass filter using their frequency of maximum power (usually between 4 and 10 Hz) ±3 Hz. Then, the instantaneous phase ϕ(t) was obtained with the Hilbert transform.
The global phase synchrony is computed as:
We measured the average signal-to-noise ratio (SNR) over all raw signals and nodes, using the power spectral density function denoted PSD(ω). This function was calculated using the Welch’s method [80], with 20 s time windows overlapped by 50%. We excluded the 2nd to 5th harmonics [81]. For a node i, the signal power, Psignal, was measured as the area under the curve of PSD(ω) within ωi ± 1Hz. Noise power, Pnoise, corresponds to the area under the curve of PSD(ω) outside the ±1Hz window. Then, the SNR was calculated as
The static Functional Connectivity (sFC) matrices were built from pairwise Pearson’s correlations of the entire BOLD-like time series. Instead of employing an absolute or proportional thresholding, we thresholded the sFC matrices using Fourier transform (FT) surrogate data [82] to avoid the problem of introducing spurious correlations [83]. The FT algorithm uses a phase randomization process to destroy pairwise correlations, preserving the spectral properties of the signals (the surrogates have the same power spectrum as the original data). We generated 500 surrogates time series of the original set of BOLD-like signals, and then built the surrogates sFC matrices. For each one of the (n2 − n)/2 possible connectivity pairs (with n = 90) we fitted a normal distribution of the surrogate values. Using these distributions we tested the hypothesis that a pairwise correlation is higher than chance (that is, the value is at the right of the surrogate distribution). To reject the null hypothesis, we selected a p-value equal to 0.05, and corrected for multiple comparisons with the FDR Benjamini-Hochberg procedure [84] to decrease the probability of make type I errors (false positives). The entries of the sFC matrix associated with a p-value greater than 0.05 were set to 0. The result is a thresholded, undirected, and weighted (with only positive values) sFC matrix.
Integration and segregation were evaluated over the thresholded sFC matrices. We employed the weighted versions of transitivity [85] and global efficiency [86] to measure segregation and integration, respectively. A detailed description of the metrics used can be found in Rubinov & Sporns [46]. The transitivity (similar to the average clustering coefficient) counts the fraction of triangular motifs surrounding the nodes (the equivalent of counting how many neighbors are also neighbors of each other), with the difference that it is normalized collectively. It is defined as
We also calculated other two measures of integration and segregation: the participation coefficient PCw and modularity Qw, respectively, both based on the detection of the network’s communities [46]. The detection of so-called communities or network modules in the thresholded sFC matrix, was based on the Louvain’s algorithm [87, 88]. The algorithm assigns a module to each node in a way that maximizes the modularity (14). We used the weighted version of the modularity [89] defined as
Finally, we computed the weighted version of the participation coefficient [91]. This metric quantifies, for each individual node, the strength of between-module connections respect to the within-module connections, and is defined as
The FCD matrix captures the evolution of FC patterns and, consequently, the dynamical richness of the network [9, 10]. We used the sliding window approach [9, 48] depicted in Fig 6. Window length was set to 100 s with a displacement of 2 s between consecutive windows (Fig 6A). The length was chosen on the basis of the lower limit of the band-pass filter (0.01 Hz), in order to minimize spurious correlations [92]. For each window, an FC matrix was calculated from the pairwise Pearson’s correlations of BOLD-like signals (neglecting negative values), thus we obtained 251 weighted and undirected FC matrices from the 600 s simulated BOLD-like signals (Fig 6B).
The upper triangle of each FC matrix is unfolded to make a vector, and the FCD is built by calculating the Clarkson angular distance [93] between each pair of FCs (Fig 6C)
The speed of the FCD was measured as described by Battaglia et al. [49]. We computed the histogram of FCD values through a straight line from FCD(τ, 0) to FCD(tmax, tmax − τ), with tmax = 600 s as the total time-length of the signals and τ = 100 s. The median of the histogram distribution corresponded to the typical FCD speed dtyp. Values closer to 1 indicate a constant switching of states, and values closer to 0 correspond to stable FC patterns.
We thank to Gustavo Deco who kindly provided the anatomical connectivity matrix used in the model. We also want to thank to Chiayu Chiu and Andrés Chávez for their feedback and suggestions about the manuscript.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
Cholinergic neuromodulation of inhibitory interneurons facilitates functional integration in whole-brain models
Facebook
Twitter
Linkedin
Whatsapp
