The cortical mass model

We consider a cortical mass model which consists of two populations of excitatory adaptive (E) and inhibitory (I) exponential integrate-and-fire (AdEx) neurons (Fig 1). Both populations are delay-coupled and the excitatory population has a somatic adaptation feedback mechanism. The low-dimensional mean-field model (Fig 1a) is derived from a large network of spiking AdEx neurons (Fig 1b).

For the construction of the mean-field model, a set of conditions need to be fulfilled: We assume the number of neurons to be very large, all neurons within a population to have equal properties, and the connectivity between neurons to be sparse and random. Additional assumptions about the mathematical nature and a detailed derivation of the mean-field model is presented in the Methods section.

Bifurcation diagrams: A map of the dynamical landscape

The E-I motif shown in Fig 1 can occupy various network states, depending on the baseline inputs to both populations. By gradually changing the inputs, we map out the state space of the system, depicted in the bifurcation diagrams in Fig 2. Small changes of the parameters of a nonlinear system can cause sudden and dramatic changes of its overall behavior, called bifurcations. Bifurcations separate the state space into distinct regions of network states between which the system can transition from one to another. In our case, the dynamical state of the E-I system depends on external inputs to both subpopulations, which are directly affected by external electrical stimulation and other driving sources, e.g. inputs from other neural populations such as other brain regions.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Bifurcation diagrams and time series.

Bifurcation diagrams depict the state space of the E-I system in terms of the mean external input currents to both subpopulations α ∈ {E, I}. (a) Bifurcation diagram of mean-field model without adaptation with up and down-states, a bistable region bi (green dashed contour) and an oscillatory region LC_EI (white solid contour). (b) Diagram of the corresponding AdEx network with N = 50 × 10³ neurons. (c) Mean-field model with somatic adaptation. The bistable region is replaced by a slow oscillatory region LC_aE. (d) Diagram of the corresponding AdEx network. The color in panels a—d indicate the maximum population rate of the excitatory population (clipped at 80 Hz). (e) Example time series of the population rates of excitatory (red) and inhibitory (blue) populations at point A2 (top row) which is located in the fast excitatory-inhibitory limit cycle LC_EI, and at point B3 (bottom row) which is located in the slow limit cycle LC_aE. (f) Time series at corresponding points for the AdEx network. All parameters are listed in Table 1. The mean input currents to the points of interest A1-A3 and B3-B4 are provided in Table 2.

https://doi.org/10.1371/journal.pcbi.1007822.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Summary of the model parameters.

Parameters apply for the Mean-Field model and the spiking AdEx network.

https://doi.org/10.1371/journal.pcbi.1007822.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Values of the mean external inputs to the excitatory () and the inhibitory population () in units of nA for points of interest in the bifurcation diagrams Fig 2.

https://doi.org/10.1371/journal.pcbi.1007822.t002

Comparing the bifurcation diagrams of the mean-field model (Fig 2a and 2c) to the ground truth spiking AdEx network (Fig 2b and 2d) demonstrates the similarity between both dynamical landscapes. Transitions between states take place at comparable baseline input values and in a well-preserved order.

Since the space of possible biophysical parameter configurations is vast, we focus on two variants of the model: one without a somatic adaptation mechanism, Fig 2a and 2b, and one with finite sub-threshold and spike-triggered adaptation in Fig 2c and 2d. Both variants feature distinct states and dynamics.

Somatic adaptation causes slow oscillations.

In Fig 2c and 2d, bifurcation diagrams of the system with somatic adaptation are shown. Compared to Fig 2a and 2b (without adaptation), the state space, including the oscillatory region LC_EI, is shifted to the right, meaning that larger excitatory input currents are necessary to compensate for the inhibiting sub-threshold adaptation currents. The most notable effect that is caused by adaptation is the appearance of a slow oscillatory region labeled LC_aE in Fig 2c and 2d. The reason for the emergence of this oscillation is the destabilizing effect the inhibiting adaptation currents have on the up-state inside the bistable region [29–31]. As the mean adaptation current builds up due to a high population firing rate, the up-state “decays” and the system transitions to the down-state. The resulting low activity causes a decrease of the adaptation currents which in turn allow the activity to increase back to the up-state, resulting in a slow oscillation. These low-frequency oscillations range from 0.5 Hz to 5 Hz for the parameters given.

The bifurcation diagrams in Fig 3 show how the emergence of the slow oscillation depends on the adaptation mechanism. Increasing the subthreshold adaptation parameter primarily shifts the state space to the right whereas a larger spike-triggered adaptation parameter value enlarges oscillatory regions. Both parameters cause the bistable region to shrink until it is eventually replaced by a slow oscillatory region LC_aE. Again, the state space of the mean-field model (Fig 3a) reflects the AdEx network (Fig 3b) accurately.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Transition from multistability to slow oscillation is caused by somatic adaptation.

