## Abstract

The mammalian brain operates in multiple spatial scales simultaneously, ranging from the microscopic scale of single neurons through the mesoscopic scale of cortical columns, to the macroscopic scale of brain areas. These levels of description are associated with distinct temporal scales, ranging from milliseconds in the case of neurons to tens of seconds in the case of brain areas. Here, we examine theoretically how these spatial and temporal scales interact in the functioning brain, by considering the coupled behaviour of two mesoscopic neural masses (NMs) that communicate with each other through a microscopic neuronal network (NN). We use the synchronization between the two NM models as a tool to probe the interaction between the mesoscopic scales of those neural populations and the microscopic scale of the mediating NN. The two NM oscillators are taken to operate in a low-frequency regime with different peak frequencies (and distinct dynamical behaviour). The microscopic neuronal population, in turn, is described by a network of several thousand excitatory and inhibitory spiking neurons operating in a synchronous irregular regime, in which the individual neurons fire very sparsely but collectively give rise to a well-defined rhythm in the gamma range. Our results show that this NN, which operates at a fast temporal scale, is indeed sufficient to mediate coupling between the two mesoscopic oscillators, which evolve dynamically at a slower scale. We also establish how this synchronization depends on the topological properties of the microscopic NN, on its size and on its oscillation frequency.

## 1. Introduction

The mammalian brain is composed of a myriad of coupled neurons that interact dynamically. It possesses a rich topological structure and exhibits complex dynamics, operating as a noisy, nonlinear and highly dimensional system. Neuronal activity evolves at temporal scales ranging from a few milliseconds to tens of seconds, and emerges from neuronal assemblies that extend from micrometres to several centimetres. Owing to a complex functional hierarchy between cell groups, the brain is able to store information for long periods of time, process multiple sensory inputs efficiently and produce coherent output in the form of actions and thoughts.

Even though the brain has been studied for centuries, a full theoretical description of its normal and pathological functioning is still missing. Owing partly to the lack of a full description of the anatomical connectivity, and partly to our incomplete knowledge of the interplay between different neural processes, the brain is still the great unknown organ. Its study is usually partitioned into different research fields devoted to distinct brain structures (such as the thalamus, amygdala, hippocampus, etc.), cortical functional areas (motor, visual, auditory cortex, etc.) or particular microscopic circuits, from the level of cortical columns down to single-neuron responses. Moreover, studies of the global activity of the brain usually focus for convenience on specific cognitive or motor tasks, in order to compare them with a control state such as spontaneous activity at rest.

The various aforementioned approaches deal with different scales of description, from the macroscopic to the microscopic level. Accordingly, different computational models have been developed to account for the activity at each scale. Single neurons, for instance, can be characterized by detailed biophysical models that consider ion-channel dynamics, as initially proposed by Hodgkin & Huxley [1,2], or by more abstract models of neural excitation such as the integrate-and-fire model [3,4] or the FitzHugh–Nagumo model [5,6]. The set of equations representing each neuron's membrane potential can be coupled in a way that mimics the synaptic junction. Thus, given a connectivity matrix, one can ideally build any neuronal network (NN) *in silico* from its individual constituents, and thereby move towards the mesoscopic level of neuronal assemblies. This allows the brain to be traditionally investigated in a reductionist way, using different simplified levels of description. This approach has been very fruitful in unveiling several mechanisms that lay at the basis of the observed neural tissue behaviour [7–10].

Another set of models, named neural mass (NM) models [11–14], avoid the single-neuron perspective and consider instead the averaged behaviour of the neuronal population. This mesoscopic description is more phenomenological than the single-neuron models, in the sense that it represents directly the collective behaviour of the network, without singling out individual cells. Moreover, single neurons operate at time scales faster than NM models. The former exhibit action potentials that last about 1 ms, while the coordinated activity of neuronal tissue, which emerges from the synchronization of multiple spikes, operates on time scales up to tens of seconds. Within a neuronal population all temporal scales work simultaneously, and the relative relevance of the different scales might change depending on the biological process. For instance, spike-timing precision is key to synaptic plasticity, and therefore to the formation of functional cell assemblies [15,16]. On the other hand, the frequency of collective oscillations is relevant for the synchronization of distant areas, and thus for their effective interaction within specific information-processing tasks [17,18].

