Neural mass models

Neural mass models represent one of the different formulations of rate models. They are a subtype of neural field models [Deco et al., 2008], that assume a spatial homogeneity within the neuronal population, thus neglecting spatial extensions. Below the formalism of neural mass models, that dates back to the work of [Lopes da Silva et al., 1974], [Nunez, 1974], and more recently [Liley et al., 1999] is introduced. The main variable is the average membrane voltageV_m, which is transformed into a firing rate via a sigmoidal mapping [Jansen et al., 1993, Robinson et al., 1997]

Q_k(V_k) = Q^max_k

1 + exp (−C(V_k−θ_k)/σ_k). (2.4) Here Q^max accounts for the maximal firing rate of the respective population, while θ describes the firing rate threshold of the population and σ the slope or steepness of the transition. C = π/√

3 acts as a scaling factor that relates the slope of the sigmoid function to the standard deviation of its derivative.

Spikes travel along the axonal fibers to the dendrites of the target population. For a

12 2.3. NEURAL MASS MODELS fixed axonal length, this can be described by a damped wave equation [Robinson et al., 1997]

φ¨_k=ν² Q_k(V_k)−φ_k+v²5²φ_k

−2νφ˙_k, (2.5)

that is related to the corresponding cable equation. Herevdescribes the mean axonal conduction speed andνthe product ofvand the mean inverse length of the axons [Steyn-Ross et al., 2005]. However, as neural mass models average over spatial extensions the spatial derivate can be omitted and the time evolution of the incoming spike rate φk

simplifies to

φ¨_k=ν²(Q_k(V_k)−φ_k)−2νφ˙_k, (2.6) Note that this is equivalent to the convolution of the firing rate with an alpha function α_ν representing the average axonal conduction delay

φ_k=αν(t)⊗Q_k, αν(t) =ν²texp(−νt).

For local connections within a neural populationν1holds true and the alpha function approaches theδ-distribution. In that caseφ_kcan be replaced by the instantaneous firing rateQ_k(V_k).

As second order differential equations are generally difficult to handle, they are usu-ally split into two first oder differential equations by introduction of an additional variable

φ˙_k =y_k,

˙

y_k =ν²(Q_k(V_k)−φ_k)−2νy_k.

(2.7) At the dendrites incoming spikes elicit transmitter release into the synaptic cleft.

The fraction of open synaptic channels s_lk, can be described by a convolution with an alpha functionα_l, that represents the average response to a single spike [Tuckwell, 1988].

s_lk= X

k⁰tok

α_l(t)⊗N_kk⁰φ_k⁰, α_l(t) =γ_l²texp(−γ_lt).

The sum is over all spikes from sources k⁰ that emerge at the synapses type l of populationk, scaled byN_kk⁰ that accounts for the mean number of synaptic connections originating from presynaptic population k⁰ to postsynaptic population k. The inverse rise time γ_l determines the shape of the response. Similarly to (2.6) the convolution integral can be transformed into a second order differential equation:

¨

Please note, that the alpha function is often written differently as

α_l= Γ_lγ_lexp(−γ_lt). (2.9) There, Γ_l covers the amplitude of the synaptic response, which allows for a better fit of the synaptic response to experimental measurements. However, while this suggests a more physiological description this notation has two disadvantages. First, it breaks the symmetry of Eq. (2.8) and second it can be shown, that only the full product of

ΓlγlNkk⁰Qk⁰ determines the amplitude of the overall synaptic response, resulting in an over-parametrization of the model. Therefore, the formulation of Eq. (2.6) has been chosen, absorbing the scaling of the synaptic response function into the connectivities.

The majority of neural mass models determines the EEG signal by directly summing up the synaptic inputs s_lk [Jansen et al., 1993, David and Friston, 2003, Moran et al., 2007, Ursino et al., 2010, Bhattacharya et al., 2011]. This is based on the assumption that there is an equilibrium state, the system is always close to [Wilson and Cowan, 1973, Robinson et al., 1997], resulting in a linear approximation of the neuronal response to synaptic input, that treats inputs independently from the state of the system. How-ever, for larger oscillations and events like evoked potentials or slow oscillations, this assumption is invalid.

A generalization was presented by [Liley et al., 1999, Liley et al., 2002] with the introduction of a weighting functionΨ, that scales the inputs with respect to the state the population is in. A comparison with the detailed Hodgkin-Huxley type models presented in Section 2.1 reveals that the scaling function is indeed a normalized formulation of the driving term (V_k −E_l) utilized in Eq. (2.3). Consequently, the Liley model can be rewritten with respect to physiological conventions from the conductance based models, leading to a formulation that exemplifies its similarity with the conductance based models Hodgkin-Huxley and provides a close link to physiology. The evolution of the mean membrane voltage of populationk then follows

τ_kV˙_k =−(V_k−E_L^k)−X

l

g_ls_lk(V_k−E_l). (2.10)

Here,gdepicts the weight andEthe Nernst potentials of the respective channel. Note that in contrast to the Hodgkin-Huxley model,gdoes not carry the unit of a conductance.

To better discriminate between the two we denote the physical contuctances withg. The¯ spike generating currentsIKandINafrom Eq. (2.1) have been replaced by the firing rate function Eq. (2.4), reducing the computational load significantly.

2.4 Limitations

The fundamental assumption of the neural mass models presented in Section 2.3 is the mapping between firing rate and the average membrane voltage of the neural population given in Eq. (2.4). While this assumption is valid for most cases, there are notable exceptions, with spike-frequency adaptation (SFA) and bursting modes being the most prominent ones [Fuhrmann et al., 2002, Contreras et al., 1992].

SFA occurs due to a slow negative feedback, that reduces the excitability of the neuron. Physiologically there are three distinct mechanisms that cause SFA [Benda and Herz, 2003]. The first is the partial inactivation of the spike generating sodium channels thus reducing spiking activity [Fleidervish et al., 1996]. The other mechanisms are based on the activation of slow hyperpolarizing currents either due to the increase in membrane voltage (M-type currents) [Brown and Adams, 1980, Pospischil et al., 2008] or through the change in certain ion concentrations (A-type currents) [Madison and Nicoll, 1984, Compte et al., 2003].

Next to tonic firing, bursting is one of the fundamental firing modes observed in

14 2.5. HYBRID MODELS vivo. During a burst, a neuron generates a short series of spikes, that is followed by a prolonged period of quiescence. That behavior is usually achieved by the interaction between a fast spiking and a slow adaptive current, that modulates the spiking, and is associated with synchronization of neural populations and the generation of motor patterns. A special case are thalamic relay neurons, where the fast spiking is generated through deinactivation of T-type calcium currents, that are then inactivated by the generated activity [Huguenard and McCormick, 1992, Destexhe et al., 1996b], but also have a slow anomalous rectifier current, that modulates the amplitude [Destexhe et al., 1996a].

While slow adaptive currents are intrinsic mechanisms that depend on the activity of the respective individual neuron, they greatly affect population dynamics, which is exemplified by slow oscillatory activity in the cortex or spindle activity in the thalamus.