Bifurcation diagrams depending on the external input currents to both populations α ∈ {E, I} for varying somatic adaptation parameters a and b. The color indicates the maximum rate of the excitatory population. Oscillatory regions have a white contour, bistable regions have a green dashed contour. (a) Bifurcation diagrams of the mean-field model. On the diagonal (bright-colored diagrams), adaptation parameters coincide with (b). (b) Bifurcation diagrams of the corresponding AdEx network, N = 20 × 10³. All parameters are listed in Table 1.

https://doi.org/10.1371/journal.pcbi.1007822.g003

Time-varying stimulation and electric field effects

To describe the time-dependent properties of the system, we study the effects of time-varying external stimulation and the interactions with ongoing oscillatory states. External stimulation is implemented by coupling an electric input current to the excitatory population. This additional input current may be a result of an externally applied electric field or synaptic input from other neural populations. For the cases without adaptation, we can calculate an equivalent extracellular electric field strength that correspond to the effects of an input current (see Methods).

Since due to the presence of apical dendrites, excitatory neurons in the neocortex are most susceptible to electric fields [32], only time-varying input to the excitatory population is considered. This choice is also motivated in the context of inter-areal brain network models where connections between brain areas are usually considered between excitatory subpopulations.

Given the multitude of possible states of the system, its response to external input critically depends on the dynamical landscape around its current state. It is important to keep in mind that the bifurcation diagrams (Figs 2 and 3) are valid only for constant external input currents. However, they provide a helpful estimation of the dynamics of the non-stationary system assuming that the bifurcation diagrams do not change too much as we vary the input parameter over time.

Fig 4a–4f show how a step current input pushes the system in and out of specific states of the E-I system. A positive step current represents a movement in the positive direction of the -axis in Fig 2. Fig 4a and 4b show input-driven transitions from the low-activity down-state to the fast oscillatory limit cycle LC_EI. Similar behavior can be observed in Fig 4c and 4d where we push the system’s state from LC_EI to the up-state, effectively being able to turn oscillations on and off with a direct input current.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. State-dependent population response to time-varying input currents.

Population rates of the excitatory population (black) with an additional external electrical stimulus (red) applied to the excitatory population. (a, b) A DC step input with an amplitude of 60 pA (equivalent E-field amplitude: 12 V/m) pushes the system from the low-activity fixed point into the fast limit cycle LC_EI. (c, d) A step input with amplitude 40 pA (8 V/m) pushes the system from LC_EI into the up-state. (e, f) In the multistable region bi, a step input with amplitude 100 pA (20 V/m) pushes the system from the down-state into the up-state and back. (g, h) Inside the slow oscillatory region LC_aE, an oscillating input current with amplitude 40 pA and a (matched) frequency of 3 phase-locks the ongoing oscillation. (i, j) A slow 4 Hz oscillatory input with amplitude 40 pA drives oscillations if the system is close to the oscillatory region LC_aE. For the AdEx network, N = 100 × 10³. All parameters are given in Table 1. The parameters of the points of interest are given in Table 2.

https://doi.org/10.1371/journal.pcbi.1007822.g004

Inside the bistable region, we can use the hysteresis effect to transition between the down-state and the up-state and vice versa. After application of an initial push in the desired direction, the system remains in that state, reflecting the system’s bistable nature.

With adaptation turned on, a slow oscillatory input current can entrain the ongoing oscillation. In Fig 4g and 4h, the oscillation is initially out of phase with the external input but is quickly phase-locked. Placed close to the boundary of the slow oscillatory region LC_aE, we show in Fig 4i and 4j how an oscillatory input with a similar frequency as the limit cycle periodically drives the system from the down-state into one oscillation period.

Close inspection of Fig 4b shows that the fast oscillation in the AdEx network has a varying amplitude, in contrast to the mean-field model Fig 4a. This difference is due to noise resulting from the finite size of the AdEx network and decreases as the network size N increases. (see S8 Fig) All of the state transitions take longer for the AdEx network, as it is visible in Fig 4d for example. Additionally, transitions to the up-state and the slow limit cycle LC_aE are accompanied by transient ringing activity (Fig 4f and 4h), which is not well-captured by the mean-field model.

Frequency entrainment with oscillatory input.

To study the frequency-dependent response of the E-I system, we vary the amplitude and frequency of an oscillatory input to the excitatory population (Fig 5). The unperturbed system is parameterized to be in the fast limit cycle LC_EI with its endogenous frequency f₀.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Frequency entrainment of the population activity in response to oscillatory input.