Recently, large-scale models of the brain have received special attention. So far, global brain activity has been modelled by dividing the brain into discrete volume elements, or voxels, and coupling them according to statistical correlations and structural information [19–21]. Both the Human Brain Project and the Brain Activity Map project propose integrated views to bridge the gap between the behaviour of single neurons and the functions of the full brain [22], but this quest is still in its infancy.

While new theoretical studies have attempted to connect the microscopic (NN) and mesoscopic (NM) descriptions of brain tissue, by directly applying mean field approaches to derive the latter from the former [23,24], these strategies are fraught with limitations and hard-to-justify assumptions. Here, we propose an alternative approach to explore scale interaction, by considering a system formed by two NMs that are coupled exclusively via an intermediate population described by a spiking NN model. Our results show that the two mesoscopic oscillators are able to lock their dynamics through the mediation of the microscopic population. In that way, we use synchronization as a tool to probe the interaction between the two scales of description. We also examine which characteristics of the NN connectivity allow the efficient cross-talk between dynamical scales, i.e. to determine whether which are the microscopic features that modulate the mesocopic activity. The paper is organized as follows. In §2, we describe the model used in this study. The results are presented in §3 and discussed in §4.

## 2. Dynamical model and methods

As mentioned above, our model combines two different levels of description (figure 1). The NM description evolves at a slow scale and represents the average dynamical evolution of a set of three different neural populations (pyramidal, excitatory interneurons and inhibitory interneurons) [11]. The fast scale, on the other hand, is represented by a conductance-based neural network formed by excitatory and inhibitory neurons. In this case, the time course of every neuron's transmembrane potential is given by the dynamics of voltage-dependent ion channels. We have merged these two levels of description in a simple dynamical structure, shown in figure 1, in which two NM models are coupled with a subpopulation of neurons belonging to the NN.

The two NMs are set to oscillate in two different well-defined frequencies, corresponding to two slightly different brain rhythms. The NN also displays a collective oscillatory dynamics with a different frequency. Here, we investigate how both the inter-scale coupling strength and the features of the NN contribute to the cross-talk between the three systems.

### (a) Dynamics of the neuronal network

The neural network is composed of 4000 neurons (80% excitatory and 20% inhibitory). Each neuron forms 400 chemical synaptic connections on average with other neurons of the network. The dynamics of the transmembrane potential of the soma of each neuron is described by the following set of differential equations:
2.1
where *g*_{K}, *g*_{Na} and *g*_{L} are the maximum conductances for the potassium, sodium and the leak currents, respectively, and *I*_{syn} is the synaptic current coming from the neighbouring neurons. The dynamics of the sodium and potassium channels is represented by the time-varying probabilities of a channel being in the open state:
2.2
where *x* stands for the activation (*m*) and inactivation (*h*) of the sodium channels and the activation of the potassium channels (*n*). The rate functions *α _{x}* and

*β*for each gating variable, together with all the NN parameters used throughout this paper are given in [10].

_{x}The synaptic current *I*_{syn} is described using a conductance-based formalism
2.3
where *g*_{syn} is the synaptic conductance, and *E*_{syn} is the reversal potential of the synapse. For *E*_{syn} greater than the resting potential *V*_{rest} the synapse is excitatory (mediated by AMPA receptors); otherwise it is inhibitory (mediated by GABA receptors). We consider two temporal time constants, *τ*_{d} and *τ*_{r} (decay and rise synaptic time) for the dynamics of the synaptic conductance
2.4
We have chosen the maximal conductances such that the postsynaptic potential (PSP) amplitudes are within physiological ranges: the excitatory postsynaptic potential (EPSP) in the range from 0.42 to 0.83 mV and the inhibitory postsynaptic potential (IPSP) in the range from 1.54 to 2.20 mV. In order to modify the activity time scale of the NN, we have changed *τ*_{d} for the GABAergic synapses, varying accordingly the inhibitory conductances in such a way that the maximum amplitude for *g*_{syn}(*t*) is maintained.