The color represents the log-normalized power of the excitatory population’s rate frequency spectrum with high power in bright yellow and low power in dark purple. (a) Spectrum of the mean-field model parameterized at point A2 with an ongoing oscillation frequency of f₀ = 22 Hz (horizontal green dashed line) in response to a stimulus with increasing frequency and an amplitude of 20 pA. An external electric field with a resonant stimulation frequency of f₀ has an equivalent strength of 1.5 V/m. The stimulus entrains the oscillation from 18 Hz to 26 Hz, represented by a dashed green diagonal line. At 27 Hz, the oscillation falls back to its original frequency f₀. At a stimulation frequency of 2f₀, the ongoing oscillation at f₀ locks again to the stimulus in a smaller range from 43 Hz to 47 Hz. (b) AdEx network with f₀ = 30 Hz. Entrainment with an input current of 40 pA is effective from 27 Hz to 33 Hz. Electric field amplitude with frequency f₀ corresponds to 2.5 V/m. (c) Mean-field model with a stimulus amplitude of 100 pA (7.5 V/m at 22 Hz). Green dashed lines mark the driving frequency f_ext and its first and second harmonics f_1H and f_2H and subharmonics f_1SH and f_2SH. Entrainment is effective from the lowest stimulation frequency up to 36 Hz at which the oscillation falls back to a frequency of 20 Hz. New diagonal lines appear due to interactions of the endogenous oscillation with the entrained harmonics and subharmonics. (d) AdEx network with stimulation amplitude of 140 pA (8.75 V/m at 30 Hz). For the AdEx network, N = 20 × 10³. All parameters are given in Tables 1 and 2.

https://doi.org/10.1371/journal.pcbi.1007822.g005

The external stimulus with frequency f_ext entrains the ongoing oscillation in a range around f₀, the resonant frequency of the system. Here, the ongoing oscillation essentially follows the external drive and adjusts its frequency to it (Fig 5a). A second (narrower) range of frequency entrainment appears as f_ext approaches 2f₀, representing the ability of the input to entrain oscillations at half of its frequency. Due to interference of the frequencies of ongoing and external oscillations, the spectrum has peaks at the difference of both frequencies which appear as X-shaped patterns in the frequency diagrams. The AdEx network shows a similar behavior (Fig 5b), albeit the range of entrainment is smaller than in the mean-field model, despite the stimulation amplitude being twice as large.

For stronger oscillatory input currents, the range of frequency entrainment is widened considerably. In Fig 5c and 5d, the input dominates the spectrum at very low frequencies. The peak of the spectrum reverts back to approximately f₀ if the external frequency f_ext is close to the first harmonic 2f₀ of the endogenous frequency. We see multiple lines emerging in the frequency spectra which correspond to the harmonics and subharmonics of the external frequency and its interaction with the endogenous frequency f₀, creating complex patterns in the diagrams. Differences between the spectrograms of the AdEx network in Fig 5d and the mean-field model Fig 5c can be largely attributed to the fact that the AdEx network consistently needs stronger inputs to obtain the same effect as in the mean-field model. This results in horizontal lines in areas where frequency entrainment is not effective and in faint and short diagonal lines between the lines that represent the (sub-)harmonics which are caused by interactions with the (sub-)harmonics. In the mean-field model, we mainly observe clear diagonal lines, indicating successful entrainment. Another source for the differences is the inherently noisy dynamics of the AdEx network, due to its finite size (see S8 Fig).

Overall, there is a good qualitative agreement of the frequency spectra of both models, reflecting that interactions of time-varying external inputs and ongoing oscillations are well-represented by the mean-field model.

Phase locking with oscillatory input.

Here we quantify the ability of an oscillating external input current to the excitatory population to synchronize an ongoing neural oscillation to itself if both frequencies, the driver and the endogenous frequency, are close to each other (frequency matching). An example time series of a stimulus entraining an ongoing slow oscillation is shown in Fig 4h.

In Fig 6, we find phase locking by measuring the time course of the phase difference between the stimulus and the population rate. If phase locking is successful, the phase difference remains constant. In Fig 6a, the region of phase locking for an external input of frequency f_ext is centered around the endogenous frequency f₀ of the unperturbed system. Increasing the stimulus amplitude widens the range around f₀ at which phase locking is effective, producing Arnold tongues in the diagram. An example time series of successful phase locking inside this region is shown in Fig 6c at point 1. If the input is not able to phase-lock the ongoing activity, a small difference between the driver frequency f_ext and f₀ can cause a slow beating of the activity with a frequency of roughly the difference |f_ext − f₀|. Thus, a small frequency mismatch can produce a very slowly oscillating activity (Fig 6c at points 2-4). Fig 6d at point 2 shows the same drifting effect in the AdEx network. Due to finite-size noise in the AdEx network, an irregular switching between synchrony and asynchrony can be observed at the edges of the phase locking region in Fig 6d at point 3. Compared to the mean-field model, the frequency of the beating activity in the AdEx network is less regular (Fig 6d at point 4).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Phase locking of ongoing oscillations via weak oscillatory inputs.

The left panels show heatmaps of the level of phase locking for (a) the mean-field model and (b) the AdEx network for different stimulus frequencies and amplitudes. Dark areas represent effective phase locking and bright yellow areas represent no phase locking. Phase locking is measured by the standard deviation of the Kuramoto order parameter R(t) which is a measure for phase synchrony. White dashed lines correspond to electric fields with equivalent strengths in V/m. (c) Time series of four points indicated in (a) with the excitatory population’s rate in black and the external input in red (upper panels). In the lower panels, the Kuramoto order parameter R(t) is shown, measuring the phase synchrony between the population rate and the external input. Constant R(t) represents effective phase locking (phase difference between rate and input is constant), fluctuating R(t) indicates dephasing of both signals, hence no phase locking. (d) Corresponding time series of points in (b). Both models are parameterized to be in point A2 inside the fast oscillatory region LC_EI. Insets show zoomed-in traces from 15 to 16 seconds. For the AdEx network, N = 20 × 10³. All parameters are given in Table 1.

https://doi.org/10.1371/journal.pcbi.1007822.g006

In the phase locking diagrams Fig 6a and 6b, the equivalent external electric field amplitudes are shown. Small amplitudes (0.2 V/m for the mean-field model, 0.5 V/m for the AdEx network) are able to phase-lock the ongoing oscillations if the frequencies roughly match.

Somatic adaptation causes slow network oscillations

The AdEx neuron model allows for incorporation of a slow potassium-mediated adaptation current, typically found in cortical pyramidal neurons [33]. Due to somatic adaptation, in the bistable region, the up-state loses its stability. The bistable region transforms into a second oscillatory regime (Fig 3) in which the population activity oscillates at low frequencies between 0.5 Hz and 5 Hz. This oscillatory region coexists with the fast excitatory-inhibitory oscillation. Other computational studies have focused on the origin of this adaptation-mediated oscillation [29–31], the interaction of adaptation with noise-induced state switching between up- and down-state [29, 34, 35] and how adaptation affects the intrinsic timescales of the network [36, 37].

Electric field effects and relation to experimental observations

Using the bifurcation diagrams (Fig 2), we mapped out several points of interest that represent different network states. The type of reaction to external stimulation depends on the current state of the system, as seen in the population time series during stimulation in Fig 4. Close to edges of attractors, direct currents can cause bifurcations and trigger a sudden change of the dynamics, such as transitions from a low activity down-state to a state with oscillatory activity.

In vitro stimulation experiments with electric fields [3] have shown that (time-constant) direct fields are able to switch on and off oscillations at field strengths of 6 V/m. In Fig 4, we could observe this at 8 − 12 V/m, when the system if placed close to the oscillatory state. This difference in amplitudes can be attributed to the chosen initial state of the system and could be reduced if the background inputs were parameterized closer to the limit cycle.

Inside oscillatory regions, oscillatory input causes phase locking and frequency entrainment. Frequency entrainment is the ability of an external stimulus to force the endogenous oscillation to follow its frequency. To study how frequency entrainment depends on the frequency and amplitude of the stimulus, we analyzed frequency spectrograms of the population activity when subject to external oscillating stimuli with increasing frequencies (Fig 5). We observed shifts of the peak frequency around the endogenous frequency at amplitudes corresponding to field strengths of 1.5 V/m in the mean-field model and 2.5 V/m in the AdEx network. Similar effects have been reported in in vitro experiments, where the frequency of the ongoing oscillations changed along with the stimulus frequency [3, 38–40].

Interestingly, the field amplitudes at which frequency entrainment is visible are on the same order of endogenous fields in the brain, generated by the neural activity itself. Electrophysiological experiments show that these fields can be as strong as 3.5 V/m [11] and that ephaptic coupling plays a significant role in the brain [12]. Our findings support this observation from a theoretical perspective, since one of our main results is that considerable effects on the population dynamics are expected at these field strengths.

We also observed frequency entrainment of the subharmonics of the endogenous oscillation as it was shown in in vitro experiments in Refs. [3, 39] and its harmonics in Ref. [38]. This effect could be valuable for experimental conditions where it is impractical to use stimulation frequencies close to the endogenous frequency of ongoing oscillations in the studied neural system. The range of frequency entrainment around the natural frequency of the endogenous oscillation widens as the stimulus amplitude increases, which was also observed in similar computational studies [41, 42].

If the stimulus frequency is close to the endogenous frequency (frequency matching), an oscillatory stimulus can force the ongoing oscillation to synchronize its phase with the stimulus, known as phase locking, phase entrainment, or coherence. Phase locking of ongoing brain activity to a stimulus has been observed in multiple noninvasive brain stimulation studies, including Refs. [43–45], and it has been shown to affect information processing properties of the brain [7, 8] as particularly sensory information processing depends on phase coherence of oscillations between distant brain regions [46, 47]. Compared to frequency entrainment, very weak input currents are able to phase lock ongoing oscillations. In agreement with these experiments, we find phase locking to be effective at electric field strengths of around 0.5 V/m (Fig 6), which is in the range of typical field strengths generated by transcranial alternating current stimulation (tACS) [48].