All neurons receive an additional train of excitatory presynaptic potentials, coming from brain areas other than those explicitly modelled by the NMs, which contributes to the external current term *I*_{ext} in equation (2.1). Those spikes follow a heterogeneous Poisson process with a mean event rate, which varies following an Ornstein–Uhlenbeck (OU) process. The instantaneous rate *λ*(*t*) of this external excitatory train of spikes is generated according to
2.5
where *σ*(*t*) is the standard deviation of the noisy process and is set to 0.6 spikes s^{–1}. The correlation time *τ* is set to 16 ms, leading to a 1/*f* power spectrum for the *λ* time series that is flat up to a cut-off frequency *f* = 1/(2*π**τ*) = 9.9 Hz. The term *η*(*t*) is a Gaussian white noise.

This NN model is able to reproduce the well-known synchronous irregular regime [25], in which recurrent activity leads to collective oscillations at the population level while single neurons fire irregularly. The emergent rhythmicity is achieved by a balance between the excitatory and inhibitory synaptic currents and can be explained by periodic changes of the excitability in the network, i.e. periodic modulation of the distance to threshold. Despite the fact that excitatory neurons are dominant in the network, the stronger synaptic inhibitory conductances and the higher firing rate of the inhibitory neurons allow the system to reach a balance between excitation and inhibition. In order to obtain collective oscillations in the alpha (gamma) band, we set the decay synaptic time to be *τ*_{d} = 15 ms (5 ms).

### (b) Dynamics of the neural mass model

The description of the mesoscopic neuronal ensemble is based on a model proposed by Jansen & Rit [11]. This model characterizes the dynamics of a cortical column by using a mean field approximation. In this sense, Jansen's model describes the average activity of three cortical populations; excitatory and inhibitory interneurons and pyramidal cells. All three populations form a feedback circuit. The main pyramidal population excites both interneuronal populations in a feedforward manner and the excitatory (inhibitory) interneurons feed back in an excitatory (inhibitory) manner into the pyramidal population. The dynamical evolution of these three populations is introduced considering two different transformations. Each population transforms the total average density of action potentials reaching their afferent synapses from different origins, ∑_{m}*p _{m}*(

*t*), into an average postsynaptic excitatory or inhibitory membrane potential

*y*(

_{i}*t*). This transformation can be introduced in the model using the differential operator 2.6 and correspondingly

*L*(

*y*)(

_{i}*t*);

*b*) for the inhibitory integration of the average density of action potentials, with

*b*and

*B*substituting

*a*and

*A*above.

*A*and

*B*are related with the maximum heights of the EPSP and IPSP, respectively, whereas

*a*and

*b*represent the inverse of the membrane time constants and the dendritic delays. Here,

*A*= 3.25 mV,

*a*= 0.1 kHz,

*B*= 22 mV and

*b*= 0.050 kHz for NM1 (figure 1), leading to oscillations at approximately 11 Hz; and

*A*= 3.25 mV,

*a*= 0.1 kHz,

*B*= 21.5 mV and

*b*= 0.048 kHz for NM2, which consequently exhibits oscillations at approximately 4.5 Hz.

The second dynamical transformation in the model is the conversion of the net average membrane potential into an average density of spikes. This conversion is done at the somas of the neurons that form the population and is described mathematically by a sigmoidal function
2.7
Here *e*_{0} = 0.0025 kHz determines the maximum firing rate of the neural population, *ν*_{0} = 0.6 mV sets the net PSP for which a 50% firing rate is achieved, *r* = 0.56 mV^{−1} is the steepness of the sigmoidal transformation, and *m*(*t*) corresponds to the net PSP input into the population being considered. The average density of action potentials produced by the presynaptic population acting upon the postsynaptic population, *p _{i}*(

*t*), turns out to be proportional to

*S*(

*m*(

*t*)), where the proportionality constant weights the contact between the populations, and gives the range of efficiency of the synaptic interaction. Combining equations (2.6) and (2.7), we obtain the complete model for the NMs 2.8 where

*y*

_{0}(

*t*) is the EPSP produced by the pyramidal population on the interneuron populations, and

*y*

_{1}(

*t*) is the EPSP acting upon the pyramidal population and arriving from (i) the excitatory interneurons, (ii) other areas of the brain and, in our model, (iii) the neural network (see equation (2.11)). Finally,

*y*

_{2}(

*t*) is the IPSP acting upon the pyramidal population and arriving from the inhibitory interneurons and, in our model, the neural network (see equation (2.12)). The intra-columnar connectivity constants values are defined in terms of

*C*, with

_{i}*i*= 1, … ,4. We use the values given in Jansen

*et al*. [11].

### (c) Inter-scale coupling terms

The effect of the mass models upon the neural network also contributes to the *I*_{ext} term of the NN (see equation (2.1)), together with the external excitatory Poissonian train of spikes. Hence, each neuron of the NN receives a train of excitatory spikes whose mean firing rate, FR, is given by
2.9
where *S*(*y*_{1}(*t*) − *y*_{2}(*t*)) translates the PSP of the pyramidal population of the NM which affects that particular neuron (or both NMs if that is the case) into a spiking rate. *γ*_{1} and *k* control the strength of this coupling. Here *γ*_{1} = 200, while *k* will be a varying parameter. EFR(*t*) corresponds to aforementioned Poissonian train of spikes:
2.10
with being the mean external firing rate and *λ*_{OU}(*t*) an OU process (see equation (2.5)) representing the fluctuations around the mean. We have set . The NN acts upon the NM models through *p*_{e}(*t*) and *p*_{i}(*t*) (see equation (2.8))
2.11
and
2.12
where is the multiunit activity coming from the neural network, i.e. the sum of spikes over the subset of neurons coupled to the NMs, calculated within a sliding window of length 1 ms. is a constant input coming from other areas of the brain distinct from those considered explicitly in our model ( for both NMs). *γ*_{2}, *γ*_{3} and *k* are scaling factors that take into account the synaptic efficiency, where, *γ*_{2} = 25 and *γ*_{3} = 3. Note that we assume that NN neurons affect only the pyramidal population in the NM. This is in accordance with the Jansen–Rit model of two coupled NMs [11], which considers that only pyramidal cells receive excitatory input from the other column.

### (d) Local field potential

In order to quantify the activity of the NN, we have defined a collective measure, the local field potential (LFP), as the average of the absolute values of AMPA and GABA synaptic currents acting upon a typical excitatory neuron [26]
2.13
where *I*_{AMPA} represents the external excitatory heterogeneous Poisson spike train and the recurrent excitatory synaptic current due to the network, *I*_{GABA} accounts for the recurrent inhibitory synaptic current and *R*_{e} represents the resistance of a typical electrode used for extracellular measurements. The symbol represents an average over all excitatory neurons.

### (e) Spectral analysis

We have computed the power spectral density (PSD) of the LFPs and of the PSPs of the pyramidal population of the NMs using the Welch method: the signal is split up into 500-point segments with 50% overlap and filtered with a Hamming window. Given that the sampling time is 1 ms, the segment length chosen leads to a frequency resolution of 2 Hz. The periodogram is calculated by computing the square magnitude of the discrete Fourier transform. These periodograms are then averaged to obtain the PSD estimate, which reduces the variance of the individual power measurements. LFP and PSP spectral quantities are averaged over 20 trials of length 3500 and 1500 ms, respectively. The code has been implemented in MATLAB. The frequency mismatch between the PSP of the two NM models is calculated as the inverse of the difference between the periods of the mass models, averaging over trials. The periods correspond to the temporal distance between two maxima of the autocorrelation function. Each trial corresponds to a different set of random initial conditions, NN architecture and realization of the OU process.

## 3. Results

The effective interaction between neuronal ensembles described at different scales can be studied by coupling mesoscopic and microscopic models. As mentioned in §1, mesoscopic models are best exemplified by NM descriptions, which are derived phenomenologically from experimental studies, and characterize the average population activity by means of a mean field approximation. In particular, NMs describe the neuronal activity happening at slow time scales, such as synaptic potentials arising from the synchronized firing of thousands of neurons. On the other hand, models of single neurons reproduce the time course of the electric currents crossing the neuronal membrane, and thus account for the individual action potentials and the postsynaptic response of each cell composing the network. In order to preserve the precision of the spiking times, these models involve fast time scales. Certainly, networks built from spiking neuron models can also provide measures of the population activity by averaging across neurons. Thus, patterns of collective activity can be observed in the synaptic current, evoked by the summation of multiple spikes on the target neurons.

To analyse the evolution of the model introduced in §2, we consider two different dynamical variables corresponding to each of the two scales. The NM model activity is given by *y*(*t*) = *y*_{1}(*t*) − *y*_{2}(*t*), where *y*_{1}(*t*) is the EPSP and *y*_{2}(*t*) the IPSP acting upon the pyramidal population (see equation (2.8)). The NN activity is quantified in terms of the LFP as defined in equation (2.13). Both types of model operate in an oscillatory regime. The NM model is an intrinsic oscillator whose frequency can be varied by changing the parameters *B* and *b* (see blue and green lines (dark grey and grey in print, respectively) in figure 2*a,b*). On the other hand, the oscillations of the NN are an emergent property of the system, reflecting the variability of the individual PSPs (i.e. the microscopic events). Hence, their frequency is less well defined (see turquoise line (light grey in print) in figure 2*a*,*b* and power spectra in figure 3*c*,*d*).

Our aim here is to find fingerprints of an effective interaction of scales. To do so, we have studied how the two NM models, one oscillating in the theta band and the other in the alpha band, synchronize their dynamics when the coupling is mediated by the NN (figure 1). The interaction mechanism is bidirectional. This architecture was used by Vicente *et al*. [27] and Gollo *et al*. [28] to demonstrate the emergence of zero-lag synchronization mediated by dynamical relaying between NN populations. In our case, the output of each NM is converted into a firing rate (see equation (2.9)) impinging on a subpopulation of 2000 neurons within the NN. In turn, the firing rate of these selected neurons contributes to both the excitatory and inhibitory PSP densities that act upon the pyramidal populations of the NMs. We also examine the effect of varying several properties of the subpopulation of neurons of the NN, including its number, involved in the coupling, in order to explore the effect of the structural properties of that network on the scale interaction efficiency.

The effect of the coupling intensity *k* on the dynamics of the interacting populations is shown in figure 2. When the NMs are uncoupled, they oscillate in different dynamical regimes that evolve at different frequencies, around 4.5 and 11 Hz, respectively. One NM oscillates in a spike-like manner, whereas the other oscillates more harmonically (figure 2*a*, compare the blue and green (dark grey and grey in print, respectively) lines). The NN, in turn, exhibits collective oscillations within the gamma range, around 45 Hz. The dynamical evolution for the coupled case, at *k* = 1, is shown in figure 2*b*. In this case, the dynamical regimes of the NMs are similar, and they become frequency locked. We have scanned *k* in order to track the transition to the frequency locked regime as coupling increases. Figure 2*c* shows the increase in the maximum cross covariance between the net PSPs of the two NM models, averaged over 20 trials, when increasing *k*. When the NMs operate at different regimes they hardly synchronize but, for sufficiently high *k*, they increase their synchronization with increasing *k*. The averaged frequency mismatch decreases sharply at *k* ≈ 0.6 (figure 2*d*). According to these results, frequency locking for the two NMs is achieved through an NN that oscillates naturally at a much faster scale.

We have further characterized the effect of the interaction through the power spectrum of the time traces. As can be expected, the power spectrum of the mass models in isolation (figure 3*a*) shows a clear peak at their natural frequencies (4.5 and 11 Hz), while the LFP shows a strong peak around 45 Hz (figure 3*c*) that exceeds the non-zero contribution of the slower frequencies approximately 4 Hz. We have seen that increasing coupling leads to a frequency locking regime between the NMs, which is reflected in their spectral behaviour. For instance, at *k* = 1 the power spectra of the two NMs overlap, with a dominant peak around 4 Hz (figure 3*b*). The local gamma peak of the NN is preserved (figure 3*d*), although the major change in amplitude occurs at smaller rhythms, around the frequency of the NMs. This increase in the NN power at the alpha band is due to the emergence of phase locking between this population and the outer NMs, as shown in figure 3*e*. This phase locking results in a sizable cross-correlation between the activities of the microscopic and mesoscopic populations for intermediate values of the size *N* of the NN subpopulation coupling the two NMs, as depicted in figure 3*f* (the difference between cross-correlations with NM1 and NM2 for small *N* is due to the different intrinsic dynamics of the two mesoscopic populations). The effect of *N* is studied in more detail later. The slower time scale of the NMs cannot follow the faster dynamics of the neural network and average out the gamma rhythm, resulting in a frequency shift towards the slower rhythm, which is also enhanced in the NN.

Since the output of the NN arises from the spiking activity of thousands of neurons, the interaction across models is mainly driven by the average dynamics of the population. Although the modelled LFP evolves in a faster time scale, NM models filter out rapid fluctuations. Therefore, the NMs mainly respond to changes of the mean input coming from the NN modulated by *k*.

The input contribution into the NMs coming from the NN dynamics increases the average excitatory and inhibitory input signal into the pyramidal population (denoted by *p*_{e} and *p*_{i}, respectively, in equations (2.11) and (2.12)). Since increasing the constant input to an NM can lead to changes in the dynamical regime (and thus the frequency) of the oscillator [13], one could argue that the role of the NN dynamics is unnecessary to mediate the synchronization transition observed. However, simulations in which the terms given in equations (2.11) and (2.12) are replaced by the temporal average of the coupling contributions indicate that the NMs are unable to synchronize their phases in these conditions (results not shown). This result shows that the NN dynamics is a key ingredient to achieve not only frequency locking, but also phase locking between the two NMs.

In order to take advantage of the microscopic description of the neural network, we have also varied two main features of its architecture: its clustering (figure 4*a*,*b*) and the size of the area involved in the coupling, determined by the number of neurons projecting onto the NMs (figure 5). In graph theory, networks composed of nodes and edges can be characterized by their clustering coefficient, which quantifies the *connectedness* or local connectivity of the network (i.e. the probability that all nodes that are connected to a given node, are also connected between them). According to the Watts & Strogatz algorithm [29], a pure regular network can be turned into a small-world network, in which few edges separate any two nodes, by rewiring the connections. A rewiring probability parameter, rp, determines the probability of replacing an existing edge by another one chosen randomly. Therefore, a rewiring probability equal to zero implies a regular network, whereas a rewiring probability equal to unity implies a completely random network. By studying these parameters, we are changing the dynamics of the NN and, therefore, its capacity to mediate the interaction between the two NMs.

Figure 4*a*,*b* outlines the dependence of the maximum cross covariance and the frequency mismatch between the two NMs on the coupling strength *k* for different rp values. Note that the case rp = 1 corresponds to the results shown in figure 2*c*,*d*. Networks with higher clustering (rp = 0.2) are less efficient in synchronizing the oscillatory output of the NMs. In this case, larger coupling strengths *k* are needed, with respect to a random network (rp = 1), to reach the frequency locking regime. Thus, the topology of the neural network affects the synchronization between the neural ensembles. Random networks have small path lengths at the expense of low clustering, and thus the average transmission time of the action potentials across the population is decreased. In this situation, synchronization arises for smaller coupling strengths. The result for a regular network, rp = 0 (which is not a realistic situation in the brain because the NN dynamics is lost), is also included.

Besides topology, the intrinsic dynamics of the NN also has an impact on the synchronization of NMs. In our NN model, we can slow down the frequency peak of the LFP by increasing the decay time constant *τ*_{d} of the inhibitory synapses (equation (2.4)), without altering the firing rate of the population (result not shown). If the peak of the NN power spectrum is shifted towards the alpha band, closer to where the NMs oscillate, then the maximum cross covariance is reduced and the frequency mismatch is increased for a given *k* value (figure 4*c* and 4*d*, respectively). Thus, even though the NN is operating closer in frequency to the NMs, and its individual neurons fire at the same rate as when the network operates in the gamma band (resulting in a similar MUA activity, result not shown), the NMs are more difficult to synchronize. In the neural network, the action potentials are transiently synchronized and paced according to the time course of inhibition, leading to a recurrent behaviour that causes the global oscillatory dynamics. Faster rhythms, like gamma, correspond to a better precise timing of the firing, i.e. the action potentials of multiple neurons are tightly bounded in time, which seems to be key for the synchronization of the NMs.

Finally, and as mentioned above, we have also studied how the synchronization of the NMs is affected by the size *N* of the subpopulation of neurons that mediate the coupling between them. In the results presented so far, this subpopulation was formed by *N* = 2000 neurons, randomly chosen from the whole population of 4000 neurons of the NN. We have scanned *N* between 1 and 4000 neurons, the latter case corresponding to all neurons in the NN contributing to the firing rate impinging on the NMs and receiving their input. Figure 5*a*,*b* shows the maximum cross covariance and the frequency mismatch for increasing coupling *k* at varying subpopulation sizes. The interaction between the NMs decreases as *N* decreases, and synchronization is only significant for *N* > 1000. *N* directly affects the strength of the coupling between the NN and the NMs, since this parameter determines the average MUA, i.e. the number of spikes elicited within the subpopulation. Hence, given a coupling strength *k* that enables an efficient interaction of the models, larger values of *N* lead to a lower frequency mismatch (figure 5*c*,*d*).

It is important to note that, although the size of the neural network is kept constant, increasing *N* boosts the coupling term, spreading the input from the NM across a larger population of neurons within the NN. Figure 5*e* shows the LFP power spectrum for increasing values of *N* for *k* = 0.9. Similar to the transition from figure 3*c* (network in isolation) to figure 3*d* (*k* = 1 for *N* = 2000), the major changes produced by the coupling occur at small frequencies, where the synchronization scale is centred, while the gamma rhythm interacts directly with the slower dynamics of the NMs. Decreasing *N* dramatically affects the dynamics of the coupling, which only takes into account the activity of this subpopulation. For sizes below *N* ∼ 1000 the interaction is carried out by the low firing and highly noisy activity of small numbers of neurons, which are unable to synchronize large ensembles.

## 4. Discussion

Modelling the dynamics of the full brain from a purely microscopic scale is computationally unfeasible. Thus, a hybrid description of the brain that encompasses multiple scales is an appealing concept. In that scenario, it would only be necessary to represent microscopically those neuronal populations involved in a particular task, and which are monitored with single-cell resolution. The rest of the brain, while modulating the activity of the population of interest, would not necessarily need to be represented with microscopic detail. Currently, this is accomplished by representing the activity of the rest of the brain by a background noisy activity, but this approach is not useful when the neuronal population of interest feeds back into those other brain regions, thereby modifying the background activity that acts upon the population itself. We consider here one way of facing this situation, based on coupling bidirectionally microscopic and mesoscopic descriptions of neuronal populations. In particular, we use synchronization in order to probe the interaction between the two scales. Specifically, we employ a scheme in which two mesoscopic populations are coupled through a third microscopic network, since the behaviour that can be expected from two coupled NM models is well known [11,12], and can be used as a reference to interpret the coordinated behaviour emerging from our hybrid scenario.

Our results do not imply that two NM oscillators can only synchronize through the mediation of an NN. In fact, if all three neuronal populations were described by NNs (or by NMs, for that matter) synchronization will also arise (see, for instance, the studies of Vicente *et al*. [27] and Gollo *et al*. [28] for the case of three coupled NNs leading to zero-lag synchronization). Neither do we claim that two brain oscillators can synchronize only through the mediation of a third (see, for instance, the study of David & Friston [12] for an example of synchronization between two coupled NMs). What our study shows is that two mesoscopic brain oscillators can synchronize *even when* they are coupled only through a mediating population that is described by a microscopic model. In that sense, we use synchronization as a tool to probe the interaction between different spatial scales of neuronal populations. Previous efforts have been devoted to analysing this interaction by performing a direct comparison of the behaviours of the microscopic and mesoscopic models. Faugeras *et al*. [23], for instance, derived the equations of evolution of NMs from the dynamics of a network of neurons described by a voltage-based model, by performing an involved mean field analysis of the network, an approach that would be very challenging to apply to spiking neuron models. In order to perform such a multiscale mapping, Rodrigues *et al*. [24] had to apply strong assumptions that included high correlation between the neurons in the microscopic populations and low-amplitude input currents. Here, we have attempted to circumvent the complexity of those approaches by using a more phenomenological strategy, whose goal is to test whether microscopic and mesoscopic descriptions of neuronal populations communicate with one another by using synchronization as a proxy of effective communication.

Even when the NN operates in a fast dynamical collective regime in the gamma range, a sufficiently large subpopulation of neurons within that network is able to mediate the communication and subsequent synchronization between two NMs that are described mesoscopically and operate at much lower frequencies. Frequency and phase locking arise even when the two NMs operate at very different frequencies (in the theta and alpha bands) and with very different dynamical features (spike-like dynamics in one case and quasi-harmonic dynamics in the other). Structural clustering within the NN diminishes the ability of the microscopic neuronal population to induce synchronization. The size of the subpopulation of neurons that directly coupled the two NMs must also be large enough to allow the intrinsically irregular neurons to reach a sufficiently strong collective regime through which the two NMs can communicate.

Two main features indicate the non-trivial contribution of the microscopic NN in mediating the synchronization between the mesoscopic models. First, the two mesoscopic populations lock not only in frequency, but also in phase, when they interact with a dynamically evolving NN. If the role of the network is played by an increased constant input to the NMs equal to the average activity of the NN, phase locking disappears. Second, if the NN is made to operate in a slower collective regime (e.g. in the alpha band) the synchronization between the NMs is decreased (while still being significant), even though the three oscillators are now closer in frequency.

The synchronization between the NMs is mediated by the locking between the NMs and the NN, which leads to an increase in the alpha band activity of the NN, as reflected in figure 3. The fact that synchronization is maintained even when the NN is operating in the alpha band (figure 4*c*,*d*) indicates that the intrinsic NN dynamics does not interfere notably in the communication between the NM populations. Furthermore, the fact that synchronization improves slightly when the NN is operating in gamma (as shown also in figure 4*c*,*d*) shows that fast and slow scales interact to a certain extent in order to drive the synchronization. We interpret this to be due to an increase in the precise timing of the firing that is associated with a faster neuronal rhythm. The results reported here point towards an alternative way to probe the interaction of scales in the brain, by using synchronization between neuronal populations as a way of testing the structural and functional conditions under which scale interaction occurs.

## Funding statement

This work was supported by the European Commission through the FP7 Marie Curie Initial Training Network 289146 (NETT: Neural Engineering Transformative Technologies), the Ministerio de Economia y Competividad (Spain, project FIS2012-37655) and the Generalitat de Catalunya (project 2009SGR1168). J.G.O. acknowledges support from the ICREA Academia programme.

## Footnotes

One contribution of 12 to a Theme Issue ‘Complex network theory and the brain’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.