To summarize: Our results confirm the interesting notion that, while weak electric fields with strengths in the order of 1 V/m that are typically applied in tACS experiments have only a small effect on the membrane potential of a single neuron [32], the effects on the network however, and therefore on the dynamics of the population as a whole, can be quite significant, which was also observed experimentally [49]. This indicates that field effects are strongly amplified in the network. Considering slightly stronger fields, our results suggest that endogenous fields, generated by the activity of the brain itself, are expected to have a considerable effect on neural oscillations, facilitating phase and frequency synchronization across neighboring cortical brain areas.

Validity of the mean-field method and limitations of our approach

We found all observed input-dependent effects in the AdEx network to be well-represented by the mean-field model, which demonstrates its accuracy also in the non-stationary case. However, partly due to the difference of parameters that define the states in the bifurcation diagrams, the AdEx network consistently requires larger input amplitudes in order to cause the same effect size as observed in the mean-field model (Figs 5 and 6). Although it was not investigated here, we hypothesize that the number of neurons might play an important role in how much external inputs are amplified within a network.

Related to this is the fact that our approximation assumes the number of neurons to be infinitely large. Therefore, differences between the mean-field model and the AdEx network in the case without external stimulation also depend on the network size. However, we find good agreement between the bifurcation diagrams for as low as N = 4 × 10³ neurons (S7 Fig). With increasing network size, the amplitudes of oscillations in the AdEx network shown in Fig 4b approach the predictions of the mean-field model and become less irregular (S8 Fig).

Comparing the bifurcation diagrams of both models (Fig 2), the shape of the oscillatory region as well as the frequencies of the oscillations differ (see S4 Fig). We suspect that the oscillatory states are where the steady-state approximations that are used to construct the mean-field model break down due to the fast temporal dynamics in this state. Hence, both models have notable differences between the oscillatory regions.

Sharp transitions between states cause transient effects that are visible as ringing oscillations in the population firing rate, typically observed in simulations of spiking networks (such as in Fig 4f and 4h) or experimentally [50]. The poor reproduction of these oscillations in our model is likely due to the use of an exponentially-decaying linear response function instead of a damped oscillation which would be a better approximation of the true response (c.f. Ref. [23]). This constitutes a possible improvement of our work. Recent advancements in mean-field models of cortical networks [51] can account for its finite size as well as reproduce transient oscillations caused by sharp input onsets.

Another important limitation of our method is the assumption of homogeneous and weak synaptic coupling in the mathematical derivation of our mean-field model. Synapses in the brain are known to be log-normally distributed [52] with long tails, implying the existence of few but strong synapses. Other computational papers have specifically focused on the effect of strong synapses on the population activity (cf. [53] and [54]). Therein, the incorporation of strong synapses causes the emergence of a new asynchronous state in which the firing rates of individual neurons fluctuate strongly, similar to chaotic states studied in networks of rate neurons [55], which are qualitatively different from the up-state that we observed (see S1 Fig). In Ref. [53], the author shows that firing rate models similar to what we consider break down and cannot capture the large fluctuations present in this state. We therefore conclude that our mean-field model is limited to describing only weak synaptic coupling.

Furthermore, a number of assumptions were made in our model of electric field interaction. Most importantly, we have chosen typical morphological and electrophysiological parameters of a ball-and-stick neuron model to represent pyramidal cells in layer 4/5 of cortex (see Methods). Despite the simplicity of the ball-and-stick model, it was shown in Ref. [28] that it can reproduce the somatic polarization of a pyramidal cell in a weak and uniform field. This was then translated into effective input currents to point neurons which lack any morphological features. In addition to the crude assumption that all neurons have the same simple morphology (effects of a more complex morphology were studied in Ref. [56]), we also assumed perfect alignment of the dendritic cable to the external electric field. While the latter might be a good approximation for a local region of the cortex it is not the case for the brain as a whole with its folded structure. It has been shown that the somatic polarization of a neuron strongly depends on the angle between the neuron’s main axis and the electric field [32].

We only focused on field effects that are caused by the dendrite. Although we expect that, in principle, axons could contribute to the somatic polarization in a weak and uniform electric field, their contribution could be relatively small, since most cortical axons are not geometrically aligned with each other the way that dendrites are organized in the columnar structure of the cortex [57].

Finally, we assumed that the field effects are only subthreshold such that our results do not generalize to stimulation scenarios with strong electric fields that can elicit action potentials by themselves.

Conclusion

Overall, our observations confirm that a sophisticated mean-field model of a neural mass is appropriate for studying the macroscopic dynamics of large populations of spiking neurons consisting of excitatory and inhibitory units. To our knowledge, such a remarkable equivalence of the dynamical states between a mean-field neural mass model and its ground-truth spiking network model has not been demonstrated before. Our analysis shows that mean-field models are useful for quickly exploring the parameter space in order to predict states and parameters of the neural network they are derived from. Since the dynamical landscapes of both models are very similar, we believe that it should be possible to reproduce a variety of stationary and time-dependent properties of large-scale network simulations using low-dimensional population models. This may help to mechanistically describe the rich and plentiful observations in real neural systems when subject to stimulation with electric currents or electric fields, such as switching between bistable up and down-state or phase locking and frequency entrainment of the population activity.

Bifurcations, as studied in dynamical systems theory, offer a plausible mechanism of how networks of neurons as well as the brain as a whole [58, 59] can change its mode of operation. Understanding the state space of real neural systems could be beneficial for developing electrical stimulation techniques and protocols, represented as trajectories in the dynamical landscape, which could be used to reach desirable states or specifically inhibit pathological dynamics.

Due to the variety of possible macroscopic network states that arise from this basic E-I architecture, it is critical to consider the state of the system in order to comprehend and predict its response to external stimuli. This might explain the numerous seemingly inconclusive experimental results from noninvasive brain stimulation studies [5, 6] where it is hard to account for the state of the brain before stimulation. In conclusion, additional to the stimulus parameters, the response of a system to external stimuli has to be understood in context of the dynamical state of the unperturbed system [3, 42, 60].

Neural population setting

In order to derive the mean-field description of an AdEx network, we consider a very large number of N → ∞ neurons for each of the two populations E and I. We assume (1) random connectivity (within and between populations), (2) sparse connectivity [61, 62], but each neuron having a large number of inputs [63] K with 1 ≪ K ≪ N, (3) and that each neuron’s input can be approximated by a Poisson spike train [64, 65] where each incoming spike causes a small (c/J ≪ 1) and quasi-continuous change of the postsynaptic potential (PSP) [66] (diffusion approximation).

The spiking neuron model

The adaptive exponential (AdEx) integrate-and-fire neuron model forms the basis for the derivation of the mean-field equations as well as the spiking network simulations. Each population α ∈ {E, I} has N_α neurons which are indexed with i ∈ [1, N_α]. The membrane voltage of neuron i in population α is governed by (1) (2) The first term of I_ion (Eq 2) describes the voltage-depended leak current, the second term the nonlinear spike initiation mechanism, and the last term I_A, the somatic adaptation current. I_i,ext(t) = μ^ext(t) + σ^ext ξ_i(t) is a noisy external input. It consists of a mean current μ^ext(t) which is equal across all neurons of a population and independent Gaussian fluctuations ξ_i(t) with standard deviation σ^ext (σ^ext is equal for all neurons of a population). For a neuron in population α, synaptic activity induces a postsynaptic current I_i which is a sum of excitatory and inhibitory contributions: (3) with C being the membrane capacitance and J_αβ the coupling strength from population β to α, representing the maximum current when all synapses are active. The synaptic dynamics is given by (4) s_i,αβ(t) represents the fraction of active synapses from population β to α and is bound between 0 and 1. G_ij is a random binary connectivity matrix with a constant row sum K_α and connects neurons j of population β to neurons i of population α. With the constraint of a constant in-degree K_α of each unit, all neurons of population β project to neurons of population α with a probability of p_αβ = K_α/N_β and α, β ∈ {E, I}. G_ij is generated independently for every simulation. The first term in Eq 4 is an exponential decay of the synaptic activity, whereas the second term integrates all incoming spikes as long as s_i,αβ < 1 (i.e. some synapses are still available). The first sum is the sum over all afferent neurons j, and the second sum is the sum over all incoming spikes k from neuron j emitted at time t^k after a delay d_αβ. If s_i,αβ = 0, the amplitude of the postsynaptic current is exactly C ⋅ c_αβ which we set to physiological values from in vitro measurements [67] (see Table 1).

For neurons i of the excitatory population, the adaptation current I_A,i(t) is given by (5) a representing the subthreshold adaptation and b the spike-dependent adaptation parameters. The inhibitory population doesn’t have an adaptation mechanism, which is equivalent to setting these parameters to 0. When the membrane voltage crosses the spiking threshold, V_i ≥ V_s, the voltage is reset, V_i ← V_r, clamped for a refractory time T_ref, and the spike-triggered adaptation increment is added to the adaptation current, I_A,i ← I_A,i + b. All parameters are given in Table 1.

Finally, we define the mean firing rate of neurons in population α as (6) which measures the number of spikes in a time window dt, set to the integration step size in our numerical simulations.

The mean-field neural mass model

For a sparsely connected random network of AdEx neurons as defined by Eqs 1–5, the distribution of membrane potentials p(V) and the mean population firing rate r can be calculated using the Fokker-Planck equation in the thermodynamic limit N → ∞ [13, 68]. Determining the distribution involves solving a partial differential equation, which is computationally demanding. A low-dimensional linear-nonlinear cascade model [24, 69] can be used to capture the steady-state and transient dynamics of a population in form of a set of simple ODEs. Briefly, for a given mean membrane current μ_α with standard deviation σ_α, the mean of the membrane potentials as well as the population firing rate r_α in the steady-state can be calculated from the Fokker-Plank equation [70] and captured by a set of simple nonlinear transfer functions Φ(μ_α, σ_α) (shown in Fig 7a and 7b).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Precomputed quantities of the linear-nonlinear cascade model.

(a) Nonlinear transfer function Φ for the mean population rate (Eq 15) (b) Transfer function for the mean membrane voltage (Eq 10) (c) Time constant τ_α of the linear filter that approximates the linear rate response function of AdEx neurons (Eq 7). The color scale represents the level of the input current variance σ_α across the population. All neuronal parameters are given in Table 1.

https://doi.org/10.1371/journal.pcbi.1007822.g007

These transfer functions can be precomputed (once) for a specific set of single AdEx neuron parameters. All other parameters, such as input currents, network parameters and synaptic coupling strengths, as well as the parameters that govern the somatic adaptation mechanism, the membrane timescale, and the synaptic timescale are identical and directly represented in the equations of the mean-field model. Thus, for any given parameter configuration of the AdEx network, there is a direct translation to the parameters that define the mean-field model. This also allows for direct comparison of both models under changes of said parameters.

The reproduction accuracy of the linear-nonlinear cascade model for a single population has been systematically reviewed in Ref. [23] and has been shown to reproduce the dynamics of an AdEx network in a range of different input regimes quite successfully, while offering significant increase in computational efficiency.

Obtaining bifurcation diagrams and determining bistability

Each point in the bifurcation diagrams in Figs 2 and 3 was simulated for pairs of external inputs and and the resulting time series of the excitatory population rate of the mean-field model and the AdEx network were analyzed and the dynamical state was classified.

To classify a point in the state space as bistable or not, in both models, we apply a negative and a subsequent positive stimulus to the excitatory population and measure the difference in activity after both stimuli are turned off again. In the AdEx network, a simple step input can cause over- and under-shoot as a reaction, which is a problem when assessing the stability of a basin of attraction around a fixed point. To overcome this problem, we constructed a slowly-decaying stimulus (in contrast to previous work [80], where bistability was identified using a step current). An inverted example of this stimulus is shown in Fig 4e and 4f. Using this stimulus, we first made sure that the population rate is in the down-state (the initial state) with an initial negative external input current that slowly decays back to zero. We then kicked the activity into the up-state (the target state) with a positive input and then let the current slowly decay to zero again. A slowly-decaying stimulus (in contrast to a step stimulus) ensures that transient effects such as over- and undershooting are minimized that would otherwise disturb the target state. As a result, the stability of the target state can be observed. We determined whether the up-state is stable or the activity has decayed back into the down-state by comparing the 1 s mean of the population rates after both stimuli have decayed. We classified a state as bistable if the rate difference after both kicks and subsequent relaxation phases was greater than 10 Hz. This threshold value was chosen to be smaller than every observed difference between the up- and the down-state. We confirmed the validity of this method for the mean-field model by using a continuation method to determine the stability of the fixed-point states, which provided the same bifurcation diagrams.

Determining frequency spectra of the population activity

In the bifurcation diagrams in Figs 2 and 3, regions were classified as oscillating if the time series showed oscillations during the last 1 s after the first (negative) stimulus pushed the system into the down-state. The power spectrum of this oscillation was computed using the implementation of Welch’s method [81] scipy.signal.welch in the Python package SciPy (1.2.1) [82]. A rolling Hanning window of length 0.5 s was used to compute the spectrum. If the dominant frequency was above 0.1 Hz and its power density was above 1 Hz we classified the state as oscillating. Visual observation of the time series confirmed that these thresholds classified the oscillating regions well. In cases where the transient of 1 s was too short such that the activity state of the population jumped from the down-state to the up-state within this period, misclassifications of these points as oscillatory states caused artifacts at the right-hand border of the bistable region to the up-state.

In Fig 5, we determined frequency entrainment by observing changes of the frequency spectrum of the population activity r_E. Each run was simulated for 6 s. We waited for 1 s for transient effects to vanish before turning on the oscillating stimulus and measured the power spectrum of the remaining 5 s. A rolling Hanning window of length 1 s was used. For better visibility, the power was normalized between 0 and 1 on a logarithmic scale and plotted with a linear colormap.

Measuring phase locking using the Kuramoto order parameter

In Fig 6 we quantified the degree of phase locking of an oscillatory input current with the E-I system’s ongoing oscillation. We calculated the Kuramoto order parameter [83] to measure phase synchrony. The Kuramoto order parameter R is given by: (22) In our case, the number of oscillators N_osc = 2 and Φ_j ∈ [0, 2π) is the instantaneous phase of the stimulus (j = 1) and the population activity r_E (j = 2). We define the instantaneous phase Φ_j at time t as (23) where t_n is the time of the last maximum and t_n−1 the penultimate one. To robustly detect the oscillation maxima of the noisy AdEx network population rate, the time series was first smoothed using the Gaussian filter scipy.ndimage.filters.gaussian_filter implemented in SciPy. The Gaussian kernel had a standard deviation of 5 ms. Then, the maxima were detected using the peak finding algorithm scipy.signal.find_peaks_cwt with a peak-width between 0.1 and 0.2 ms.

For R = 1, perfect (zero-lag) phase synchronization is reached; if R ≈ 0, the oscillations are maximally desynchronized. To measure phase locking in Fig 6a and 6c, we calculated the standard deviation in time of R(t) after transient effects vanished for t > 1.5 s. A low standard deviation means that the phase difference between the input and the ongoing oscillation stays constant.

Calculating equivalent electric field strengths

Our results can be used to estimate the necessary amplitude of an external electric field to reproduce the effects of electrical input currents. An external field at the location of a neural population might be produced by endogenous electric fields due to the activity of a neural population or external stimulation techniques such as transcranial electrical stimulation (tES) with direct (tDCS) or alternating (tACS) currents. The lack of a spatial extension of point neuron models such as the AdEx neuron makes it impossible to directly couple an external electric field that could affect the internal membrane voltages. Following Ref. [28], we obtain an equivalent electrical input current to a point neuron by matching it to reproduce the effects of an oscillating extracellular electric field on a spatially extended ball-and-stick (BS) model neuron of a given morphology. We calculated the equivalent current amplitudes for the exponential integrate-and-fire (EIF) neuron, which is the same as the AdEx neuron without somatic adaptation (a = b = 0). In the case with adaptation, the translation from current to field works for high frequency inputs only and the approximation breaks down for slowly oscillating inputs. Thus, we have limited our estimated field strengths to the case without adaptation.

The amplitude of the equivalent input current that causes the same subthreshold depolarization of a (linearized) EIF neuron as the somatic depolarization caused by the effects of an oscillatory electric field on the BS neuron’s dendrite is then calculated using (24) A is the amplitude of the electric field in V/m, U_BS(f) is the frequency-dependent polarization transfer function of the BS neuron and z_EIF(f) is the impedance of the EIF neuron which are both given by (25) (26) with the following substitutions: (27) (28) (29) (30) The BS neuron we used to estimate electric field strengths has the following parameters: The soma has a diameter of d_s = 10 μm, a specific membrane capacitance of C_m = 10 mF/m² and a membrane resistivity of ρ_s = 2.8 Ωm². The dendritic cable has a length of l_d = 1200 μm, a diameter of d_d = 2 μm (as in the typical range of cortical pyramidal cells [32, 84]), a membrane resistivity of ρ_m = 2.8 Ωm² and an axial resistivity of ρ_a = 1.5 Ωm.

Using these parameters, an step input using an electric field with an amplitude of 1 V/m changes the somatic membrane potential by about 0.5 mV from its resting potential of −65 mV which is in agreement with in vitro measurements [32]. The curves shown in Fig 8 translate an electric field of a given amplitude and frequency to a corresponding input current and vice versa. An increase of the mean membrane current by 0.1 nA corresponds to an increase of the static electric field strength by 20 V/m.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Conversion between electric field amplitudes and equivalent input currents.

Each curve shows the frequency-dependent amplitude of an equivalent input current in pA to an exponential integrate-and-fire neuron with parameters as defined in Table 1 when the electric field amplitude acting on an equivalent ball-and-stick neuron is held constant. Electric field amplitudes in V/m are annotated for each curve.

https://doi.org/10.1371/journal.pcbi.1007822.g008

Numerical simulations

The mean-field equations were integrated using the forward Euler method. In Fig 2, each time series for a set of external inputs and in the bifurcation diagrams was obtained after t = 5 s simulation with an integration timestep of dt = 0.05 ms. In Fig 3 we simulated each point for t = 10 s with dt = 0.01 ms. For Fig 5 we simulated for t = 30 s with dt = 0.05 ms.

The spiking network model was implemented using BRIAN2 [85] (2.1.3.1) in Python. The equations were integrated using the implemented Heun’s integration method. An integration step size of 1 ms was used. In all network simulations, N_e = N_i. For the bifurcation diagrams in Fig 2, we used N = 50 × 10³ (i.e. 25 × 10³ per population), a total simulation time of t = 6 s and in Fig 3, N = 20 × 10³ and t = 6 s. The stimulation experiments in Fig 4 used N = 100 × 10³, t = 3 s. The spectra in Fig 5 used N = 20 × 10³, t = 6 s. Phase locking plots in Fig 6 used N = 20 × 10³, t = 20 s.

Benchmarking the AdEx network with N = 100 × 10³ on a single core took around 10⁴ times longer to run than the corresponding mean-field simulation. This does not include the time required for initializing the simulation, such as setting up all synapses, which can also require a comparable amount of time. The computation time scales nearly linearly with the number of neurons.

An implementation of the mean-field model is available as a Python library in our GitHub repository https://github.com/neurolib-dev/neurolib. The Python code of the comparison to AdEx network simulations, the stimulation experiments, as well as the data analysis and the ability to reproduce all presented figures in this paper can be found at https://github.com/caglarcakan/stimulus_neural_populations.