# Modeling Neural Activity

Authors: Wim van Drongelen
Journal: ISRN Biomathematics
 (2013)

---

## Abstract

This paper provides an overview of different types of models for studying activity of nerve cells and their networks with a special emphasis on neural oscillations. One part describes the neuronal models based on the Hodgkin and Huxley formalism first described in the 1950s. It is discussed how further simplifications of this formalism enable mathematical analysis of the process of neural excitability. The focus of the paper’s second component is on network activity. Understanding network function is one of the important frontiers remaining in neuroscience. At present, experimental techniques can only provide global recordings or samples of the activity of the huge networks that form the nervous system. Models in neuroscience can therefore play a critical role by providing a framework for integration of necessarily incomplete datasets, thereby providing insight into the mechanisms of neural function. Network models can either explicitly contain individual network nodes that model the neurons, or they can be based on representations of compound population activity. The latter approach was pioneered by Wilson and Cowan in the 1970s. Finally I provide an overview and discuss how network models are employed in the study of neuronal network pathology such as epilepsy.


---

## Body

## 1. Introduction

A fundamental and famous model in neuroscience was described by Hodgkin and Huxley in 1952 [1]. They generated a formalism for the dynamics of the membrane potential of the giant axon of squid describing measurements of sodium, potassium, and leakage currents. Their model can be represented by an equivalent electronic circuit of the axon’s membrane in which a capacitor models the membrane’s phospholipids and several voltage-dependent resistors represent its ion channels. Much later, after computer technology became readily available, their formalism has been widely employed and also includes other types of channels. In addition, it is used to create detailed cell models in which the cell is divided into a set of coupled compartments each represented by a separate membrane model. In many studies these computational models are embedded in networks (e.g., [2–7]). In this approach, the individual nodes of the network and their connections are simulated with the purpose of finding the properties associated with emergent network activities such as generation of synchronized bursts and oscillations. The enormous advantage of these models is that they are close to the experimental domain; that is, they may include recorded ion conductance values, morphology of observed cell types, and known intercellular connectivity. In spite of the many parameters in these models, an overwhelming amount of detail is still missing.  As a consequence, due to nonlinearities in neuronal function (e.g., membrane conductance dynamics, the coupling function between neurons), small inaccuracies in the model parameters may lead to large prediction errors of the model. In spite of this potential shortcoming, some models perform surprisingly well in predicting or mimicking experimental outcomes (e.g., [8]).An additional problem with these detailed models is that, because of their complexity, they cannot be analyzed mathematically and may not therefore lead to deeper understanding of the process it models. Further model reduction is required before mathematical analysis can be employed. For example, the Hodgkin and Huxley model is a four dimensional nonlinear model (the dynamic variables are membrane potentialV, and gating variables m, h, and n). Both FitzHugh [9] and Nagumo et al. [10], inspired by the nonlinear van der Pol oscillator [11], reduced the Hodgkin and Huxley formalism to one pair of dynamic variables (membrane potential V and a recovery variable w) to gain insight in the dynamics of the excitable membrane. Later a similar simplified approach was published for the excitation in the barnacle giant muscle fiber by Morris and Lecar [12].A further simplification for simulation of nerve cell activity is the leaky integrate-and-fire (IF) model [13]. The basic IF model describes a neuron’s subthreshold behavior as a leaky integrator (a resistor and capacitor) and adds an ad hoc superthreshold spike when a preset threshold is exceeded. In between the spike events, the IF circuit is a one-dimensional linear model for the cell’s subthreshold activity. There are additions to this model, for instance, resonance properties are included by making the model two-dimensional and a quadratic function has been proposed to add a nonlinear component that is capable of generating the spike upward deflection [14]. More recently, Izhikevich [15] extended and combined these approaches into what he defines as the simple model of choice (SMC). In this model there is no threshold for the spike generator and only the downward deflection (reset) of the spike is generated in an ad hoc fashion.When modeling neuronal networks, one can start from the network nodes, that is, the neurons, or make direct assumptions at a higher organization level, that is, the population activity. In the first case, the bottom-up approach, one must decide how to model the nodes. Simple on/off switches are the most abstract models of neurons in networks. The underlying hypothesis of this approach is that, for some of the aggregate behavior, the network properties themselves are of principal interest while the details of the cellular properties (beyond having active and inactive states) are irrelevant [16–20]. Many network models that use the bottom-up approach include more complex representations for the nerve cells and their connections. These representations may vary from a simple IF model to compartmental models including biophysically realistic channels (e.g., [2–7, 21–27]). Recently, a network model of IF neurons was connected to an artificial, simulated eye and arm, thereby creating a system that is capable of replicating human behavior [28].A different family of network models describes aggregate behavior of neuronal populations. Although the components of these models are derived from or inspired by the neuronal function, they do not explicitly consist of a network of individual agents. This approach was pioneered by Wilson and Cowan [29, 30]. Their result is a nonlinear, two dimensional model describing the dynamics of excitatory and inhibitory activity. The expressions for the population activity are based on assumptions about cellular statistics such as the distribution of neuronal excitability. The initial 1972 Wilson-Cowan model describes only the time-dependent dynamics of the neural populations; whereas their 1973 paper includes a spatial component to this model. In the same decade, Nunez proposed models of cortical networks to explain and describe the electroencephalogram (EEG) representing compound activity of millions of nerve cells [31–34]. This work focused on the synaptic activities in the cortical pyramidal cell population since they are considered to be the major contributors to the mammalian EEG. To describe cortical activity, Nunez included inhibitory and excitatory populations and used linearized relationships to describe their interactions. Since the early 1970s, many models have been inspired by the Wilson and Cowan approach and the models by Nunez. Examples are the models by Lopes da Silva et al. [35] and van Rotterdam et al. [36], the Jansen neural mass model [37, 38], and Freeman’s models of the olfactory system [39, 40]. More recent examples are the models by Wendling and coworkers (e.g., [41]), Liley and Bojak [42], van Albada and Robinson [43], and Hindriks and van Putten [44]. A model developed by Jirsa and Haken [45] produces a field equation for the network activity generated by a cortical surface and can be considered a synthesis of the different population model approaches. In part due to increased access to computer technology facilitating the analysis of nonlinear dynamics, the original Wilson-Cowan approach and its derived models are being “rediscovered” [46]. Due to the deterministic character of these models, they fail to model fluctuations. Stochastic approaches are capable of including higher order moments that characterize the neural activity such as synaptic input and spike trains (e.g., [47, 48]). Recently, Buice and Cowan [49] used a statistical mechanics approach to include fluctuations in the population model.Understanding the function of neuronal networks is a frontier in current neuroscience. Our lack of fundamental understanding of network function in the nervous system prevents us from linking neuronal activity to higher level physiologic and pathologic behavior. Although compound activity of huge networks can be captured by the EEG or possibly by functional magnetic resonance imaging (fMRI), current experimental techniques in neuroscience cannot record the individual activity of cells forming larger networks (>1000 neurons) with sufficient resolution in the time domain (<1 ms). Due to the lack of such techniques, there is a real need for theoretical and computational models in neuroscience because they can provide a framework for the integration of experimental data.  In addition, models provide an efficient tool for generating and testing series of hypotheses; to cite Abbott [50], “A skillful theoretician can formulate, explore, and often reject models at a pace that no experimental program can match.”The purpose of this paper is to provide a summary of different models for neural activity with a focus on oscillatory patterns, to discuss the relationship between modeling approaches, and to place them in a logical and a historical context. Section2 summarizes the Hodgkin-Huxley formalism and its simplifications. The next part, Section 3, describes the network models with nodes that are based on the models for the individual neurons, followed by Section 4 about population cell models without explicit individual neurons. Section 5 discusses a procedure and effects of adding a stochastic component. The application of neuronal models in studying pathology such as epilepsy is evaluated in Section 6. The final section contains the concluding remarks.


## 2. Cellular Models

The Hodgkin and Huxley (H&H) equations [1] describe the nonlinear dynamics of the electrical properties of an excitable membrane of the giant nerve fiber that innervates the muscles involved in water jet propulsion in the squid. The inside and outside of the nerve fiber is separated by a biomembrane (Figure 1(a)). The principal component of the biomembrane is a double layer of phospholipids that serve as an insulator between the conductive media in- and outside of the fiber; this component can be modeled as a capacitor (C in Figure 1(a)). Ion pumps in the membrane create gradients of different ion species, such as sodium and potassium, across the membrane, and ion channels allow passive leak of the ions.(a) Diagram of a biomembrane that separates the inside and outside of the nerve cell. See text for further explanation. (b) Membrane equivalent circuit withk ion channels. E1,E2,…,Ek are the Nernst potentials for the individual ion species; R1,R2,…,Rk and i1,i2,…,ik represent the resistance of the different ion channels and their currents, respectively. IInject is the current injected. The membrane capacitance and the associated current are C and IC. V is the membrane potential. The equation IC+I=IInject is Kirchhoff’s law applied to the circuit.
(a)
(b)Because ions are charged, an ion pump creates a potential difference and can therefore be modeled as a potential source (E in Figure 1(a)). For individual ion species, these potentials are known as Nernst potentials that can be computed by the Nernst equation: for example, at 25°C, the Nernst potential for the positively charged sodium and potassium ions is determined from their intra- and extracellular concentrations ([Cin] and [Cout], resp.) as 58log10([Cout]/[Cin]) mV. An ion channel conducts one or more ion species and may be modeled as a resistor (R in Figure 1(a)). The resistances of the sodium and potassium channels depend on the potential across the biomembrane. This membrane potential dependence is modeled by the H&H equations using activation and inactivation processes that govern the degree of opening of the channel and consequently its conductance. It is of interest to note that Hodgkin and Huxley proposed their formalism in the 1950s as an empirical model, and that molecular evidence for the existence of the different membrane components came decades later!In the H&H formalism, the potential difference between the fiber’s in- and outside, the membrane potential (V) is related to activation and inactivation parameters of sodium and potassium ion channels (m, n,  h). The differential equations that H&H used to describe the dynamics are
(1)CdVdt=-ILeak-INa-IK+IInject
or
(2)CdVdt=-gLeak(V-ELeak)-g-Nam3h(V-ENa)-g-Kn4(V-EK)+IInject.
The parameter’s dynamics are governed by three first-order differential equations
(3)dmdt=m∞(V)-mτm(V),(4)dhdt=h∞(V)-hτh(V),(5)dndt=n∞(V)-nτn(V).
Here, the variables τm, τh, and τn  are the time constants of parameters m, h, and n. Equation (2) is based on Kirchhoff’s law stating that the sum of all membrane currents is equal to zero. An equivalent membrane circuit is depicted in Figure 1(b). Based on the equation Q=CV (with Q—charge of the capacitor, C—its capacitance, and V—potential across the capacitor), we determine that current IC equals C(dV/dt). Further we determine that each ion current ik is described by Ohm’s law gk(V-Ek), with C—membrane capacitance; V—membrane potential; gk—conductance (=1/Rk) of ion species k, which can be presented in (2) by the product of the maximum conductance g-k and activation and inactivation parameters; Ek—equilibrium potential for ion species k. Equations (3)-(4) describe the dynamics for the activation (m) and inactivation (h) parameters of the sodium channel. Equation (5) describes the dynamics of the activation (n) of potassium. In the following, we often represent activation/inactivation variables generically by x. It can be seen in (3)–(5) that the parameters  x∞ and τx depend on the membrane potential V. These V-dependent relationships are all nonlinear as well as the activation parameters m and n since H&H determined experimentally that the model fits best if they used m3 and n4 in (2). The dynamics of the parameters in (3)–(5) were originally presented by H&H in the form: dx/dt=αx(1-x)-βxx, where αx and βx are rate constants that are membrane potential (voltage) sensitive. These presentations can be linked to the one in (3)–(5) by τx=1/(αx+βx) and x∞=αx/(αx+βx). The voltage-dependent functions for αx and βx were determined by H&H using the so-called voltage-clamp technique. This technique enables the experimenter to keep V constant, the holding potential, while measuring membrane current. By repeating this measurement for a range of holding potentials, H&H determined the I-V relationship. This relationship was determined for sodium and potassium currents (a specific ion current can be determined by disabling other ion currents pharmacologically). In their model, H&H used these I-V relationships to determine αx and βx. As it can be seen in (2), the total membrane current is proportional to the time derivative of the membrane potential dV/dt=V˙; therefore, the experimentally determined I-V relationship can be interpreted as the V˙-V phase space for membrane potential dynamics (e.g., [15, 51]).


### 2.1. Linearization of the Hodgkin and Huxley Equations

A linearized version of the H&H equations is useful to examine subthreshold oscillatory properties of the biomembrane around a resting or holding potential (sayV*). The mathematical details of the following can be found in Appendix A. Here I describe the details of the linearization procedure so that the physiological parameters in the full model can be related to its linearized version; readers not interested in this aspect can skip to Section 2.1.3.In order to linearize the nonlinear equations, we rewrite (2)–(5) in a more compact form as C(dV/dt)=f(V,m,h,n,IInject) and dx/dt=f(V,x), with x representing an activation/inactivation variable m, h, or n. Then we linearize about an equilibrium potential V*. We define small variations in the parameters about this equilibrium as δV, δm, δh, δn, δI. Using this, we obtain a linearized equation (2)
(6)CdδVdt=(∂f∂V)*δV+(∂f∂m)*δm+(∂f∂h)*δh+(∂f∂n)*δn+(∂f∂IInject)*δIInject,
with: f=f(V,m,h,n,IInject), (∂f/∂V)*=-(gLeak+g-Nam3h+g-Kn4)|*, (∂f/∂m)*=-g-Na3m2h(V-ENa)|*, (∂f/∂h)*=-g-Nam3(V-ENa)|*, (∂f/∂n)*=-g-K4n3(V-EK)|*, (∂f/∂IInject)*=1.For the gating parameters ((3)–(5)), using dx/dt=(x∞(V)-x)/τx(V)=f(V,x), we get the linearized expression
(7)dδxdt=(1τx(V))*(dx∞(V)dV)*δV-(1τx(V))*δx.
In the literature, it is common to further simplify the notation of (6) and (7). For example in Richardson et al. [52], (6) and (7) are simplified to (Appendix A)
(8)Cdvdt=-gMv-∑x=1Ngxwx+δIInject︸Iapp,τx(V)dwxdt=v-wx.


#### 2.1.1. Linearization of the Potassium Channel

Now we consider the potassium current and model small perturbations and show we can model this by using a two branch electrical circuit with resistorsR, RK, and inductance L (Figure 2(a)). A small perturbation in the circuit with the two parallel branches is governed by
(9)δV=δI1/RK+1/(R+L(d/dt))︸ZRL.
Equation (9) is just Ohm’s law (V=IZ) applied to the circuit in Figure 2(a). Here Z, the inverse of 1/RK+1/(R+L(d/dt)), represents the impedance of the circuit that links the fluctuations in potential and current (R and L are in series and they are parallel to RK). Now if we consider a single ion current, the K+ current, we show that we can write this in the same form as (9). (Experimentally, a single ion current can be measured by disabling all other currents pharmacologically.)Figure 2
Diagrams of equivalent circuits representing ion channels. (a) A two branch equivalent circuit for a linearized potassium channel. The potassium channel has properties similar to an inductor. A change in current flowing through an inductor creates a time-varying magnetic field inside the coil; this induces a potential that opposes the change in current that created it. In the potassium channel, a change of K+ current is opposed by the change in membrane potential and the associated change in K+-conductance caused by the current. (b) Equivalent circuit for linearized sodium activation and inactivation. (c) A part of the linearized membrane model for the sodium channel with negative values for both Ln and rn (c1) and its equivalent RC circuit (c2).First we restate, from (2), that the potassium current IK is given by
(10)IK=g-Kn4(V-EK).
A small perturbation approximated by a linearization around potential V* is
(11)δIK=(∂IK∂V)*δV+(∂IK∂n)*δn,
with (∂IK/∂V)*=g-Kn4|*=GK*, and (∂IK/∂n)*=4g-Kn3(V-EK)|*.We can now use (7), substituting n for x, to relate δn to δV(12)δn=(1/τn)*(dn∞/dV)*(1/τn)*+d/dtδV.
Substituting this result into (11) gives
(13)δIK=GK*δV+(4g-Kn3(V-EK))*(1/τn)*(dn∞/dV)*(1/τn)*+d/dtδV.
This result can be written as
(14)δV=δIK(GK*+[A*/(B*+d/dt)])
with A*=(4g-Kn3(V-EK))*(1/τn)*(dn∞/dV)* and B*=(1/τn)*.Now we can directly relate the linearized potassium channel to the circuit in Figure2(a). Comparing (9) and (14), we can determine that both expressions are equivalent if RK=1/GK*, R=B*/A*, and L=1/A*.Combining the above we find that(15)L=RB*=Rτn*,R=B*A*=1(4g-Kn3(V-EK))*(dn∞/dV)*.
Thus we can successfully model small perturbations in the potassium channel with the equivalent circuit in Figure 2(a). Note that A* is positive because dn∞/dV and V-EK are both positive at physiological values of V*. Therefore L is positive and (because B* is also >0) R is positive.


#### 2.1.2. Linearization of the Sodium Channel

Now we apply the linearization to inward currentINa(16)INa=g-Nam3h(V-ENa).
A small perturbation approximated by a linearization around potential V* is
(17)δINa=(∂INa∂V)*δV+(∂INa∂m)*δm+(∂INa∂h)*δh.
We simplify notation
(18)δINa=(GNa*+[A*B*+d/dt]+[D*E*+d/dt])δV
using: (∂INa/∂V)*=g-Nam3h|*=GNa*, A*=(3g-Nam2h(V-ENa))*(1/τm)*(dm∞/dV)*, B*=(1/τm)*, D*=(g-Nam3(V-ENa))*(1/τh)*(dh∞/dV)*, and E*=(1/τh)*.From the three terms in (17), it is obvious that we need three branches in an equivalent electrical circuit for the sodium channel. One resistor RNa and two branches each with a resistor and inductor, one for activation m and one for inactivation h (Figure 2(b)). From the circuit properties in Figure 2(b), we determine that a perturbation of the potential is governed by
(19)δV=δI1/RNa+1/(Rm+Lm(d/dt))︸ZRLm+1/(Rh+Lh(d/dt))︸ZRLh.
We find that the expressions in (18) and (19) are equal if RNa=1/GNa*, Rm=B*/A*, and Lm=1/A* with Lm=Rm/B*=Rmτm*, Rm=B*/A*=1/(3g-Km2h(V-ENa))*(dm∞/dV)*, Rh=E*/D*, and Lh=1/D*.Combining the above we findLh=Rh/E*=Rhτh* and Rh=E*/D*=1/(g-Km3(V-ENa))*(dh∞/dV)*.So we can successfully model small perturbations in the sodium channel with the equivalent circuit in Figure2(b). The model can easily be evaluated by simulation; however, building the model using hardware components cannot be done as simply because some of the conductance and inductance values are negative over the range we are interested in. Note that in the expressions above A* is negative because dm∞/dV>0 and V-ENa is negative at physiological values of V*. Therefore Lm is negative and (because B* is >0) Rm is also negative. In contrast, note that D* is positive because both dh∞/dV and V-ENa<0 at V*. Therefore Lh is positive and (because E* is >0) Rh is also positive. Mauro et al. [53] and Koch [54] proposed to replace the negative inductance by a positive capacitance and correct for negative resistance in one of the other branches of the circuit. They show that a branch with negative conductance and inductance parallel to a true (positive valued) resistor can be replaced with an RC circuit without changing the effective impedance of the two branches. The equivalent diagrams of the circuitry are shown in Figure 2(c). In the original circuit (Figure 2(c1)) with the positive resistor RK, negative resistance rn, and negative inductor Ln, we have impedance Zo in the original circuit
(20)Zo=(1RK+1rn+jωLn)-1=((RK+rn)+jωLnRKrn+jωRKLn)-1.
In the replacement circuit (Figure 2(c2)) with the capacitor we have impedance Zr(21)Zr=(RK+rnRKrn+1(-rn+1/jωCn))-1=(-RKrn+RK/jωCn-rn2+rn/jωCn+RKrn(-RKrn2+RKrn/jωCn))-1.
After substituting Cn=-Ln/rn2 and multiplying numerator and denominator by -jωLn/rn2 we get
(22)Zr=((RK+rn)+jωLnRKrn+jωRKLn)-1.
This result is identical to the original impedance Zo; it will therefore also have the same frequency response: absolute value and phase (i.e., the same Bode plot).


#### 2.1.3. Generalized Equivalent Circuit

In the previous sections, it was shown that the channels in a linearized version of the H&H equations can be represented by equivalent circuits consisting of a network of resistors and inductors (Figure2). Even in the case where negative values occur for these components (due to positive feedback between sodium conductance and membrane depolarization), they can be represented by positive valued alternative electronic components (Figure 2(c)). When combining these findings, we simplify the equivalent circuit in Figure 1 and consider the RCL network in Figure 3 as a linearized representation of the membrane.Figure 3
Electrical equivalent circuit of the linearized H&H equations.The inductor channel in the circuit represents the stabilizing effects of the outward current activation (e.g., potassium) and inward current inactivation (e.g., sodium). The resistor and capacitor represent membrane resistance and capacitance plus corrections for negative impedances. Since we look at small perturbations around an equilibrium, which we conveniently set at zero, the driving forces of the difference between Nernst potential and resting potential can be ignored.The pair of differential equations that govern the circuit in Figure3 are
(23)CdVdt+VR+IL=IInject+Synapse,LdILdt=V-RLIL.
We now simplify and we nondimensionalize the system. First we substitute w=RLIL and I0=IInject+Synapse and get C(dV/dt)+V/R+w/RL=I0 and (L/RL)(dw/dt)=V-w. Then we change to a timescale τ such that τ=(RL/L)t→dt=(L/RL)dτ; substituting this we get the following pair of ODEs:
(24)CdV(L/RL)dτ+VR+wRL=I0→dVdτ+(LCRL)VR+(LCRL)wRL=(LCRL)I0,LRLdw(L/RL)dτ=V-w→dwdτ=V-w.
Finally, using V=w′+w and V′=w′′+w′ (note that w′=dw/dτ), we obtain the following 2nd order ODE:
(25)w′′+w′+(LCRRL)w′+(LCRRL)w+(LCRL2)w=(LCRL)I0→w′′+[1+(LCRRL)︸α]w′+[(LCRRL)+(LCRL2)︸α+β]w=(LCRL)I0.
This result for α and β is identical to the equations analyzed by Erchova et al. [55] (their equations (A16) and (A17)). An alternative route to create a single equation from the two original ones is by direct substitution while leaving dimensionality intact. If we differentiate the first expression in (23) with respect to time we have CV¨+V˙/R+I˙L=I˙0. Here V¨, V˙, and I˙ indicate the second and first time derivatives of V and I. We can use the second expression in (23) to substitute I˙L=(V-RLIL)/L where IL as a function of V can be found from CV˙+V/R+IL=I0. This yield
(26)CV¨+V˙R+V-RL(I0-CV˙-V/R)L=I˙0→CV¨+(1R+RLCL)︸γV˙+1L(1+RLR)︸εV=RLLI0+I˙0.
This result is identical to the expression used by Erchova et al. [55] (their Equation (A1)). Interestingly this equation is also similar to the network equation by Jirsa and Haken [45] in which the spatial derivative is set to zero (see Section 4.5). This similarity across the scale from neuron to network is seen frequently. In one example, investigators even use a cellular model for action potential generation to represent activity of a population [56]. Although Curto et al. [56] do not explicitly describe this, in this case the inward and outward currents of the cell membrane can be considered to represent the network’s excitation and inhibition, respectively.


#### 2.1.4. The Frequency Response of the Linearized Model

In electrophysiology, cellular properties are often probed with hyperpolarizing and depolarizing current steps. The unit impulse response, the derivative of the unit step, is the inverse Fourier transform or inverse Laplace transform of the system’s frequency response or transfer function, respectively. Since we deal with a linearized system with currentI as input and membrane potential V as output, the frequency response is just the impedance Z (i.e., Z(ω)=V(ω)/I(ω), with ω—frequency in rad/s). We can use the Laplace transform of (26) to obtain the transfer function H(s) of the linearized model: H(s)=V(s)/I(s)=(ε+γs+Cs2)/((RL/L)+s). In this form we can apply the Matlab freqs command to depict the Bode plot of the model (Figure 4).Bode plot of the linearized model depicted in Figure3. The parameters used to prepare the plots are C=100 pF; R=200 MΩ; RL=20 MΩ; L=2 MH. Note in the upper panel that a resonance peak is present at 75 rad/s.
(a)
(b)An interesting observation is that the Bode plot in Figure4 shows a peak in the magnitude of the impedance, indicating that neurons and their networks display resonance [55, 57–60]. With an appropriate, realistic choice of parameters, this peak is located  in the alpha-frequency range of the EEG at 75 rad/s (~12 Hz).


### 2.2. Models That Can Be Derived from the Hodgkin and Huxley Formalism

The H&H representation and its linearized version can be seen as a basis for many different neuronal models. Of course one could either add complexity to the H&H formalism or further reduce it. The following sections summarize representative examples of both types.


#### 2.2.1. Complex Cell Models

Using the H&H formalism one might extend the sodium and potassium channels of the model neurons to also include different types of potassium channels (such as the A-current, a transient potassium current), persistent sodium channels, and different types of Ca++ channels. In addition to the voltage sensitive channels, it is common practice to also include Ca++ and/or Na+ sensitive channels. In addition to the addition of channels, the cell model can also be divided into many multiple coupled compartments (e.g., [4]). Each compartment represents a small part of the cell that can be approximated by a compartment with a uniform membrane potential. The rationale for the use of multiple equipotential compartments is to mimic the spatial separation of different functions: synaptic input, passive and active propagation of signals, and generation of action potentials. The more complex cell models are also frequently used to create networks in order to determine emergent behavior of neuronal populations (e.g., [3, 5, 6, 26, 61]). Special neural simulation packages such as Genesis and Neuron have been developed to simulate such models, even on parallel machines [62, 63]. A recent example in which many details are incorporated in the individual cell models is the Blue Brain Project [7]. The advantage of these models is that they are based on experimental data enabling the modeler to link the model’s behavior to measurements. This close relationship can also be a disadvantage of these complex models; due to the complexity, the simulation results can be as difficult to interpret as real recorded data.


#### 2.2.2. Simplified Cell Models

One can simplify the H&H equations while preserving the nonlinearity, or one might use the linearization results described in Section2.1. In fact both these approaches have been employed. In one case, the four dimensional H&H model can be reduced to a lower dimensional one and in the other case, the equivalent circuit of a linearized membrane model (Figure 3) is coupled to an ad hoc spike generator (Figure 5). The advantage of preserving the nonlinearity is that the spike generation mechanisms remain (at least in part) intact, but mathematical analysis of such a model can be fairly complicated. The linearized cell models are more accessible for analysis but require an ad hoc spike generator (Table 1).Table 1
Overview of the dimensionality of cellular models.
Amplifying
Resonant
Membrane potential(V)
Inward current activation
Inward current inactivation
Outward current activation
Leak current capacitance(not time dependent)
Physiological AP mechanism
Hodgkin and Huxley  [1]
1
(m) 2
(h) 3
(n) 4
Yes
Yes
FitzHugh-Nagumo [9, 10]
1
Instantaneous(not time dependent)
Generic recovery
(w) 2
Yes
Yes
Morris and Lecar [12] (Ca instead of Na)
1
Instantaneous(not time dependent)
No
(w) 2
Yes
Yes
IF [13]
1
No
No
No
Yes
Generator
RF [15]
1
No
No
Simplified (L) 2
Yes
Generator
QIF [14]
1
Generic amplification(quadratic)
No
No
Yes
Reset
Izhikevich SMC [15]
1
Generic amplification(quadratic)
No
Simplified (L) 2
Yes
Reset
AP: action potential; L: inductor; IF: integrate and fire; QIF: quadratic integrate and fire; RF: resonate and Fire; SMC: simple model of choice;V, m, h, n: parameters of the Hodgkin and Huxley equations; w: recovery variable. Numbers 1–4 count the dimensions of the model.Figure 5
Example of an electronic circuit of an IF neuron model using analog components. The membrane is modeled by the R and C components (orange), the threshold is implemented by the OP-Amp as a comparator (the threshold is the value at the potentiometer, orange). The reset function is performed by the Reed relay. The diodes and 10 kΩ resistors at the in-and output are included to rectify the signals. This model is linear for the subthreshold part of the activity. (b) shows a measurement of the subthreshold activity of the membrane including the reset at the threshold; (c) depicts a measurement of the spike output.(1)FitzHugh-Nagumo/Morris and Lecar.The simplification used independently by FitzHugh [9], Nagumo et al., [10], and Morris and Lecar [12], is all based on the observation that the cell has a relatively fast activating inward current and a slower activating outward current. This allows one to simplify (2)–(5) from a four dimensional model into a model with only two dynamic variables. In the FitzHugh-Nagumo approach, the inward current is sodium whereas in the Morris and Lecar model for a barnacle muscle cell, the inward current is calcium. For the following discussion, it is not critical which ion species carries the inward current, so I will focus on simplification of the H&H formalism with a Na+ current. In the original equations (2)–(5) we have V, m, h, and n. The simplified version ignores activation of the inward current m because it is so fast it is considered instantaneous. The effects of h and n are similar; that is, they repolarize the membrane, so they can either be lumped or one of the two can be ignored so that only one repolarization parameter (often indicated by w) remains. This generates a nonlinear dynamical system with two variables. Although a few variations exist, the original equations in FitzHugh [9] were
(27)x˙=c(y+x-x33+z),y˙=-(x-a+by)c.
Here x, y, and z are membrane potential, recovery variable, and injected current, respectively; a, b, and c are constants.(2)Reduced Spiking Models.From the linearization procedure we can appreciate that the neural membrane can be represented by the RCL-circuit (Figure 3). Of course, this representation only works for studying subthreshold behavior such as integration of inputs or resonance properties. This type of study explains the propensity of some neurons to resonate with specific input frequencies. If, as an additional simplification, the ion channel properties are completely ignored, we may even further simplify the membrane to an RC-circuit; in this configuration there are no resonant properties. Both these simplifications, the RCL- and RC-circuits, ignore the nonlinear properties of the excitable membrane. In the so-called integrate-and-fire (IF) models the nonlinearity is again introduced, but in an ad hoc fashion by employing a firing threshold. When this threshold is exceeded, a spike + reset are pasted into the signal. Although such an ad-hoc firing process can be included in both RCL- and RC-circuits, it is commonly applied to the latter.Historically, the RC-based integrate-and-fire model was introduced long before the H&H formalism in 1907 by Lapicque [13]. Because of the leakage current in the resistor, it is also often called the leaky integrate-and-fire neuron (LIF). Much later a nonlinear variant of the LIF was introduced where the I-V relationship is not linear as in (23), but quadratic, the QIF model (e.g., [14]). Because this quadratic relationship, the I-V curve contains a stable and an unstable equilibrium, the spike’s upstroke can be generated by the model and only an ad-hoc spike reset is required (Figure 6).Figure 6
Phase space representation for the QIF neuron. The relationship of membrane current, proportional to the derivative ofV, is plotted against membrane potential, V. Due to the quadratic relationship, there are two fixed points: one stable (the resting potential) and one unstable one, the threshold, at the border between the subthreshold (blue) and superthreshold (red) regions. The model generates the upstroke of the action potential when V exceeds the threshold. To avoid instability, an ad hoc reset is provided when the arbitrary reset value Vmax is reached.Finally, one can combine the QIF scenario with the RCL-circuit to produce a model that is capable of both resonance and spike upstroke generation. This is defined as the simple model of choice (SMC) by Izhikevich and was shown to be able to produce a plethora of neuronal behavior (e.g., [15]).


### 2.3. Summary of the Cell Models

We can consider the different models and summarize their relationships by counting the dimensions (Table1). The H&H model is four dimensional whereas the subthreshold component of the IF model is only one dimensional. Of course, dimensions are added when the ad hoc spike generator is included but most studies using this model will focus on the subthreshold behavior. Table 1 summarizes the model dimensions and provides an overview of the spiking mechanisms of the models we evaluated above.


## 2.1. Linearization of the Hodgkin and Huxley Equations

A linearized version of the H&H equations is useful to examine subthreshold oscillatory properties of the biomembrane around a resting or holding potential (sayV*). The mathematical details of the following can be found in Appendix A. Here I describe the details of the linearization procedure so that the physiological parameters in the full model can be related to its linearized version; readers not interested in this aspect can skip to Section 2.1.3.In order to linearize the nonlinear equations, we rewrite (2)–(5) in a more compact form as C(dV/dt)=f(V,m,h,n,IInject) and dx/dt=f(V,x), with x representing an activation/inactivation variable m, h, or n. Then we linearize about an equilibrium potential V*. We define small variations in the parameters about this equilibrium as δV, δm, δh, δn, δI. Using this, we obtain a linearized equation (2)
(6)CdδVdt=(∂f∂V)*δV+(∂f∂m)*δm+(∂f∂h)*δh+(∂f∂n)*δn+(∂f∂IInject)*δIInject,
with: f=f(V,m,h,n,IInject), (∂f/∂V)*=-(gLeak+g-Nam3h+g-Kn4)|*, (∂f/∂m)*=-g-Na3m2h(V-ENa)|*, (∂f/∂h)*=-g-Nam3(V-ENa)|*, (∂f/∂n)*=-g-K4n3(V-EK)|*, (∂f/∂IInject)*=1.For the gating parameters ((3)–(5)), using dx/dt=(x∞(V)-x)/τx(V)=f(V,x), we get the linearized expression
(7)dδxdt=(1τx(V))*(dx∞(V)dV)*δV-(1τx(V))*δx.
In the literature, it is common to further simplify the notation of (6) and (7). For example in Richardson et al. [52], (6) and (7) are simplified to (Appendix A)
(8)Cdvdt=-gMv-∑x=1Ngxwx+δIInject︸Iapp,τx(V)dwxdt=v-wx.


### 2.1.1. Linearization of the Potassium Channel

Now we consider the potassium current and model small perturbations and show we can model this by using a two branch electrical circuit with resistorsR, RK, and inductance L (Figure 2(a)). A small perturbation in the circuit with the two parallel branches is governed by
(9)δV=δI1/RK+1/(R+L(d/dt))︸ZRL.
Equation (9) is just Ohm’s law (V=IZ) applied to the circuit in Figure 2(a). Here Z, the inverse of 1/RK+1/(R+L(d/dt)), represents the impedance of the circuit that links the fluctuations in potential and current (R and L are in series and they are parallel to RK). Now if we consider a single ion current, the K+ current, we show that we can write this in the same form as (9). (Experimentally, a single ion current can be measured by disabling all other currents pharmacologically.)Figure 2
Diagrams of equivalent circuits representing ion channels. (a) A two branch equivalent circuit for a linearized potassium channel. The potassium channel has properties similar to an inductor. A change in current flowing through an inductor creates a time-varying magnetic field inside the coil; this induces a potential that opposes the change in current that created it. In the potassium channel, a change of K+ current is opposed by the change in membrane potential and the associated change in K+-conductance caused by the current. (b) Equivalent circuit for linearized sodium activation and inactivation. (c) A part of the linearized membrane model for the sodium channel with negative values for both Ln and rn (c1) and its equivalent RC circuit (c2).First we restate, from (2), that the potassium current IK is given by
(10)IK=g-Kn4(V-EK).
A small perturbation approximated by a linearization around potential V* is
(11)δIK=(∂IK∂V)*δV+(∂IK∂n)*δn,
with (∂IK/∂V)*=g-Kn4|*=GK*, and (∂IK/∂n)*=4g-Kn3(V-EK)|*.We can now use (7), substituting n for x, to relate δn to δV(12)δn=(1/τn)*(dn∞/dV)*(1/τn)*+d/dtδV.
Substituting this result into (11) gives
(13)δIK=GK*δV+(4g-Kn3(V-EK))*(1/τn)*(dn∞/dV)*(1/τn)*+d/dtδV.
This result can be written as
(14)δV=δIK(GK*+[A*/(B*+d/dt)])
with A*=(4g-Kn3(V-EK))*(1/τn)*(dn∞/dV)* and B*=(1/τn)*.Now we can directly relate the linearized potassium channel to the circuit in Figure2(a). Comparing (9) and (14), we can determine that both expressions are equivalent if RK=1/GK*, R=B*/A*, and L=1/A*.Combining the above we find that(15)L=RB*=Rτn*,R=B*A*=1(4g-Kn3(V-EK))*(dn∞/dV)*.
Thus we can successfully model small perturbations in the potassium channel with the equivalent circuit in Figure 2(a). Note that A* is positive because dn∞/dV and V-EK are both positive at physiological values of V*. Therefore L is positive and (because B* is also >0) R is positive.


### 2.1.2. Linearization of the Sodium Channel

Now we apply the linearization to inward currentINa(16)INa=g-Nam3h(V-ENa).
A small perturbation approximated by a linearization around potential V* is
(17)δINa=(∂INa∂V)*δV+(∂INa∂m)*δm+(∂INa∂h)*δh.
We simplify notation
(18)δINa=(GNa*+[A*B*+d/dt]+[D*E*+d/dt])δV
using: (∂INa/∂V)*=g-Nam3h|*=GNa*, A*=(3g-Nam2h(V-ENa))*(1/τm)*(dm∞/dV)*, B*=(1/τm)*, D*=(g-Nam3(V-ENa))*(1/τh)*(dh∞/dV)*, and E*=(1/τh)*.From the three terms in (17), it is obvious that we need three branches in an equivalent electrical circuit for the sodium channel. One resistor RNa and two branches each with a resistor and inductor, one for activation m and one for inactivation h (Figure 2(b)). From the circuit properties in Figure 2(b), we determine that a perturbation of the potential is governed by
(19)δV=δI1/RNa+1/(Rm+Lm(d/dt))︸ZRLm+1/(Rh+Lh(d/dt))︸ZRLh.
We find that the expressions in (18) and (19) are equal if RNa=1/GNa*, Rm=B*/A*, and Lm=1/A* with Lm=Rm/B*=Rmτm*, Rm=B*/A*=1/(3g-Km2h(V-ENa))*(dm∞/dV)*, Rh=E*/D*, and Lh=1/D*.Combining the above we findLh=Rh/E*=Rhτh* and Rh=E*/D*=1/(g-Km3(V-ENa))*(dh∞/dV)*.So we can successfully model small perturbations in the sodium channel with the equivalent circuit in Figure2(b). The model can easily be evaluated by simulation; however, building the model using hardware components cannot be done as simply because some of the conductance and inductance values are negative over the range we are interested in. Note that in the expressions above A* is negative because dm∞/dV>0 and V-ENa is negative at physiological values of V*. Therefore Lm is negative and (because B* is >0) Rm is also negative. In contrast, note that D* is positive because both dh∞/dV and V-ENa<0 at V*. Therefore Lh is positive and (because E* is >0) Rh is also positive. Mauro et al. [53] and Koch [54] proposed to replace the negative inductance by a positive capacitance and correct for negative resistance in one of the other branches of the circuit. They show that a branch with negative conductance and inductance parallel to a true (positive valued) resistor can be replaced with an RC circuit without changing the effective impedance of the two branches. The equivalent diagrams of the circuitry are shown in Figure 2(c). In the original circuit (Figure 2(c1)) with the positive resistor RK, negative resistance rn, and negative inductor Ln, we have impedance Zo in the original circuit
(20)Zo=(1RK+1rn+jωLn)-1=((RK+rn)+jωLnRKrn+jωRKLn)-1.
In the replacement circuit (Figure 2(c2)) with the capacitor we have impedance Zr(21)Zr=(RK+rnRKrn+1(-rn+1/jωCn))-1=(-RKrn+RK/jωCn-rn2+rn/jωCn+RKrn(-RKrn2+RKrn/jωCn))-1.
After substituting Cn=-Ln/rn2 and multiplying numerator and denominator by -jωLn/rn2 we get
(22)Zr=((RK+rn)+jωLnRKrn+jωRKLn)-1.
This result is identical to the original impedance Zo; it will therefore also have the same frequency response: absolute value and phase (i.e., the same Bode plot).


### 2.1.3. Generalized Equivalent Circuit

In the previous sections, it was shown that the channels in a linearized version of the H&H equations can be represented by equivalent circuits consisting of a network of resistors and inductors (Figure2). Even in the case where negative values occur for these components (due to positive feedback between sodium conductance and membrane depolarization), they can be represented by positive valued alternative electronic components (Figure 2(c)). When combining these findings, we simplify the equivalent circuit in Figure 1 and consider the RCL network in Figure 3 as a linearized representation of the membrane.Figure 3
Electrical equivalent circuit of the linearized H&H equations.The inductor channel in the circuit represents the stabilizing effects of the outward current activation (e.g., potassium) and inward current inactivation (e.g., sodium). The resistor and capacitor represent membrane resistance and capacitance plus corrections for negative impedances. Since we look at small perturbations around an equilibrium, which we conveniently set at zero, the driving forces of the difference between Nernst potential and resting potential can be ignored.The pair of differential equations that govern the circuit in Figure3 are
(23)CdVdt+VR+IL=IInject+Synapse,LdILdt=V-RLIL.
We now simplify and we nondimensionalize the system. First we substitute w=RLIL and I0=IInject+Synapse and get C(dV/dt)+V/R+w/RL=I0 and (L/RL)(dw/dt)=V-w. Then we change to a timescale τ such that τ=(RL/L)t→dt=(L/RL)dτ; substituting this we get the following pair of ODEs:
(24)CdV(L/RL)dτ+VR+wRL=I0→dVdτ+(LCRL)VR+(LCRL)wRL=(LCRL)I0,LRLdw(L/RL)dτ=V-w→dwdτ=V-w.
Finally, using V=w′+w and V′=w′′+w′ (note that w′=dw/dτ), we obtain the following 2nd order ODE:
(25)w′′+w′+(LCRRL)w′+(LCRRL)w+(LCRL2)w=(LCRL)I0→w′′+[1+(LCRRL)︸α]w′+[(LCRRL)+(LCRL2)︸α+β]w=(LCRL)I0.
This result for α and β is identical to the equations analyzed by Erchova et al. [55] (their equations (A16) and (A17)). An alternative route to create a single equation from the two original ones is by direct substitution while leaving dimensionality intact. If we differentiate the first expression in (23) with respect to time we have CV¨+V˙/R+I˙L=I˙0. Here V¨, V˙, and I˙ indicate the second and first time derivatives of V and I. We can use the second expression in (23) to substitute I˙L=(V-RLIL)/L where IL as a function of V can be found from CV˙+V/R+IL=I0. This yield
(26)CV¨+V˙R+V-RL(I0-CV˙-V/R)L=I˙0→CV¨+(1R+RLCL)︸γV˙+1L(1+RLR)︸εV=RLLI0+I˙0.
This result is identical to the expression used by Erchova et al. [55] (their Equation (A1)). Interestingly this equation is also similar to the network equation by Jirsa and Haken [45] in which the spatial derivative is set to zero (see Section 4.5). This similarity across the scale from neuron to network is seen frequently. In one example, investigators even use a cellular model for action potential generation to represent activity of a population [56]. Although Curto et al. [56] do not explicitly describe this, in this case the inward and outward currents of the cell membrane can be considered to represent the network’s excitation and inhibition, respectively.


### 2.1.4. The Frequency Response of the Linearized Model

In electrophysiology, cellular properties are often probed with hyperpolarizing and depolarizing current steps. The unit impulse response, the derivative of the unit step, is the inverse Fourier transform or inverse Laplace transform of the system’s frequency response or transfer function, respectively. Since we deal with a linearized system with currentI as input and membrane potential V as output, the frequency response is just the impedance Z (i.e., Z(ω)=V(ω)/I(ω), with ω—frequency in rad/s). We can use the Laplace transform of (26) to obtain the transfer function H(s) of the linearized model: H(s)=V(s)/I(s)=(ε+γs+Cs2)/((RL/L)+s). In this form we can apply the Matlab freqs command to depict the Bode plot of the model (Figure 4).Bode plot of the linearized model depicted in Figure3. The parameters used to prepare the plots are C=100 pF; R=200 MΩ; RL=20 MΩ; L=2 MH. Note in the upper panel that a resonance peak is present at 75 rad/s.
(a)
(b)An interesting observation is that the Bode plot in Figure4 shows a peak in the magnitude of the impedance, indicating that neurons and their networks display resonance [55, 57–60]. With an appropriate, realistic choice of parameters, this peak is located  in the alpha-frequency range of the EEG at 75 rad/s (~12 Hz).


## 2.1.1. Linearization of the Potassium Channel

Now we consider the potassium current and model small perturbations and show we can model this by using a two branch electrical circuit with resistorsR, RK, and inductance L (Figure 2(a)). A small perturbation in the circuit with the two parallel branches is governed by
(9)δV=δI1/RK+1/(R+L(d/dt))︸ZRL.
Equation (9) is just Ohm’s law (V=IZ) applied to the circuit in Figure 2(a). Here Z, the inverse of 1/RK+1/(R+L(d/dt)), represents the impedance of the circuit that links the fluctuations in potential and current (R and L are in series and they are parallel to RK). Now if we consider a single ion current, the K+ current, we show that we can write this in the same form as (9). (Experimentally, a single ion current can be measured by disabling all other currents pharmacologically.)Figure 2
Diagrams of equivalent circuits representing ion channels. (a) A two branch equivalent circuit for a linearized potassium channel. The potassium channel has properties similar to an inductor. A change in current flowing through an inductor creates a time-varying magnetic field inside the coil; this induces a potential that opposes the change in current that created it. In the potassium channel, a change of K+ current is opposed by the change in membrane potential and the associated change in K+-conductance caused by the current. (b) Equivalent circuit for linearized sodium activation and inactivation. (c) A part of the linearized membrane model for the sodium channel with negative values for both Ln and rn (c1) and its equivalent RC circuit (c2).First we restate, from (2), that the potassium current IK is given by
(10)IK=g-Kn4(V-EK).
A small perturbation approximated by a linearization around potential V* is
(11)δIK=(∂IK∂V)*δV+(∂IK∂n)*δn,
with (∂IK/∂V)*=g-Kn4|*=GK*, and (∂IK/∂n)*=4g-Kn3(V-EK)|*.We can now use (7), substituting n for x, to relate δn to δV(12)δn=(1/τn)*(dn∞/dV)*(1/τn)*+d/dtδV.
Substituting this result into (11) gives
(13)δIK=GK*δV+(4g-Kn3(V-EK))*(1/τn)*(dn∞/dV)*(1/τn)*+d/dtδV.
This result can be written as
(14)δV=δIK(GK*+[A*/(B*+d/dt)])
with A*=(4g-Kn3(V-EK))*(1/τn)*(dn∞/dV)* and B*=(1/τn)*.Now we can directly relate the linearized potassium channel to the circuit in Figure2(a). Comparing (9) and (14), we can determine that both expressions are equivalent if RK=1/GK*, R=B*/A*, and L=1/A*.Combining the above we find that(15)L=RB*=Rτn*,R=B*A*=1(4g-Kn3(V-EK))*(dn∞/dV)*.
Thus we can successfully model small perturbations in the potassium channel with the equivalent circuit in Figure 2(a). Note that A* is positive because dn∞/dV and V-EK are both positive at physiological values of V*. Therefore L is positive and (because B* is also >0) R is positive.


## 2.1.2. Linearization of the Sodium Channel

Now we apply the linearization to inward currentINa(16)INa=g-Nam3h(V-ENa).
A small perturbation approximated by a linearization around potential V* is
(17)δINa=(∂INa∂V)*δV+(∂INa∂m)*δm+(∂INa∂h)*δh.
We simplify notation
(18)δINa=(GNa*+[A*B*+d/dt]+[D*E*+d/dt])δV
using: (∂INa/∂V)*=g-Nam3h|*=GNa*, A*=(3g-Nam2h(V-ENa))*(1/τm)*(dm∞/dV)*, B*=(1/τm)*, D*=(g-Nam3(V-ENa))*(1/τh)*(dh∞/dV)*, and E*=(1/τh)*.From the three terms in (17), it is obvious that we need three branches in an equivalent electrical circuit for the sodium channel. One resistor RNa and two branches each with a resistor and inductor, one for activation m and one for inactivation h (Figure 2(b)). From the circuit properties in Figure 2(b), we determine that a perturbation of the potential is governed by
(19)δV=δI1/RNa+1/(Rm+Lm(d/dt))︸ZRLm+1/(Rh+Lh(d/dt))︸ZRLh.
We find that the expressions in (18) and (19) are equal if RNa=1/GNa*, Rm=B*/A*, and Lm=1/A* with Lm=Rm/B*=Rmτm*, Rm=B*/A*=1/(3g-Km2h(V-ENa))*(dm∞/dV)*, Rh=E*/D*, and Lh=1/D*.Combining the above we findLh=Rh/E*=Rhτh* and Rh=E*/D*=1/(g-Km3(V-ENa))*(dh∞/dV)*.So we can successfully model small perturbations in the sodium channel with the equivalent circuit in Figure2(b). The model can easily be evaluated by simulation; however, building the model using hardware components cannot be done as simply because some of the conductance and inductance values are negative over the range we are interested in. Note that in the expressions above A* is negative because dm∞/dV>0 and V-ENa is negative at physiological values of V*. Therefore Lm is negative and (because B* is >0) Rm is also negative. In contrast, note that D* is positive because both dh∞/dV and V-ENa<0 at V*. Therefore Lh is positive and (because E* is >0) Rh is also positive. Mauro et al. [53] and Koch [54] proposed to replace the negative inductance by a positive capacitance and correct for negative resistance in one of the other branches of the circuit. They show that a branch with negative conductance and inductance parallel to a true (positive valued) resistor can be replaced with an RC circuit without changing the effective impedance of the two branches. The equivalent diagrams of the circuitry are shown in Figure 2(c). In the original circuit (Figure 2(c1)) with the positive resistor RK, negative resistance rn, and negative inductor Ln, we have impedance Zo in the original circuit
(20)Zo=(1RK+1rn+jωLn)-1=((RK+rn)+jωLnRKrn+jωRKLn)-1.
In the replacement circuit (Figure 2(c2)) with the capacitor we have impedance Zr(21)Zr=(RK+rnRKrn+1(-rn+1/jωCn))-1=(-RKrn+RK/jωCn-rn2+rn/jωCn+RKrn(-RKrn2+RKrn/jωCn))-1.
After substituting Cn=-Ln/rn2 and multiplying numerator and denominator by -jωLn/rn2 we get
(22)Zr=((RK+rn)+jωLnRKrn+jωRKLn)-1.
This result is identical to the original impedance Zo; it will therefore also have the same frequency response: absolute value and phase (i.e., the same Bode plot).


## 2.1.3. Generalized Equivalent Circuit

In the previous sections, it was shown that the channels in a linearized version of the H&H equations can be represented by equivalent circuits consisting of a network of resistors and inductors (Figure2). Even in the case where negative values occur for these components (due to positive feedback between sodium conductance and membrane depolarization), they can be represented by positive valued alternative electronic components (Figure 2(c)). When combining these findings, we simplify the equivalent circuit in Figure 1 and consider the RCL network in Figure 3 as a linearized representation of the membrane.Figure 3
Electrical equivalent circuit of the linearized H&H equations.The inductor channel in the circuit represents the stabilizing effects of the outward current activation (e.g., potassium) and inward current inactivation (e.g., sodium). The resistor and capacitor represent membrane resistance and capacitance plus corrections for negative impedances. Since we look at small perturbations around an equilibrium, which we conveniently set at zero, the driving forces of the difference between Nernst potential and resting potential can be ignored.The pair of differential equations that govern the circuit in Figure3 are
(23)CdVdt+VR+IL=IInject+Synapse,LdILdt=V-RLIL.
We now simplify and we nondimensionalize the system. First we substitute w=RLIL and I0=IInject+Synapse and get C(dV/dt)+V/R+w/RL=I0 and (L/RL)(dw/dt)=V-w. Then we change to a timescale τ such that τ=(RL/L)t→dt=(L/RL)dτ; substituting this we get the following pair of ODEs:
(24)CdV(L/RL)dτ+VR+wRL=I0→dVdτ+(LCRL)VR+(LCRL)wRL=(LCRL)I0,LRLdw(L/RL)dτ=V-w→dwdτ=V-w.
Finally, using V=w′+w and V′=w′′+w′ (note that w′=dw/dτ), we obtain the following 2nd order ODE:
(25)w′′+w′+(LCRRL)w′+(LCRRL)w+(LCRL2)w=(LCRL)I0→w′′+[1+(LCRRL)︸α]w′+[(LCRRL)+(LCRL2)︸α+β]w=(LCRL)I0.
This result for α and β is identical to the equations analyzed by Erchova et al. [55] (their equations (A16) and (A17)). An alternative route to create a single equation from the two original ones is by direct substitution while leaving dimensionality intact. If we differentiate the first expression in (23) with respect to time we have CV¨+V˙/R+I˙L=I˙0. Here V¨, V˙, and I˙ indicate the second and first time derivatives of V and I. We can use the second expression in (23) to substitute I˙L=(V-RLIL)/L where IL as a function of V can be found from CV˙+V/R+IL=I0. This yield
(26)CV¨+V˙R+V-RL(I0-CV˙-V/R)L=I˙0→CV¨+(1R+RLCL)︸γV˙+1L(1+RLR)︸εV=RLLI0+I˙0.
This result is identical to the expression used by Erchova et al. [55] (their Equation (A1)). Interestingly this equation is also similar to the network equation by Jirsa and Haken [45] in which the spatial derivative is set to zero (see Section 4.5). This similarity across the scale from neuron to network is seen frequently. In one example, investigators even use a cellular model for action potential generation to represent activity of a population [56]. Although Curto et al. [56] do not explicitly describe this, in this case the inward and outward currents of the cell membrane can be considered to represent the network’s excitation and inhibition, respectively.


## 2.1.4. The Frequency Response of the Linearized Model

In electrophysiology, cellular properties are often probed with hyperpolarizing and depolarizing current steps. The unit impulse response, the derivative of the unit step, is the inverse Fourier transform or inverse Laplace transform of the system’s frequency response or transfer function, respectively. Since we deal with a linearized system with currentI as input and membrane potential V as output, the frequency response is just the impedance Z (i.e., Z(ω)=V(ω)/I(ω), with ω—frequency in rad/s). We can use the Laplace transform of (26) to obtain the transfer function H(s) of the linearized model: H(s)=V(s)/I(s)=(ε+γs+Cs2)/((RL/L)+s). In this form we can apply the Matlab freqs command to depict the Bode plot of the model (Figure 4).Bode plot of the linearized model depicted in Figure3. The parameters used to prepare the plots are C=100 pF; R=200 MΩ; RL=20 MΩ; L=2 MH. Note in the upper panel that a resonance peak is present at 75 rad/s.
(a)
(b)An interesting observation is that the Bode plot in Figure4 shows a peak in the magnitude of the impedance, indicating that neurons and their networks display resonance [55, 57–60]. With an appropriate, realistic choice of parameters, this peak is located  in the alpha-frequency range of the EEG at 75 rad/s (~12 Hz).


## 2.2. Models That Can Be Derived from the Hodgkin and Huxley Formalism

The H&H representation and its linearized version can be seen as a basis for many different neuronal models. Of course one could either add complexity to the H&H formalism or further reduce it. The following sections summarize representative examples of both types.


### 2.2.1. Complex Cell Models

Using the H&H formalism one might extend the sodium and potassium channels of the model neurons to also include different types of potassium channels (such as the A-current, a transient potassium current), persistent sodium channels, and different types of Ca++ channels. In addition to the voltage sensitive channels, it is common practice to also include Ca++ and/or Na+ sensitive channels. In addition to the addition of channels, the cell model can also be divided into many multiple coupled compartments (e.g., [4]). Each compartment represents a small part of the cell that can be approximated by a compartment with a uniform membrane potential. The rationale for the use of multiple equipotential compartments is to mimic the spatial separation of different functions: synaptic input, passive and active propagation of signals, and generation of action potentials. The more complex cell models are also frequently used to create networks in order to determine emergent behavior of neuronal populations (e.g., [3, 5, 6, 26, 61]). Special neural simulation packages such as Genesis and Neuron have been developed to simulate such models, even on parallel machines [62, 63]. A recent example in which many details are incorporated in the individual cell models is the Blue Brain Project [7]. The advantage of these models is that they are based on experimental data enabling the modeler to link the model’s behavior to measurements. This close relationship can also be a disadvantage of these complex models; due to the complexity, the simulation results can be as difficult to interpret as real recorded data.


### 2.2.2. Simplified Cell Models

One can simplify the H&H equations while preserving the nonlinearity, or one might use the linearization results described in Section2.1. In fact both these approaches have been employed. In one case, the four dimensional H&H model can be reduced to a lower dimensional one and in the other case, the equivalent circuit of a linearized membrane model (Figure 3) is coupled to an ad hoc spike generator (Figure 5). The advantage of preserving the nonlinearity is that the spike generation mechanisms remain (at least in part) intact, but mathematical analysis of such a model can be fairly complicated. The linearized cell models are more accessible for analysis but require an ad hoc spike generator (Table 1).Table 1
Overview of the dimensionality of cellular models.
Amplifying
Resonant
Membrane potential(V)
Inward current activation
Inward current inactivation
Outward current activation
Leak current capacitance(not time dependent)
Physiological AP mechanism
Hodgkin and Huxley  [1]
1
(m) 2
(h) 3
(n) 4
Yes
Yes
FitzHugh-Nagumo [9, 10]
1
Instantaneous(not time dependent)
Generic recovery
(w) 2
Yes
Yes
Morris and Lecar [12] (Ca instead of Na)
1
Instantaneous(not time dependent)
No
(w) 2
Yes
Yes
IF [13]
1
No
No
No
Yes
Generator
RF [15]
1
No
No
Simplified (L) 2
Yes
Generator
QIF [14]
1
Generic amplification(quadratic)
No
No
Yes
Reset
Izhikevich SMC [15]
1
Generic amplification(quadratic)
No
Simplified (L) 2
Yes
Reset
AP: action potential; L: inductor; IF: integrate and fire; QIF: quadratic integrate and fire; RF: resonate and Fire; SMC: simple model of choice;V, m, h, n: parameters of the Hodgkin and Huxley equations; w: recovery variable. Numbers 1–4 count the dimensions of the model.Figure 5
Example of an electronic circuit of an IF neuron model using analog components. The membrane is modeled by the R and C components (orange), the threshold is implemented by the OP-Amp as a comparator (the threshold is the value at the potentiometer, orange). The reset function is performed by the Reed relay. The diodes and 10 kΩ resistors at the in-and output are included to rectify the signals. This model is linear for the subthreshold part of the activity. (b) shows a measurement of the subthreshold activity of the membrane including the reset at the threshold; (c) depicts a measurement of the spike output.(1)FitzHugh-Nagumo/Morris and Lecar.The simplification used independently by FitzHugh [9], Nagumo et al., [10], and Morris and Lecar [12], is all based on the observation that the cell has a relatively fast activating inward current and a slower activating outward current. This allows one to simplify (2)–(5) from a four dimensional model into a model with only two dynamic variables. In the FitzHugh-Nagumo approach, the inward current is sodium whereas in the Morris and Lecar model for a barnacle muscle cell, the inward current is calcium. For the following discussion, it is not critical which ion species carries the inward current, so I will focus on simplification of the H&H formalism with a Na+ current. In the original equations (2)–(5) we have V, m, h, and n. The simplified version ignores activation of the inward current m because it is so fast it is considered instantaneous. The effects of h and n are similar; that is, they repolarize the membrane, so they can either be lumped or one of the two can be ignored so that only one repolarization parameter (often indicated by w) remains. This generates a nonlinear dynamical system with two variables. Although a few variations exist, the original equations in FitzHugh [9] were
(27)x˙=c(y+x-x33+z),y˙=-(x-a+by)c.
Here x, y, and z are membrane potential, recovery variable, and injected current, respectively; a, b, and c are constants.(2)Reduced Spiking Models.From the linearization procedure we can appreciate that the neural membrane can be represented by the RCL-circuit (Figure 3). Of course, this representation only works for studying subthreshold behavior such as integration of inputs or resonance properties. This type of study explains the propensity of some neurons to resonate with specific input frequencies. If, as an additional simplification, the ion channel properties are completely ignored, we may even further simplify the membrane to an RC-circuit; in this configuration there are no resonant properties. Both these simplifications, the RCL- and RC-circuits, ignore the nonlinear properties of the excitable membrane. In the so-called integrate-and-fire (IF) models the nonlinearity is again introduced, but in an ad hoc fashion by employing a firing threshold. When this threshold is exceeded, a spike + reset are pasted into the signal. Although such an ad-hoc firing process can be included in both RCL- and RC-circuits, it is commonly applied to the latter.Historically, the RC-based integrate-and-fire model was introduced long before the H&H formalism in 1907 by Lapicque [13]. Because of the leakage current in the resistor, it is also often called the leaky integrate-and-fire neuron (LIF). Much later a nonlinear variant of the LIF was introduced where the I-V relationship is not linear as in (23), but quadratic, the QIF model (e.g., [14]). Because this quadratic relationship, the I-V curve contains a stable and an unstable equilibrium, the spike’s upstroke can be generated by the model and only an ad-hoc spike reset is required (Figure 6).Figure 6
Phase space representation for the QIF neuron. The relationship of membrane current, proportional to the derivative ofV, is plotted against membrane potential, V. Due to the quadratic relationship, there are two fixed points: one stable (the resting potential) and one unstable one, the threshold, at the border between the subthreshold (blue) and superthreshold (red) regions. The model generates the upstroke of the action potential when V exceeds the threshold. To avoid instability, an ad hoc reset is provided when the arbitrary reset value Vmax is reached.Finally, one can combine the QIF scenario with the RCL-circuit to produce a model that is capable of both resonance and spike upstroke generation. This is defined as the simple model of choice (SMC) by Izhikevich and was shown to be able to produce a plethora of neuronal behavior (e.g., [15]).


## 2.2.1. Complex Cell Models

Using the H&H formalism one might extend the sodium and potassium channels of the model neurons to also include different types of potassium channels (such as the A-current, a transient potassium current), persistent sodium channels, and different types of Ca++ channels. In addition to the voltage sensitive channels, it is common practice to also include Ca++ and/or Na+ sensitive channels. In addition to the addition of channels, the cell model can also be divided into many multiple coupled compartments (e.g., [4]). Each compartment represents a small part of the cell that can be approximated by a compartment with a uniform membrane potential. The rationale for the use of multiple equipotential compartments is to mimic the spatial separation of different functions: synaptic input, passive and active propagation of signals, and generation of action potentials. The more complex cell models are also frequently used to create networks in order to determine emergent behavior of neuronal populations (e.g., [3, 5, 6, 26, 61]). Special neural simulation packages such as Genesis and Neuron have been developed to simulate such models, even on parallel machines [62, 63]. A recent example in which many details are incorporated in the individual cell models is the Blue Brain Project [7]. The advantage of these models is that they are based on experimental data enabling the modeler to link the model’s behavior to measurements. This close relationship can also be a disadvantage of these complex models; due to the complexity, the simulation results can be as difficult to interpret as real recorded data.


## 2.2.2. Simplified Cell Models

One can simplify the H&H equations while preserving the nonlinearity, or one might use the linearization results described in Section2.1. In fact both these approaches have been employed. In one case, the four dimensional H&H model can be reduced to a lower dimensional one and in the other case, the equivalent circuit of a linearized membrane model (Figure 3) is coupled to an ad hoc spike generator (Figure 5). The advantage of preserving the nonlinearity is that the spike generation mechanisms remain (at least in part) intact, but mathematical analysis of such a model can be fairly complicated. The linearized cell models are more accessible for analysis but require an ad hoc spike generator (Table 1).Table 1
Overview of the dimensionality of cellular models.
Amplifying
Resonant
Membrane potential(V)
Inward current activation
Inward current inactivation
Outward current activation
Leak current capacitance(not time dependent)
Physiological AP mechanism
Hodgkin and Huxley  [1]
1
(m) 2
(h) 3
(n) 4
Yes
Yes
FitzHugh-Nagumo [9, 10]
1
Instantaneous(not time dependent)
Generic recovery
(w) 2
Yes
Yes
Morris and Lecar [12] (Ca instead of Na)
1
Instantaneous(not time dependent)
No
(w) 2
Yes
Yes
IF [13]
1
No
No
No
Yes
Generator
RF [15]
1
No
No
Simplified (L) 2
Yes
Generator
QIF [14]
1
Generic amplification(quadratic)
No
No
Yes
Reset
Izhikevich SMC [15]
1
Generic amplification(quadratic)
No
Simplified (L) 2
Yes
Reset
AP: action potential; L: inductor; IF: integrate and fire; QIF: quadratic integrate and fire; RF: resonate and Fire; SMC: simple model of choice;V, m, h, n: parameters of the Hodgkin and Huxley equations; w: recovery variable. Numbers 1–4 count the dimensions of the model.Figure 5
Example of an electronic circuit of an IF neuron model using analog components. The membrane is modeled by the R and C components (orange), the threshold is implemented by the OP-Amp as a comparator (the threshold is the value at the potentiometer, orange). The reset function is performed by the Reed relay. The diodes and 10 kΩ resistors at the in-and output are included to rectify the signals. This model is linear for the subthreshold part of the activity. (b) shows a measurement of the subthreshold activity of the membrane including the reset at the threshold; (c) depicts a measurement of the spike output.(1)FitzHugh-Nagumo/Morris and Lecar.The simplification used independently by FitzHugh [9], Nagumo et al., [10], and Morris and Lecar [12], is all based on the observation that the cell has a relatively fast activating inward current and a slower activating outward current. This allows one to simplify (2)–(5) from a four dimensional model into a model with only two dynamic variables. In the FitzHugh-Nagumo approach, the inward current is sodium whereas in the Morris and Lecar model for a barnacle muscle cell, the inward current is calcium. For the following discussion, it is not critical which ion species carries the inward current, so I will focus on simplification of the H&H formalism with a Na+ current. In the original equations (2)–(5) we have V, m, h, and n. The simplified version ignores activation of the inward current m because it is so fast it is considered instantaneous. The effects of h and n are similar; that is, they repolarize the membrane, so they can either be lumped or one of the two can be ignored so that only one repolarization parameter (often indicated by w) remains. This generates a nonlinear dynamical system with two variables. Although a few variations exist, the original equations in FitzHugh [9] were
(27)x˙=c(y+x-x33+z),y˙=-(x-a+by)c.
Here x, y, and z are membrane potential, recovery variable, and injected current, respectively; a, b, and c are constants.(2)Reduced Spiking Models.From the linearization procedure we can appreciate that the neural membrane can be represented by the RCL-circuit (Figure 3). Of course, this representation only works for studying subthreshold behavior such as integration of inputs or resonance properties. This type of study explains the propensity of some neurons to resonate with specific input frequencies. If, as an additional simplification, the ion channel properties are completely ignored, we may even further simplify the membrane to an RC-circuit; in this configuration there are no resonant properties. Both these simplifications, the RCL- and RC-circuits, ignore the nonlinear properties of the excitable membrane. In the so-called integrate-and-fire (IF) models the nonlinearity is again introduced, but in an ad hoc fashion by employing a firing threshold. When this threshold is exceeded, a spike + reset are pasted into the signal. Although such an ad-hoc firing process can be included in both RCL- and RC-circuits, it is commonly applied to the latter.Historically, the RC-based integrate-and-fire model was introduced long before the H&H formalism in 1907 by Lapicque [13]. Because of the leakage current in the resistor, it is also often called the leaky integrate-and-fire neuron (LIF). Much later a nonlinear variant of the LIF was introduced where the I-V relationship is not linear as in (23), but quadratic, the QIF model (e.g., [14]). Because this quadratic relationship, the I-V curve contains a stable and an unstable equilibrium, the spike’s upstroke can be generated by the model and only an ad-hoc spike reset is required (Figure 6).Figure 6
Phase space representation for the QIF neuron. The relationship of membrane current, proportional to the derivative ofV, is plotted against membrane potential, V. Due to the quadratic relationship, there are two fixed points: one stable (the resting potential) and one unstable one, the threshold, at the border between the subthreshold (blue) and superthreshold (red) regions. The model generates the upstroke of the action potential when V exceeds the threshold. To avoid instability, an ad hoc reset is provided when the arbitrary reset value Vmax is reached.Finally, one can combine the QIF scenario with the RCL-circuit to produce a model that is capable of both resonance and spike upstroke generation. This is defined as the simple model of choice (SMC) by Izhikevich and was shown to be able to produce a plethora of neuronal behavior (e.g., [15]).


## 2.3. Summary of the Cell Models

We can consider the different models and summarize their relationships by counting the dimensions (Table1). The H&H model is four dimensional whereas the subthreshold component of the IF model is only one dimensional. Of course, dimensions are added when the ad hoc spike generator is included but most studies using this model will focus on the subthreshold behavior. Table 1 summarizes the model dimensions and provides an overview of the spiking mechanisms of the models we evaluated above.


## 3. Populations of Cell Models

Network models containing concrete nodes come in a wide variety, ranging from nodes represented by an on/off switch to a detailed multicompartmental cell model with multiple ion channels. Similarly, the connections between the nodes can be modeled by relatively simple rules or they can be simulated by biophysically realistic synaptic ion channels. In general, the purpose of these models is to provide simulation results that can be compared to a result from mathematical analysis and/or an experimental measurement.


### 3.1. On/Off Switches

The first simple network model using coupled switches was described by McCulloch and Pitts [16]. Their approach was to represent network elements as gates in a logical circuit (e.g., AND, OR, etc.). Others have used a similar approach and also compared the network on/off node with the spins in a spin-glass lattice (e.g., [17]). The so-called Ising spins can be up or down, indicated by +1 and −1, respectively [64]. Each spin Si tends to align with its neighbors. However, at high temperatures, the probability of this alignment decreases to 1/2. We model the cortical surface by a two dimensional lattice of Ising spins. The overall state of this lattice can be characterized by its magnetization m; the average of the individual spins: m=(1/N)∑i=1NSi. Below a certain critical temperature Tc, the lattice will retain global magnetization by an external magnetic field h. Above the critical temperature (also called the Curie temperature), thermal fluctuations will cause demagnetization. Although this is a mathematical framework in statistical mechanics, it has been used to represent neuronal networks where the spins are the neuronal nodes in the network; in these models, the up and down states are often presented by 1 (active) and 0 (inactive), respectively. These networks with nodes Si are typically simulated in discrete time t and follow the following update rules.(1)
At each time step, the state ofSi is determined by the state of its set (Sk) of its connected neighbors, reflected by function Hi=∑j∈Skwj,iSj, with wj,i representing a weighting function from node j to node i (in many models this weight is set to one).(2)
The update ofSi now takes place in a probabilistic fashion; the probability P that Si will be in the up state in the next time step t+1 is: P[Si(t+1)=1]=1/(1+e-Hi/T); variable T is the temperature; it can be seen that P=1/2 at T→∞ and that the update is a deterministic process (P=1) at T=0.This update process is a stochastic variant of one of the rules described by McCulloch and Pitts [16] where a neuron’s output y is updated by the state of its input from neurons xi in a deterministic fashion, that is, output y is set to 1 when its input exceeds a threshold and to 0 otherwise. A mean field approach for the analysis of the Ising spin model is described in Section 4.1.A model similar to the Ising spin lattice was described by Hopfield [18]. He also considered switches interconnected in a network in which the strength of the connections between any neuron pair i and j is determined by the synaptic weight wij of their connection (no unit has connection with itself wii=0 and, unrealistically, he also assumed that the weight matrix is symmetric, wij=wji). The principal difference between Hopfield’s approach and the ones described above is that the updates occur randomly and asynchronously; that is, for every time step a node is randomly selected and then an update rule is applied for that node only. This update rule is similar to the deterministic one employed by McCulloch and Pitts [16], that is, the selected neuron is active if the sum of activities of its connected neurons exceeds a threshold. More on stochastic approaches is discussed in Section 5.


### 3.2. Networks of Reduced Spiking Models

As described in Section2.2.2 (2), there is a variety of simplified models ranging from simple IF to rather complete SMC representations (Table 1). Especially the IF-model is frequently used for network simulations. One advantage of this model is that it is simply linear in between the spikes. By employing this approach, one can not only simulate the network but also, to some extent, analyze its properties mathematically. Using the simplification in Figure 3 and a further simplification by neglecting induction effects, we can simplify (23): CV˙=-V/R+IInject+Synapse. By simplifying this notation, one can describe the network dynamics by membrane potential vi, the activity level of each neuron, by
(28)v˙i(t)=A-vi(t)+Ei(t),
(e.g., [21, 22]). The injected current is separated into A and E(t). Parameter A represents spontaneous drive and E is the time-dependent input. The presynaptic spike train can be represented by a train of unit impulses ∑nδ(t-tn). Each incoming impulse generates a postsynaptic signal that can be modeled by a so-called alpha function; that is, for an impulse arriving at time tn, its postsynaptic effect is given by ∑nθ[t-tn](t-tn)α2e-α(t-tn). The Heaviside function θ is included because the alpha-function does not have a physical meaning for negative arguments, (t-tn)<0. Combining these considerations, one can conclude that input E is the weighted sum of all alpha functions resulting from the incoming spike train (e.g., [23, 24])
(29)Ei(t)=1ki∑nWj,iθ[t-tn](t-tn)α2e-α(t-tn).
Here we determine the input by using the synaptic weight matrix W, in which each entry wj,i is the connectivity strength of neuron j to neuron i. Scaling factor ki represents the number of presynaptic inputs of neuron i. We now determine the Laplace transform of (29) and rearrange terms (s2+2αs+α2)Ei(t)=(α2/ki)∑nWj,ie-tns, which can also be inverse transformed into the time domain in the form of a differential equation
(30)E¨i(t)+2αE˙i(t)+α2Ei(t)=α2ki∑nWj,iδ(t-tn).
Since the first order differential equation (28) and the second order one (30) are coupled, via input E, they create a three dimensional nonlinear system; a system capable of exhibiting oscillatory and chaotic behavior. To study the activity in the model, Olmi et al. [23] rewrite the system as three 1st-order differential equations and integrate between pulses tn and tn+1. This creates an event-driven map that can be further investigated. Tattini et al. [24] examine the propensity for chaotic behavior by evaluating the attractors characterizing the network activities and they show that both connectivity and network size are critical properties—unsurprisingly, the level of chaos, reflected by the maximum Lyapunov exponent, decreases with network size, although less so for sparsely connected ones.More complicated cell models have also been employed to investigate network activity. Recently, Izhikevich and Edelman [25] used the SMC model to create a large network while using connectivity based on the tracts determined by a diffusion tensor image of the brain.


### 3.3. Networks of Ion Channel Models Based on Biophysical Considerations

Using the H&H formalism combined with approximation of the cell morphology by segments, one can create realistic representations of neurons. Addition of ligand-sensitive, synaptic channels allows one to create networks. Equation (2) is then adapted to include currents associated with other voltage sensitive channels, the synaptic currents, and intercompartmental currents
(31)CdVdt=-ILeak-IV-Sensitive-ISynapse+IIntercompartmental+IInject.
Interestingly, by applying Kirchhoff’s first law, it is apparent that intracellular intercompartmental currents must be equal to the extracellular current between the compartments, thereby providing a tool to model extracellular measurements. The dynamics of the voltage-sensitive channels is governed by expressions similar to (3)–(5). The intercompartmental current, following Ohm’s law, is determined by the potential difference between the isopotential compartments divided by the axial resistance between them. Traub and co-workers were amongst the first to develop this approach (e.g., [2, 26]), the blue brain project [7] is a famous recent attempt to determine properties that govern network activity at mesoscopic levels. The advantage of the approach is that the neurons in the network are modeled with realistic ion channels, which allows a direct comparison with experimental data. However, well-known problems with the approach are associated with the complexity and large number of parameters, many of which are not experimentally determined yet. The large parameter space renders this type of model ill-posed and the complexity can make interpretation of the results very difficult. A combined approach with associated detailed and global models can be very useful to understand the principles underpinning the simulation results of the detailed models [5, 6, 27, 65, 66].The detailed model of neocortex outlined in Figure7(a) shows pyramidal cell populations and the inhibitory basket and chandelier cells. The network generates an EEG signal based on the weighted sum of the neuronal currents. Surprisingly, during reduced excitatory coupling, the model starts to oscillate (top trace in Figures 7(a) and 7(b)) [5]. The mechanism generating these oscillations is difficult to grasp when studying the individual rasters of the superficial and deep pyramidal neurons (grey panels in Figure 7(b)). However when lumping each of the superficial and deep pyramidal activities (Figure 7(c)), it can be seen that the simulated EEG oscillation is associated with alternating oscillations in the two populations. The mechanism of the EEG oscillation can be much better understood if one considers population activity patterns while ignoring the details of the individual neural activity. If one considers individual neuron’s activity as depicted in the rasters in Figure 7(b), one might not see the relationship between neuron and network activity since only few neurons are active during each cycle; that is, most neurons show cycle skipping (Wallace et al., 2012).Model of Neocortex. The detailed model (after [6]) is outlined in diagram (a). The excitatory populations consist of superficial (S_PYR) and deep (D_PYR) pyramidal cells and the inhibitory neuron type (I) consists of populations of basket cells (B) and chandelier cells (CH). The weighted sum of all soma currents is used as the “EEG” signal (top, (a)). The excitatory connections of the microcircuitry are shown in the left part, and the inhibitory synapses in the right part of (a). (b) depicts a detail of the EEG during seizure-like oscillation in the “EEG” and the associated raster plots of the S_PYR and D_PYR populations. The simplified network version, after Visser et al. [65], where only excitatory pyramidal cells and inhibitory neuronal populations are considered, is shown in the right part of (b). (c) The average membrane potentials of both pyramidal cell populations show that the activities of the deep S_PYR and D_PYR populations alternate during seizure-like oscillation.
(a)
(b)
(c)Not only the oscillatory activity pattern itself but also the sudden transition between nonoscillatory and oscillatory network behavior can be understood from bifurcation analysis of a global model of two coupled excitatory-inhibitory populations as depicted in Figure7(b) [65, 66].


## 3.1. On/Off Switches

The first simple network model using coupled switches was described by McCulloch and Pitts [16]. Their approach was to represent network elements as gates in a logical circuit (e.g., AND, OR, etc.). Others have used a similar approach and also compared the network on/off node with the spins in a spin-glass lattice (e.g., [17]). The so-called Ising spins can be up or down, indicated by +1 and −1, respectively [64]. Each spin Si tends to align with its neighbors. However, at high temperatures, the probability of this alignment decreases to 1/2. We model the cortical surface by a two dimensional lattice of Ising spins. The overall state of this lattice can be characterized by its magnetization m; the average of the individual spins: m=(1/N)∑i=1NSi. Below a certain critical temperature Tc, the lattice will retain global magnetization by an external magnetic field h. Above the critical temperature (also called the Curie temperature), thermal fluctuations will cause demagnetization. Although this is a mathematical framework in statistical mechanics, it has been used to represent neuronal networks where the spins are the neuronal nodes in the network; in these models, the up and down states are often presented by 1 (active) and 0 (inactive), respectively. These networks with nodes Si are typically simulated in discrete time t and follow the following update rules.(1)
At each time step, the state ofSi is determined by the state of its set (Sk) of its connected neighbors, reflected by function Hi=∑j∈Skwj,iSj, with wj,i representing a weighting function from node j to node i (in many models this weight is set to one).(2)
The update ofSi now takes place in a probabilistic fashion; the probability P that Si will be in the up state in the next time step t+1 is: P[Si(t+1)=1]=1/(1+e-Hi/T); variable T is the temperature; it can be seen that P=1/2 at T→∞ and that the update is a deterministic process (P=1) at T=0.This update process is a stochastic variant of one of the rules described by McCulloch and Pitts [16] where a neuron’s output y is updated by the state of its input from neurons xi in a deterministic fashion, that is, output y is set to 1 when its input exceeds a threshold and to 0 otherwise. A mean field approach for the analysis of the Ising spin model is described in Section 4.1.A model similar to the Ising spin lattice was described by Hopfield [18]. He also considered switches interconnected in a network in which the strength of the connections between any neuron pair i and j is determined by the synaptic weight wij of their connection (no unit has connection with itself wii=0 and, unrealistically, he also assumed that the weight matrix is symmetric, wij=wji). The principal difference between Hopfield’s approach and the ones described above is that the updates occur randomly and asynchronously; that is, for every time step a node is randomly selected and then an update rule is applied for that node only. This update rule is similar to the deterministic one employed by McCulloch and Pitts [16], that is, the selected neuron is active if the sum of activities of its connected neurons exceeds a threshold. More on stochastic approaches is discussed in Section 5.


## 3.2. Networks of Reduced Spiking Models

As described in Section2.2.2 (2), there is a variety of simplified models ranging from simple IF to rather complete SMC representations (Table 1). Especially the IF-model is frequently used for network simulations. One advantage of this model is that it is simply linear in between the spikes. By employing this approach, one can not only simulate the network but also, to some extent, analyze its properties mathematically. Using the simplification in Figure 3 and a further simplification by neglecting induction effects, we can simplify (23): CV˙=-V/R+IInject+Synapse. By simplifying this notation, one can describe the network dynamics by membrane potential vi, the activity level of each neuron, by
(28)v˙i(t)=A-vi(t)+Ei(t),
(e.g., [21, 22]). The injected current is separated into A and E(t). Parameter A represents spontaneous drive and E is the time-dependent input. The presynaptic spike train can be represented by a train of unit impulses ∑nδ(t-tn). Each incoming impulse generates a postsynaptic signal that can be modeled by a so-called alpha function; that is, for an impulse arriving at time tn, its postsynaptic effect is given by ∑nθ[t-tn](t-tn)α2e-α(t-tn). The Heaviside function θ is included because the alpha-function does not have a physical meaning for negative arguments, (t-tn)<0. Combining these considerations, one can conclude that input E is the weighted sum of all alpha functions resulting from the incoming spike train (e.g., [23, 24])
(29)Ei(t)=1ki∑nWj,iθ[t-tn](t-tn)α2e-α(t-tn).
Here we determine the input by using the synaptic weight matrix W, in which each entry wj,i is the connectivity strength of neuron j to neuron i. Scaling factor ki represents the number of presynaptic inputs of neuron i. We now determine the Laplace transform of (29) and rearrange terms (s2+2αs+α2)Ei(t)=(α2/ki)∑nWj,ie-tns, which can also be inverse transformed into the time domain in the form of a differential equation
(30)E¨i(t)+2αE˙i(t)+α2Ei(t)=α2ki∑nWj,iδ(t-tn).
Since the first order differential equation (28) and the second order one (30) are coupled, via input E, they create a three dimensional nonlinear system; a system capable of exhibiting oscillatory and chaotic behavior. To study the activity in the model, Olmi et al. [23] rewrite the system as three 1st-order differential equations and integrate between pulses tn and tn+1. This creates an event-driven map that can be further investigated. Tattini et al. [24] examine the propensity for chaotic behavior by evaluating the attractors characterizing the network activities and they show that both connectivity and network size are critical properties—unsurprisingly, the level of chaos, reflected by the maximum Lyapunov exponent, decreases with network size, although less so for sparsely connected ones.More complicated cell models have also been employed to investigate network activity. Recently, Izhikevich and Edelman [25] used the SMC model to create a large network while using connectivity based on the tracts determined by a diffusion tensor image of the brain.


## 3.3. Networks of Ion Channel Models Based on Biophysical Considerations

Using the H&H formalism combined with approximation of the cell morphology by segments, one can create realistic representations of neurons. Addition of ligand-sensitive, synaptic channels allows one to create networks. Equation (2) is then adapted to include currents associated with other voltage sensitive channels, the synaptic currents, and intercompartmental currents
(31)CdVdt=-ILeak-IV-Sensitive-ISynapse+IIntercompartmental+IInject.
Interestingly, by applying Kirchhoff’s first law, it is apparent that intracellular intercompartmental currents must be equal to the extracellular current between the compartments, thereby providing a tool to model extracellular measurements. The dynamics of the voltage-sensitive channels is governed by expressions similar to (3)–(5). The intercompartmental current, following Ohm’s law, is determined by the potential difference between the isopotential compartments divided by the axial resistance between them. Traub and co-workers were amongst the first to develop this approach (e.g., [2, 26]), the blue brain project [7] is a famous recent attempt to determine properties that govern network activity at mesoscopic levels. The advantage of the approach is that the neurons in the network are modeled with realistic ion channels, which allows a direct comparison with experimental data. However, well-known problems with the approach are associated with the complexity and large number of parameters, many of which are not experimentally determined yet. The large parameter space renders this type of model ill-posed and the complexity can make interpretation of the results very difficult. A combined approach with associated detailed and global models can be very useful to understand the principles underpinning the simulation results of the detailed models [5, 6, 27, 65, 66].The detailed model of neocortex outlined in Figure7(a) shows pyramidal cell populations and the inhibitory basket and chandelier cells. The network generates an EEG signal based on the weighted sum of the neuronal currents. Surprisingly, during reduced excitatory coupling, the model starts to oscillate (top trace in Figures 7(a) and 7(b)) [5]. The mechanism generating these oscillations is difficult to grasp when studying the individual rasters of the superficial and deep pyramidal neurons (grey panels in Figure 7(b)). However when lumping each of the superficial and deep pyramidal activities (Figure 7(c)), it can be seen that the simulated EEG oscillation is associated with alternating oscillations in the two populations. The mechanism of the EEG oscillation can be much better understood if one considers population activity patterns while ignoring the details of the individual neural activity. If one considers individual neuron’s activity as depicted in the rasters in Figure 7(b), one might not see the relationship between neuron and network activity since only few neurons are active during each cycle; that is, most neurons show cycle skipping (Wallace et al., 2012).Model of Neocortex. The detailed model (after [6]) is outlined in diagram (a). The excitatory populations consist of superficial (S_PYR) and deep (D_PYR) pyramidal cells and the inhibitory neuron type (I) consists of populations of basket cells (B) and chandelier cells (CH). The weighted sum of all soma currents is used as the “EEG” signal (top, (a)). The excitatory connections of the microcircuitry are shown in the left part, and the inhibitory synapses in the right part of (a). (b) depicts a detail of the EEG during seizure-like oscillation in the “EEG” and the associated raster plots of the S_PYR and D_PYR populations. The simplified network version, after Visser et al. [65], where only excitatory pyramidal cells and inhibitory neuronal populations are considered, is shown in the right part of (b). (c) The average membrane potentials of both pyramidal cell populations show that the activities of the deep S_PYR and D_PYR populations alternate during seizure-like oscillation.
(a)
(b)
(c)Not only the oscillatory activity pattern itself but also the sudden transition between nonoscillatory and oscillatory network behavior can be understood from bifurcation analysis of a global model of two coupled excitatory-inhibitory populations as depicted in Figure7(b) [65, 66].


## 4. Mean Field Network Models

The network models discussed in this section do not explicitly simulate each network node, but start directly by modeling populations of nerve cells. Some of the models explicitly use the mean field approach commonly employed in statistical mechanics, others create neural populations with functionality inspired by the single neuron function. The final result of these approaches can be fairly similar since they both usually contain excitatory and inhibitory components (the Ising spin model is an exception). As I will point out repeatedly, this results in a strong similarity between the equations that result from the various network models and even between these and the equations describing single neurons.


### 4.1. The Ising Spin Model

Two- or three-dimensional models of the Ising-spin lattice, where the spins represent the neurons, are not easy to solve. A mean field approach is the simplest approximation to describe such a system. In the lattice, each node experiences the sum of the magnetic forces created by the other nodes and, if present, an external magnetic field. In the mean field approach, we replace the internal field generated by the spins by its average value. It is beyond the scope of this paper to derive the mean field relationship, which can be found in many physics textbooks (e.g., [67, Chapter 7]). Using the mean field expression, we can describe the equilibrium state of the lattice of spins—or network of model neurons—by the expression
(32)(Jnm+h)-kTtanh-1(m)=0,
with k—Boltzmann constant (in the following k is set to unity), T—temperature, J—magnetic coupling strength, n—number of neighbors, m—magnetization of the lattice, and h—an externally applied magnetic field. As shown in the following, (32) shows sudden transitions in behavior according to the so-called cusp catastrophe [68]. First we can determine the fixed points for the expression graphically by finding the intersections of both terms, Ttanh-1(m) and Jnm+h as a function of m. When depicting these functions for different T and h, it can be seen that they intersect at either 3 or at 1 point(s). In case of the intersection at 3 points, the middle one is an unstable fixed point, the other two are stable ones; therefore the system is bistable in one area of the T-h plane and monostable elsewhere.The transition between bi- and mono-stability occurs at the bifurcation where the straight lineJnm+h intersects the sigmoid Ttanh-1(m) tangentially; this indicates that two conditions must be satisfied at the transition in point Tc, h*:
(33)Tctanh-1(m)=Jnm+h*,
and since they must be tangent
(34)ddm[Tctanh-1(m)]=ddm[Jnm+h*]→Tc1-m2=Jn.
From these two expressions, one can find the expression for Tc and h* (as a function of m). The plot of these functions in the T-h plane is the so-called stability diagram (Figure 8(a)). These lines mark the border where the system transitions from mono- (1 fixed point) to bi-stability (3 fixed points with one of the points unstable). Another common presentation is a bifurcation diagram; here we have a codimension-2 case with two independent variables T and h. Note that for h=0, m=0 is always a solution (fixed point) for (32). If we use Jn=1, we have a stable fixed point at m=0 for T>1 and an unstable one plus two stable fixed points for T<1. This is depicted in Figure 8(b); we can see that we have the characteristic of a supercritical pitchfork bifurcation. This plot can be considered as cross-sections of a 3D plot of m(h,T) shown in Figure 8(c); this surface depicts the so-called cusp catastrophe.Properties of the mean field equation of the Ising model. (a) The stability diagram. (b) Bifurcation diagram. (c) The cusp catastrophe shows the dependence of magnetizationm on both the temperature (T) and an external magnetic field (h).
(a)
Stability diagram
(b)
Bifurcation diagram ath=0
(c)
Cusp catastrophe


### 4.2. The Wilson and Cowan Model

The 1st model by Wilson and Cowan [29] is exclusively temporal and describes the behavior of two neuronal populations—one excitatory and one inhibitory—embedded in neural tissue. It is more realistic than using the Ising model that only considers the equivalent of excitation. In the Wilson-Cowan model, population properties are derived from assumed statistics of the neurons [29]. Each population receives input from itself, the other population, and external sources. The input to a population affects the neurons of that population that are available (the inactive proportion of the population). The model is nonlinear because the incoming activity for the population is related to its output through a nonlinear sigmoid function S. In a 2nd version of the model that will not be further discussed here, Wilson and Cowan [30] extended their model to include the spatial component as well. Instead of formulating the model from underlying statistical properties, a coarse grained version of the model can be directly derived from the diagram in Figure 9. Two populations, the excitatory and inhibitory, are coupled with strengths c1–c4. The external population affects the excitatory and inhibitory networks with activities P and Q. We denote excitatory activity by E and inhibitory activity by I, the excitatory population in Figure 9 receives input that can be quantified as c1E-c2I+P. This input of the population is converted to an output expression via a sigmoid function Se{⋯}.  In most followers of their approach, such a sigmoid curve is employed to represent conversion from membrane potential to spike rate and/or vice versa. In the Wilson and Cowan’s approach, the sigmoid curve represents the cumulative distribution of the probability that a neuron will fire at a given input, that is, the distance between membrane potential and firing threshold. Further, assuming that the unavailable part of the population, that is, the active + refractory portion is proportional (by a fraction re) to the activity level E then the receptive population (available portion of E) is 1-reE.Figure 9
Diagram of the Wilson-Cowan Model containing excitatory and inhibitory populations. The coupling strength between the populations is denoted with constantsc1–c4 and external activity is symbolized by P and Q. Note that the model contains nonlinear functions to relate input and output for each population (not shown in the diagram).In Wilson and Cowan [29], the authors use a value a little bit less than unity for ke, reflecting that the maximum values for the sigmoid function Se will be a little bit less than one (see pages 9-10 in the original paper). Therefore, the available population becomes ke-reE. The formalism for population activity also must describe spontaneous decay proportional to E, reflecting that a nonstimulated group of neurons converge to a base level of activity (zero, reflecting a low level of spontaneous activity in this model). Finally all changes dE/dt are weighted by the time constant of the neuronal population τe. Putting this all together the coarse grained equation for excitatory activity becomes
(35)τe︸EXCITATORYTIME CONSTANTdEdt=-E︸DECAY+(ke-reE)︸AVAILABLEPOPULATIONSe{c1E-c2I+P}︸.INPUT︷SIGMOID
A similar approach for the inhibitory population yields
(36)τidIdt=-I+(ki-riI)Si{c3E-c4I+Q}.
In their 1972 paper [29], Wilson and Cowan derive these relationships from first principles followed by a coarse graining procedure (their equations (11) and (12)). When looking into the two dimensional nonlinear Wilson-Cowan model, there should be a variety of behaviors including a steady activity scenario, oscillatory activity with damping, and limit cycles (Figure 10).Phase plane representations of the Wilson-Cowan model during equilibrium, oscillatory, and limit-cycle activities. These plots correspond to the parameters in Figures 4, 10, and 11 of the original 1972 paper [29]. Note that in the equilibrium case there are two stable fixed points: one in a “down” and one in an “up” state. Also note that in these examples, oscillations are obtained by flattening the E null-cline (magenta) and making the I null-cline (orange) steeper. Sample trajectories are depicted in blue. (a) Three fixed points, two of them are stable and one (in between) is unstable. The eigenvalues of these points from left to right are as follows: −0.59, −1.13; 0.73, −1.6; −1.43, −2.88. (b) Two fixed points, one is clearly associated with a damped oscillation. The eigenvalues from left to right are as follows: −0.34, −0.98; -0.09±0.86j. (c) The limit cycle with a fixed point inside. The eigenvalues of the fixed point are 0.12±1.88j.
(a)
(b)
(c)


#### 4.2.1. Linearization of the Wilson-Cowan Equations

To study resonance phenomena, we will consider the effect onE and I by small sinusoidal perturbations (e.g., evoked by electrical stimulation) close to an equilibrium state. When stimulating electrically in the extracellular space and assuming similar spatial distribution of excitatory and inhibitory cells, it is reasonable to assume that both the excitatory and inhibitory populations receive the same stimulus current. Thus the external inputs P and Q are considered equal. Further simplifying the sigmoid curves Se and Si to lines with slopes α and β, the linearized versions of (35) and (36) become
(37)E˙=-Eτe+(ke-reE)τeα(c1E-c2I+P),I˙=-Iτi+(ki-riI)τiβ(c3E-c4I+P).
Rewriting these results and using ατ=α/τe and βτ=β/τi, we get
(38)E˙=f(E,I,P)=-Eτe+keατc1E-keατc2I+keατP-reατc1E2+reατc2EI-reατEP,I˙=g(E,I,P)=-Iτi+kiβτc3E-kiβτc4I+kiβτP-riβτc3EI+riβτc4I2-riβτIP.
We now further linearize and consider small perturbations δE, δI, δP around the equilibrium indicated by *(39)δE˙=(∂f∂E)*δE+(∂f∂I)*δI+(∂f∂P)*δP,δI˙=(∂g∂E)*δE+(∂g∂I)*δI+(∂g∂P)*δP.
The partial derivatives determined at equilibrium * are:
(40)εE=(∂f∂E)*=(-1τe+keατc1-2reατc1E+reατc2I-reατP)*,εI=(∂f∂I)*=(-keατc2+reατc2E)*,εP=(∂f∂P)*=(keατ-reατE)*,γE=(∂g∂E)*=(kiβτc3-riβτc3I)*,γI=(∂g∂I)*=(-1τi-kiβτc4-riβτc3E+2riβτc4I-riβτP)*,γP=(∂g∂P)*=(kiβτ-riβτI)*.
Combining these expressions and substituting e and i for δE and δI, and making the perturbation δP sinusoidal psinωt (where amplitude p is a small value), we can simplify notation
(41)e˙=εEe+εIi+εPpsinωt,(42)i˙=γEe+γIi+γPpsinωt.
These two 1st order differential equations can also be written as a single 2nd order one. Taking the derivative of (42) and substituting the expression in (41) for i gives
(43)e¨=εEe˙+εI(γEe+γIi+γPpsinωt)+εPωpcosωt.
Using (41), we can substitute i=(e˙-εEe-εPpsinωt)/εI, and we get the standard expression for a 2nd order differential equation with a driving force
(44)e¨-(εE+γI)︸be˙+(εEγI-εIγE)︸ce=εPpdsinωtdt+(εIγP-εPγI)psinωt︸driving"force".
This equation is again of the same form as the equation that governs a parallel RCL circuit as shown in Figure 3. The 2nd-order ODE governing this RCL circuit is (26). This equation is also very similar to the network equation by Jirsa and Haken [45] in which the spatial derivative is set to zero (Section 4.5). Thus the linearized network model is capable of resonance behavior, described in Section 2.1.4. Since we deal with a linearized version of a nonlinear system, unstable solutions of the linearized case may relate to bifurcations in the full nonlinear version of the model (see Section 6.1 for further discussion).


### 4.3. Models Inspired by the Wilson-Cowan Approach

The model described by Lopes da Silva et al. [35] is a model of thalamus and one of the first that  is inspired by the Wilson-Cowan approach. A similar treatment of the mesoscopic model, with altered synaptic function, was given by van Rotterdam et al. [36]. The base model describes a thalamic network; however, the populations are identical to the ones in the cortical Wilson-Cowan model: one excitatory and one inhibitory cell group. Similar to the Wilson-Cowan approach they assume an invertible function f that translates membrane potential, Ve(t) and Vi(t)  (i.e., the input) for the excitatory and inhibitory populations, into a rate of action potentials for each population E(t) and I(t) (i.e., the output). For the excitatory population we then get the following relationships: E(t)=fe[Ve(t)]  and Ve(t)=fe-1[E(t)]. Here fe and fe-1 are the function and its inverse that translates membrane potential into spike rate and vice versa. Both the functions fe and its equivalent for the inhibitory population fi are basically the static nonlinearities in their model. Each nonlinear component is subsequently approached by a Taylor series about the mean value of E denoted as E-: fe-1[E(t)]=ae0E-+ae1(E(t)-E-)+ae2(E(t)-E-)2+⋯. Subsequently Lopes da Silva et al. simplify notation by only considering the deviation from the mean only, e(t)=E(t)-E-. This is done for all variables in the model, and next they linearize about the mean by neglecting the 2nd and higher order terms in the Taylor series so that the original equations can be approximated by ve(t)=fe-1[e(t)]=ae1e(t) and vi(t)=fi-1[i(t)]=ai1i(t). External input P(t) or its demeaned version p(t) perturbs the system (top Left in Figure 11(a)). For the excitatory component in the diagram in the Laplace domain we get Ve(s)=He(s)P(s), with the following Laplace transform pairs: Ve(s)⇔ve(t),  He(s)⇔he(t), and P(s)⇔p(t). Here he and He are the unit impulse response and transfer function of the excitatory network.(a) Diagram of the Model (after [35]). P(t) is the external input; he and hi are the population’s linear unit impulse responses reflecting the synaptic function. Variables Ve and E are the population’s mean membrane potential and firing rate respectively; Vi and I are the same for the inhibitory population. Functions fe and fi are the nonlinear functions that convert population membrane potential into population firing rate. C1 and C2 are coupling constants. (b) Result for the spectrum following (49).
(a)
(b)In a similar fashion we include the effect of the inhibition in the diagram in Figure11(a) and get
(45)Ve(s)=ae1E(s)=He(s)P(s)-C2I(s)Hi(s).
Following the same procedure for the linearized inhibitory component (bottom part in Figure 11(a)) gives
(46)Vi(s)=ai1I(s)=C1E(s)He(s).
If we substitute the expression for I(s) (obtained from (46)) into (45) and use E(s)=Ve(s)/ae1 (from (45)), we obtain the following expression for Ve(s):
(47)Ve(s)=He(s)P(s)1+C1C2Hi(s)He(s)/ai1ae1.
Then the authors use a linear system to model synaptic input with unit impulse responses (UIRs) he and hi. In Lopes da Silva et al. [35], the synaptic UIR is a dual exponential A[e-a1t-e-a2t] and in van Rotterdam et al. [36] it is an alpha-function Ate-at; in both cases, these functions are only considered for t≥0.The associated Laplace and Fourier transforms are, of course, slightly different. Here we follow the paper of Lopes da Silva et al. [35] and use the dual exponential unit impulse response for both the excitatory and inhibitory populations. These functions and their Laplace transforms are
(48)he(t)=A[e-a1t-e-a2t]⟺He(s)=A(a2-a1)(s+a1)(s+a2),hi(t)=B[e-b1t-e-b2t]⟺Hi(s)=B(b2-b1)(s+b1)(s+b2).
We can substitute this into (47) and simplify
(49)Ve(s)=A(a2-a1)(s+b1)(s+b2)(s+a1)(s+a2)(s+b1)(s+b2)+KP(s),
with K=(C1C2/ai1ae1)(a2-a1)(b2-b1)AB.We can substitutes by jω in Ve(s) to obtain an expression for the spectrum of ve. The spectrum can be further quantified by assuming that the thalamic input p(t) is white noise and therefore P(ω) is constant. The numerical values we use for Figure 11(b) are the same as the ones in the example given in Lopes da Silva et al. [35]: A=1.65 mV, B=32 mV, C1=32 cells, C2=3 cells, a1=55 s−1, a2=605 s−1, b1=27.5 s−1, b2=55 s−1, qe1=1/ae1, qi1=1/ai1, and qi1qe1=4.55×106 (K≈3.5e8). It can be seen that the spectrum peaks around 10 Hz, the alpha rhythm frequently encountered in the EEG of subjects with eyes closed.


#### 4.3.1. Neural Mass Models

The Wilson-Cowan Equations [29] and Lopes da Silva et al. approach [35] inspired many other models. Instead of considering individual neurons, these models determine the activity of coupled populations of nerve cells, each population representing a neural mass. A very basic model, using two populations representing a coupled excitatory and inhibitory network is described by Dayan and Abbott [69, pages 265–270]. The rate of change for each network is inversely proportional to its activity level, and proportional to the input it receives from itself, the other population, and an external input. This results in a set of two coupled differential equations. Although the model seems to be linear, their model is actually nonlinear because the population activities are limited to positive values only. Therefore their model is capable of generating limit cycles.Jansen’s neural mass model [37] is applied by Grimbert and Faugeras [38] to create a cortical network unit (Figure 12(a)). It is a temporal model of neuronal populations inspired by the Lopes da Silva et al. [35] approach where each population is represented by a Wiener cascade ([70], Chapters 2–5): a dynamic linear component, modeling the synaptic conversion process, and a nonlinear static component that represents the conversion of membrane potential to spike rate (Figure 12(b)). As it can be seen in Figure 12(b), input x1 (input spike rate) is first convolved with the unit impulse response (UIR) h of the linear dynamic component, the result x2=x1⊗h (population membrane potential; ⊗-convolution) is subsequently converted with a static sigmoid to produce output x3=σ(x2) (output spike rate). The structure of the model by Grimbert and Faugeras [38] is similar to the canonical cortical model of Douglas and Martin [71] in the sense that it also includes two excitatory populations E, one inhibitory population I, and an external input P. The development of the equations that govern the activities is straightforward. First we use the UIR to write a convolution expression. Again, the synapse is considered a linear system and in this model, the synaptic impulse response function is defined as an alpha function (as in [36]):
(50)h(t)={αβte-βt,t≥0,0,otherwise,
in which α and β are population-specific constants. The Laplace transform of the unit impulse response is: H(s)=αβ/(s+β)2. In the time/Laplace domain we may link the input x(t)/transformed input X(s) to the output y(t)/transformed output Y(s) of any linear system by the transform pair: Y(s)=H(s)X(s)⇔y(t)=h(t)⊗x(t). Second, we write the convolution relationship as an ODE. Substituting the above expression for H in the Laplace transformed expression and we get (s+β)2Y(s)=αβX(s). This result can be re written as s2Y(s)+2βsY(s)+β2Y(s)=αβX(s). Recall that sY(s) and s2Y(s) denote the Laplace transforms of the 1st and 2nd derivatives of y(t) in the time domain; using this, we can transform the above expression into a 2nd order ODE for the linear synaptic process: y¨(t)=αβx(t)-2βy˙(t)-β2y(t), which can be written as two 1st order ODEs
(51)y˙(t)=z(t),z˙(t)=αβx(t)-2βz(t)-β2y(t).
Now we add the static nonlinearity employed for the membrane potential to pulse rate conversion, a sigmoid curve σ(v) (Figure 12(b)). In this study, it is defined as a Boltzmann function with maximum rate vmax, slope r and a 50% level threshold v0: σ(v)=vmax/(1+e-r(v-v0)). This completes the formalism for one population, a single Wiener cascade in Figure 12(b). The pair of ODEs plus the static nonlinearity can be defined for each of the three populations in the neural mass model (two excitatory and one inhibitory, Figure 12(a)). This means we will have a pair of 1st order differential equations per population, generating a total of six equations (Equation (3) in [38])
(52)y˙0(t)=y3(t),y˙3(t)=Aaσ(y1(t)-y2(t))-2ay3(t)-a2y0(t),y˙1(t)=y4(t),y˙4(t)=Aa{p(t)+C2σ(C1y0)}-2ay4(t)-a2y1(t),y˙2(t)=y5(t),y˙5(t)=BbC4σ(C3y0)-2by5(t)-b2y2(t).
The parameters y0–y2 denote the membrane potentials, and A, a, B, b are the constants in the population specific alpha functions. Note that in these equations, some outputs y and sigmoid functions are weighted by connectivity strengths C1–C4. The focus of the paper is on the variable y, the aggregate membrane potential of the main population, which is thought to be proportional with the local field potential. A simulation using the following parameters, A=3.25, a=100, B=22, b=50, C1=135, C2=0.8*C1, C3=0.25*C1, C4=0.25*C1, vmax=5, r=0.56, and v0=6, shows a signal that is similar to the EEG alpha rhythm (Figure 12(a)).Neural mass models. (a) A diagram of the model of a cortical unit employed by Grimbert and Faugeras [38] and a sample trace showing an oscillation similar to an EEG alpha rhythm. (b) Representation of a Wiener cascade used to model a population of nerve cells. The Wiener cascade consists of a linear dynamical module and a static nonlinear one. (c) Schematic representation of the model described by Freeman [39] and sample traces showing chaotic behavior of the impulse responses (impulse input—arrows).
(a)
(b)
(c)


#### 4.3.2. Modeling the Olfactory System

The models created by Freeman in the 1980s used a combined network modeling and experimental approach to study the electrical activity in rabbit olfactory system. Although Freeman’s approach may not have been directly inspired by the Wilson-Cowan equations and/or the models by Lopes da Silva and co-workers, the fundamental building blocks of his model are neural masses. The appeal of Freeman’s approach is that it is strongly based on the anatomy of the olfactory system, and that he relates modeling results to electrophysiological experiments. As compared to the previous models, a unique addition is that it contains latencies L1–L4, modeling the conduction delays occurring in between populations [39]. Differential equations without delays need initial values, whereas equations with delays require initial functions to determine the delayed response of the system. Equations with delays can exhibit behavior that cannot be displayed by equations without delay of the same order. It is given that delays in the nervous system occur because of the finite conduction velocity of nerve fiber activity. In addition to Freeman’s model, several others have explored the role of these delays or included the delays in the formalism (e.g., [65, 66, 72–74]).In his models, Freeman describes different population levels. Within the most basic population type, which is defined as K0, he includes pulse-to-wave conversion, linear spatiotemporal integration, wave-to-pulse conversion, and conduction of the output. This basic population is modeled with a 2nd order linear ODE coupled to a nonlinear static sigmoid function, the Wiener cascade (Figure12(b)). In general, Freeman [39] (and many other papers from Freeman and co-workers) defines the following neuronal population levels (Figure 12(c)).(1) K0: a subset of noninteracting neurons, either all excitatory (K0e) or inhibitory (K0i). The smallest unit in Figure 12(c) signifies a K0 subset. In each K0 subset j, the wave variable vj is governed by a 2nd order ODE: abF(vj)=v¨j+(a+b)v˙j+abvj, with a and b constants. For the output there is a sigmoid function to convert the membrane potential to pulse rate: pj=A-Bexp(C-Cevj), with A, B, and C constants.(2) KI: coupled K0 networks create the KI level. In Freeman’s terminology, twoK0e modules give a KIe module and a pair of K0i gives one KIi (Figure 12(c)). For each population j, the coupling is governed by F(vj), a function representing the incoming pulses originating from the connected populations. As shown above, the dynamics of the system’s basic K0 population is described by a 2nd order differential equation with a driving term (i.e., F(vj)). If we Laplace transform the ODE while assuming that the driving force F(vj) is a unit impulse δ(t), we get ab=[s2+(a+b)s+ab]H(s), with H(s) - transfer function. The inverse transform of the transfer function, the UIR is the dual exponential (also used in other models such as Lopes da Silva et al. [35]).(3) KII: this type of population arises when aKIe and KIi are coupled (Figure 12(c)). The distance between KII populations can be large enough that conduction delays cannot be neglected. These delays are indicated by L in Figure 12(c).(4) KIII: when KII sets are interconnected, a KIII set is created. Note that there is just one such a set in Figure12(c). In Freeman [39] there are three coupled KII sets: the olfactory bulb (OB), the anterior olfactory nucleus (AON), and the prepyriform cortex (PC). The connections in this set are longer range so that lags (latencies) were incorporated. These latencies are included in his expressions of the driving force by a backward counter J.  For example, the effect of a subpopulation from AON to one in OB is included as
(53)(19)∑J=1422PE1(J).
This term represents the weighted sum of the AON population’s pulse values PE1 in the past 14–22 steps. Note that the coefficient is divided by 9, the number of backward steps that are included in the sum.Simulations of Freeman’s model were implemented in Matlab and typical results are depicted in the bottom panel in Figure12(c). To show that the system is initially quiet, we deliver a first stimulus a little after onset of the simulation and to demonstrate that the response of the active intact network is slightly different, we deliver another stimulus at about the halfway point of the simulation (arrows in Figure 12(c)).


### 4.4. Model for the Electroencephalogram

The electrical signal that can be recorded from the scalp, the electroencephalogram (EEG), reflects the currents generated by the brain’s nerve cells. If we ignore electrical signals from other more remote sources (e.g., muscle, heart) and nonbiological artifacts, we can assume that the EEG is the weighted sum of the underlying neural network in which the weights are determined by the geometry. The EEG signal is used both in research and clinical settings. Because a single EEG signal includes the activity of millions of nerve cells, the relationship between smaller networks of the brain and the EEG signal is not necessarily a trivial one. Therefore this relationship has been the target of many modeling studies; the framework developed by physicist Paul Nunez is an example. Because it is generally thought that the slower synaptic potentials, especially those of the neocortical pyramidal cells, contribute most to the EEG signals on the surface of cortex and scalp, Nunez’ model focused on describing synaptic activity across cortex (e.g., [31, 32]). The spatiotemporal model of cortical activity is based on the schematic of neocortex shown in Figure 13. This diagram shows the effect of action potential firing activity g produced by the left cube on the synaptic activity h of the right cube. The synaptic activity has subscripts e and i for the excitatory and inhibitory synapses, respectively. Both activity variables g and h are a function of location: r→1 (left cube), r→ (right cube), and time t. The remaining variables in Figure 13 are v for conduction velocity and u(r→,t) for external, subcortical input. Consequently the delay for action potentials arriving at the right cube that were generated in the left cube is |r→-r→1|/v; therefore, the action potential firing function from the left cube arriving at the right cube is indicated as g(r→1,t-|r→-r→1|/v). Furthermore, the connectivity within the cortex is governed by functions RE and RI for the excitatory and inhibitory connections, respectively. These functions also depend on the locations r→, r→1, and fiber system (with conduction velocity v), that is, RE=RE(r→,r→1,v). The equation that describes the excitatory synaptic activity as a function of the input to the right cube can now be formulated as
(54)hE(r→,t)=u(r→,t)+∫0∞dv∫SRE(r→,r→1,v)g(r→1,t-|r→-r→1|v)d2r1.
For the inhibitory synaptic activity a simpler expression can be used since there are only local inhibitory fibers all with similar conduction velocity, that is, RI=RI(r→,r→1). Further, because of the local character of inhibition, one can neglect the delays for the inhibitory fibers: that is, now the action potential firing function is simply g(r→1,t). For the inhibitory expression, we get
(55)hI(r→,t)=u0+∫SRI(r→,r→1)g(r→1,t)d2r1.Figure 13
The model presented by Nunez is based on interaction between cortical volume units. The outputg of one volume at position r→1 arrives at its target at position r→ with a delay determined by conduction velocity (v). Every volume unit is characterized by its synaptic activity, both for excitatory (he) and inhibitory synapses (hi). In the mammalian brain, the excitatory cortical synaptic activity is assumed to make the principal contributions to the electroencephalogram (EEG).The external inputu0 represents subcortical input which is assumed constant for a given state. Because the physiological basis for this input is not obvious, Nunez dropped this term in later versions of the model. In the following, it is not very important since the effect of u0 vanishes when linearizing the expression. Note that both in (54) and (55) the connectivity functions RE and RI and/or activity function g must contain a component that translates the action potential rate into a level of synaptic activity. In his development, Nunez assumes this component to be linear by using a simple gain/attenuation factor. The same can be said for external input functions u(r→,t) and u0. They implicitly include a conversion from action potential rate to synaptic activity; the external input in terms of the action potential rate must be corrected to generate synaptic activity and (again) if this relationship is linear, one can use a simple factor. The above expressions in (54) and (55) relate spike input to synaptic activity (pulse-to-wave conversion); of course there is also the effect of the synaptic activity on the spike output (wave-to-pulse conversion). The nonlinear relationship for both synaptic activities hE and hI on action potential rate function g is commonly modeled by a sigmoid function (e.g., [32]). Nunez linearizes this relationship about an assumed fixed state g0. If we consider a small change of the action potential rate δg around state g0 and relate this to small changes in the synaptic activities δhE and δhI, we get δg=(∂g/∂hE)g=g0δhE+(∂g/∂hI)g=g0δhI. Defining HE=δhE and HI=δhI, QE=(∂g/∂hE)g=g0, and QI=-(∂g/∂hI)g=g0, we can simplify notation and get δg=QEHE-QIHI. If we now use (54) to write an expression for δhE, we get(56)δhE(r→,t)︸HE(r→,t)=δu(r→,t)︸U(r→,t)+∫0∞dv×∫SRE(r→,r→1,v)δg(r→1,t-|r→-r→1|v)︸QEHE(r→1,t-|r→-r→1|/v)-QIHI(r→1,t-|r→-r→1|/v)d2r1,
with the changes of notation indicated by the horizontal curly brackets, we obtain (Equation (11.5) in [32]):
(57)HE(r→,t)=U(r→,t)+∫0∞dv∫SRE(r→,r→1,v)[QEHE(r→1,t-|r→-r→1|v)-QIHI(r→1,t-|r→-r→1|v)]d2r1.
Following the same procedure for (55) gives (Nunez’ Equation (11.6))
(58)HI(r→,t)=∫SRI(r→,r→1)[QEHE(r→1,t)-QIHI(r→1,t)]d2r1.
One of the strong aspects of Nunez’ model is that it acknowledges different conduction velocities, the weak part is that it linearizes the relationship between synaptic activity and action potential rate, and vice versa.


### 4.5. A Field Equation for the EEG

An excellent summary of macroscopic modeling is provided in Jirsa and Haken [45]. These authors integrate temporal and spatial aspects, the firing rate models first described by Wilson and Cowan [29, 30], and the synaptic fields described by Nunez (e.g., [31, 32]). Their approach is straightforward. First they develop expressions for the synaptic and action potential activities. Then they simplify by zooming out, determine the expression for the synaptic field, and write that expression in the form of a convolution operation. Finally, they take the Fourier transform of the convolution and transform this result back into the spatiotemporal domain as a partial differential equation for the synaptic field of the excitatory neurons. Just as Nunez (Section 4.4), Jirsa and Haken assume that this field is a representative for the EEG.


#### 4.5.1. Spatiotemporal Expression for the Synaptic Field

The excitatory (e) and inhibitory (i) populations each have pulse rate E and I  (E(x,t) and I(x,t)) similar to the approach by Wilson and Cowan [29], and a wave expression ψ(ψe(x,t) and ψi(x,t)) similar to the development by Nunez. In the above, the spatial component x is three dimensional; in the remainder of the text, the spatial variable is reduced to one dimension represented by x. In the following text, we follow the derivation of the macroscopic field equation as it is visualized in Figure 14. In order to find the expression for the synaptic activity, we determine the synaptic component (top, Figure 14), while including the spike rate generator (bottom, Figure 14). The left part of Figure 14 shows the details whereas the right part shows the zoomed out version representing a macroscopic view where some detail is lost.Figure 14
Diagrams of the model described by Jirsa and Haken [45]. They explicitly model the conversion between cellular output and the synaptic input (transformation from pulse-to-wave, upper panels) and the conversion between synaptic activity and spike rate (transformation from wave-to-pulse, bottom panels). Top-left panel shows how the excitatory synaptic field ψe (as a function of time t and place x) is generated by the action potential field E by integration over the 1-dimensional cortical “area” Γ using integration variable X. The outcome (59) depends on the distribution of connectivity fe and the conversion operation He (60). The other panels in the top row show the same picture in various degrees of zooming out. The Bottom-left panel shows the generation of the action potential field E (as a function of time t and place x) by the excitatory synaptic field ψe (61). The outcome depends on the distribution of connectivity fE and the conversion operation HE (62). The other panels in the bottom row show the same in various degrees of zooming out. Note that the distribution functionfE approaches a delta function δ(x-X). The depicted relationships only show the excitatory population; for the inhibitors, similar diagrams could be constructed.(1)Synaptic Component. The synaptic activity ψ  resulting from the incoming spike activity, for the excitatory population is:
(59)ψe(x,t)︸synapticactivityexcitatorypopulationat  locationxand  timet=∫Γ︸SurfaceAreaof  thebraindXfe(x,X)︸Distribution  Functionof  Spatial  Connectivityof  axons(output)  atXtoxHe(x,X,t)︸Output  ofConversionPulse→Wave.
For the inhibitory population a similar expression can be presented. In the equations for synaptic activity, the contribution from the pulse activity over the whole cortex is included: at location x, this contribution is the product of all converted cortical activity H multiplied by the connection distribution f for x.  Synaptic activity ψ  is then computed by integrating (using integration variable X) over cortical surface Γ (Figure 14). The conversion function He determines how incoming action potentials originating from excitatory cells are converted into excitatory synaptic activity (excitatory post synaptic potentials). This is represented by a sigmoid function S that can be linearized over most of the trajectory using the constant slope ae(60)He(x,X,t)≈aeE(X,t-|x-X|ve).
Note that the E(⋯) function describes all action potentials generated in the past that could have traveled to location x originating from location X (Figure 14, top row) in a time interval |x-X|/ve, with ve—conduction velocity. The ones that have really arrived at x depend on the connectivity between x and X; this is described by function fe(x,X).(2)Spike Rate Component. The expression for the spike rate for the excitatory population E is similar to (59)
(61)E(x,t)︸spikedensityexcitatorypopulationat  locationxand  timet=∫Γ︸SurfaceAreaof  thebraindXfE(x,X)︸Distribution  Functionof  Spatial  Connectivityof  dendrites (input)atXtoxHE(x,X,t)︸Output  ofConversionWave→Pulse.
The output conversion HE: wave-to-pulse becomes an expression that reflects the conversion from: (a) the excitatory synapses, (b) the inhibitory synapses, and (c) a term for the extracortical activity. The synaptic activities for the three wave sources are simply added and a (static) nonlinearity Se converts it to pulse rate. Like in the expression in (60), there is a temporal delay term |x-X|/v that corrects for conduction velocity. The expression for the wave-to-pulse conversion becomes
(62)HE(x,X,t)︸ConversionWave→Pulse=Se(ψe(X,t-|x-X|ve)︸Contribution  of  theExcitatory  SynapticActivities  fromlocationXat timest-delay-ψi(X,t-|x-X|vi)︸Contribution  of  theInhibitory  SynapticActivities+pe(X,t-|x-X|v)︸Contribution  of  theExternal  SynapticActivities{ψe(X,t-|x-X|ve)︸Contribution  of  theExcitatory  SynapticActivities  fromlocationXat timest-delay}).The inclusion of a conduction delay is correct when assuming active propagation of synaptic input, for passive propagation (at the speed of light) such a delay is not needed. For the remainder of the derivation by Jirsa and Haken [45] this point is irrelevant because these delays will be neglected. The conversion function HE describes the contributions from all cortical areas X to location x and the distribution function fE in (61) determines the strength of the connection.


#### 4.5.2. Zooming out

In the model, cortical connectivity is described by exponential distributions of the form(1/2σ)e-|x|/σ, functions where the area under the curve evaluates to unity: ∫-∞∞(1/2σ)e-|x|/σdx=1. The advantage is that when we zoom out to a global overview of the cortex, these distributions can be replaced by delta functions. In the model, the distributions fi, fI, and fE are represented by a delta function because they are local (see for example fE in Figure 14 bottom row). This procedure for fE is depicted in Figure 14, bottom row: fE(x,X)≈δ(X-x). Distribution fe, describing the excitatory connections over a longer range, remains in the form of the exponential expression (Figure 14, top row).(1)The Expression forψe. Now we use the above to formulate the expression for the synaptic activity because this variable is considered proportional to the EEG. In this approach, we focus on solving for the excitatory component since this is considered the major contributor to the EEG signal. First we substitute the linearized version of HE (60) into ψe (59):
(63)ψe(x,t)=∫ΓdXfe(x,X)aeE(X,t-|x-X|ve).
In this expression, we substitute (61) for E(⋯), while changing integration variable X to X1, and we get ψe(x,t)=∫ΓdXfe(x,X)ae∫ΓdX1fE(X,X1)HE(X,X1,t-|x-X|/ve). Subsequently we substitute the expression for HE from (62) and replace fE with a delta function
(64)ψe(x,t)=∫ΓdXfe(x,X)ae∫ΓdX1δ(X-X1)×Se(ψe(X1,t-|X-X1|ve-|x-X|ve)-ψi(X1,t-|X-X1|ve-|x-X|ve)+pe(X1,t-|X-X1|ve-|x-X|ve)).
Now we apply the sifting property to evaluate the 2nd integral and get
(65)ψe(x,t)=ae∫ΓdXfe(x,X)Se(ψe(X,t-|x-X|ve)-ψi(X,t-|x-X|ve)+pe(X,t-|x-X|ve)).
For the inhibitory population we follow a similar procedure and get
(66)ψi(x,t)=ai∫ΓdXfi(x,X)Si(ψe(X,t-|x-X|vi)-ψi(X,t-|x-X|vi)+pi(X,t-|x-X|vi)).
If we now evaluate the integral in this expression, while replacing the distribution function fi with the delta function δ(X-x), so we can use the sifting property, we obtain the result for the inhibitory synaptic activity: ψi(x,t)=aiSi(ψe(x,t)-ψi(x,t)+pi(x,t)). Subsequently we follow the procedure by Jirsa and Haken [45] and only take the linear aspect of the Si function, using a slope of αi and finally solve for ψi(67)ψi(x,t)=aiαi1+aiαi(ψe(x,t)+pi(x,t)).
Substitution of (67) into (65) gives(68)ψe(x,t)=ae∫ΓdXfe(x,X)Se×{ψe(X,t-|x-X|ve)-aiαi1+aiαiψe(X,t-|x-X|ve)︷ρ~ψe+pe(X,t-|x-X|ve)-aiαi1+aiαipi(X,t-|x-X|ve)︸p(X,t-|x-X|/ve)}.Using the simplifications indicated by the horizontal curly brackets above we have the expression for the excitatory synaptic activity:(69)ψe(x,t)=ae∫ΓdXfe(x,X)Se{ρ~ψe(X,t-|x-X|ve)+p(X,t-|x-X|ve)}.


#### 4.5.3. Expression forψe  as a Convolution

The next step is to write the result forψe in the form of a convolution. First, the notation for the integrand in (69) can be simplified
(70)fe(x,X)aeSe{ρ~ψe(X,t-|x-X|ve)+p(X,t-|x-X|ve)}︸ρ(X,t-|x-X|/ve)=fe(x,X)ρ(X,t-|x-X|ve).
Because the excitatory connections cover a larger area its distribution function cannot be replaced by δ. Therefore we use the exponential distribution for fe and get: (1/2σe)e-|x-X|/σeρ(X,t-|x-X|/ve). Although it seems that we make things more complicated, we rewrite this as an integral: (1/2σe)e-|x-X|/σe∫-∞∞ρ(X,T)δ(T-(t-|x-X|/ve))dT. Now we substitute this result and the expression in (70) in (69) and rearrange the order of the terms
(71)ψe(x,t)=∫Γ∫-∞∞δ(t-T-|x-X|ve)12σee-|x-X|/σe︸G(x-X,t-T)×ρ(X,T)dXdT.
The part (1/2σe)e-|x-X|/σeδ((t-|x-X|/ve)-T) can be interpreted as a unit impulse response or Green’s function G.Now we complete this step and write the expression in the form of a convolution with respect to space and time(72)ψe(x,t)=∫Γ∫-∞∞G(x-X,t-T)ρ(X,T)dXdT.


#### 4.5.4. From Convolution to the Field Equation

Now we take advantage of the convolution form of (72). With the spatio-temporal Fourier transforms for ψ(x,t), ρ(x,t) and G(x,t), we can write the transform of the convolution as a product in the frequency domain. To summarize we have the following transform pairs:
(73)ψe(x,t)⟺ψe(k,ω),ρ(x,t)⟺ρ(k,ω),G(x,t)⟺ω02+jω0ωve2k2+(ω0+jω)2,withω0=veσe.
Using ψe(x,t)=G(x,t)⊗ρ(x,t)⇔ψe(k,ω)=G(k,ω)×ρ(k,ω), we get:
(74)ψe(k,ω)=ω02+jω0ωve2k2+(ω0+jω)2ρ(k,ω).
Note that in the temporal Fourier transform jω and -ω2 represent the transformed first and second derivatives in the time domain, while the representatives for the first and second spatial derivatives are +jk and -k2. Transforming (74), we complete our last step and get the equation for the spatio-temporal excitatory synaptic activity
(75)(ω02-ve2∂2∂x2)ψe(x,t)+2ω0ψ˙e(x,t)+ψ¨e(x,t)=(ω02+ω0∂∂t)ρ(x,t).
An example of simulated activity is depicted in Figure 15. In this example, we show the model’s response to a half sine wave pulse in the center of a cortical surface. As expected, the evoked activity is damped over both space and time. The result in (75) is also similar to earlier results obtained for linearized cellular and network models. If we ignore the spatial component in (75), we get a second order ODE with a force term, similar to (26) and (44). Although in this case, a nonlinearity is included in the factor ρ of the force term.Simulation of (75) for times 5, 7, 10, and 25 showing the response of a cortical surface to a stimulus pulse in the middle of the surface. The pulse was half a period of a sine wave between times 0 and 10. Parameters used are ω0=0.15; ve=1.52; ρ~=0.9; ae=1.2; αe=1.2.
(a)
t
=
5
(b)
t
=
7
(c)
t
=
10
(d)
t
=
25


## 4.1. The Ising Spin Model

Two- or three-dimensional models of the Ising-spin lattice, where the spins represent the neurons, are not easy to solve. A mean field approach is the simplest approximation to describe such a system. In the lattice, each node experiences the sum of the magnetic forces created by the other nodes and, if present, an external magnetic field. In the mean field approach, we replace the internal field generated by the spins by its average value. It is beyond the scope of this paper to derive the mean field relationship, which can be found in many physics textbooks (e.g., [67, Chapter 7]). Using the mean field expression, we can describe the equilibrium state of the lattice of spins—or network of model neurons—by the expression
(32)(Jnm+h)-kTtanh-1(m)=0,
with k—Boltzmann constant (in the following k is set to unity), T—temperature, J—magnetic coupling strength, n—number of neighbors, m—magnetization of the lattice, and h—an externally applied magnetic field. As shown in the following, (32) shows sudden transitions in behavior according to the so-called cusp catastrophe [68]. First we can determine the fixed points for the expression graphically by finding the intersections of both terms, Ttanh-1(m) and Jnm+h as a function of m. When depicting these functions for different T and h, it can be seen that they intersect at either 3 or at 1 point(s). In case of the intersection at 3 points, the middle one is an unstable fixed point, the other two are stable ones; therefore the system is bistable in one area of the T-h plane and monostable elsewhere.The transition between bi- and mono-stability occurs at the bifurcation where the straight lineJnm+h intersects the sigmoid Ttanh-1(m) tangentially; this indicates that two conditions must be satisfied at the transition in point Tc, h*:
(33)Tctanh-1(m)=Jnm+h*,
and since they must be tangent
(34)ddm[Tctanh-1(m)]=ddm[Jnm+h*]→Tc1-m2=Jn.
From these two expressions, one can find the expression for Tc and h* (as a function of m). The plot of these functions in the T-h plane is the so-called stability diagram (Figure 8(a)). These lines mark the border where the system transitions from mono- (1 fixed point) to bi-stability (3 fixed points with one of the points unstable). Another common presentation is a bifurcation diagram; here we have a codimension-2 case with two independent variables T and h. Note that for h=0, m=0 is always a solution (fixed point) for (32). If we use Jn=1, we have a stable fixed point at m=0 for T>1 and an unstable one plus two stable fixed points for T<1. This is depicted in Figure 8(b); we can see that we have the characteristic of a supercritical pitchfork bifurcation. This plot can be considered as cross-sections of a 3D plot of m(h,T) shown in Figure 8(c); this surface depicts the so-called cusp catastrophe.Properties of the mean field equation of the Ising model. (a) The stability diagram. (b) Bifurcation diagram. (c) The cusp catastrophe shows the dependence of magnetizationm on both the temperature (T) and an external magnetic field (h).
(a)
Stability diagram
(b)
Bifurcation diagram ath=0
(c)
Cusp catastrophe


## 4.2. The Wilson and Cowan Model

The 1st model by Wilson and Cowan [29] is exclusively temporal and describes the behavior of two neuronal populations—one excitatory and one inhibitory—embedded in neural tissue. It is more realistic than using the Ising model that only considers the equivalent of excitation. In the Wilson-Cowan model, population properties are derived from assumed statistics of the neurons [29]. Each population receives input from itself, the other population, and external sources. The input to a population affects the neurons of that population that are available (the inactive proportion of the population). The model is nonlinear because the incoming activity for the population is related to its output through a nonlinear sigmoid function S. In a 2nd version of the model that will not be further discussed here, Wilson and Cowan [30] extended their model to include the spatial component as well. Instead of formulating the model from underlying statistical properties, a coarse grained version of the model can be directly derived from the diagram in Figure 9. Two populations, the excitatory and inhibitory, are coupled with strengths c1–c4. The external population affects the excitatory and inhibitory networks with activities P and Q. We denote excitatory activity by E and inhibitory activity by I, the excitatory population in Figure 9 receives input that can be quantified as c1E-c2I+P. This input of the population is converted to an output expression via a sigmoid function Se{⋯}.  In most followers of their approach, such a sigmoid curve is employed to represent conversion from membrane potential to spike rate and/or vice versa. In the Wilson and Cowan’s approach, the sigmoid curve represents the cumulative distribution of the probability that a neuron will fire at a given input, that is, the distance between membrane potential and firing threshold. Further, assuming that the unavailable part of the population, that is, the active + refractory portion is proportional (by a fraction re) to the activity level E then the receptive population (available portion of E) is 1-reE.Figure 9
Diagram of the Wilson-Cowan Model containing excitatory and inhibitory populations. The coupling strength between the populations is denoted with constantsc1–c4 and external activity is symbolized by P and Q. Note that the model contains nonlinear functions to relate input and output for each population (not shown in the diagram).In Wilson and Cowan [29], the authors use a value a little bit less than unity for ke, reflecting that the maximum values for the sigmoid function Se will be a little bit less than one (see pages 9-10 in the original paper). Therefore, the available population becomes ke-reE. The formalism for population activity also must describe spontaneous decay proportional to E, reflecting that a nonstimulated group of neurons converge to a base level of activity (zero, reflecting a low level of spontaneous activity in this model). Finally all changes dE/dt are weighted by the time constant of the neuronal population τe. Putting this all together the coarse grained equation for excitatory activity becomes
(35)τe︸EXCITATORYTIME CONSTANTdEdt=-E︸DECAY+(ke-reE)︸AVAILABLEPOPULATIONSe{c1E-c2I+P}︸.INPUT︷SIGMOID
A similar approach for the inhibitory population yields
(36)τidIdt=-I+(ki-riI)Si{c3E-c4I+Q}.
In their 1972 paper [29], Wilson and Cowan derive these relationships from first principles followed by a coarse graining procedure (their equations (11) and (12)). When looking into the two dimensional nonlinear Wilson-Cowan model, there should be a variety of behaviors including a steady activity scenario, oscillatory activity with damping, and limit cycles (Figure 10).Phase plane representations of the Wilson-Cowan model during equilibrium, oscillatory, and limit-cycle activities. These plots correspond to the parameters in Figures 4, 10, and 11 of the original 1972 paper [29]. Note that in the equilibrium case there are two stable fixed points: one in a “down” and one in an “up” state. Also note that in these examples, oscillations are obtained by flattening the E null-cline (magenta) and making the I null-cline (orange) steeper. Sample trajectories are depicted in blue. (a) Three fixed points, two of them are stable and one (in between) is unstable. The eigenvalues of these points from left to right are as follows: −0.59, −1.13; 0.73, −1.6; −1.43, −2.88. (b) Two fixed points, one is clearly associated with a damped oscillation. The eigenvalues from left to right are as follows: −0.34, −0.98; -0.09±0.86j. (c) The limit cycle with a fixed point inside. The eigenvalues of the fixed point are 0.12±1.88j.
(a)
(b)
(c)


### 4.2.1. Linearization of the Wilson-Cowan Equations

To study resonance phenomena, we will consider the effect onE and I by small sinusoidal perturbations (e.g., evoked by electrical stimulation) close to an equilibrium state. When stimulating electrically in the extracellular space and assuming similar spatial distribution of excitatory and inhibitory cells, it is reasonable to assume that both the excitatory and inhibitory populations receive the same stimulus current. Thus the external inputs P and Q are considered equal. Further simplifying the sigmoid curves Se and Si to lines with slopes α and β, the linearized versions of (35) and (36) become
(37)E˙=-Eτe+(ke-reE)τeα(c1E-c2I+P),I˙=-Iτi+(ki-riI)τiβ(c3E-c4I+P).
Rewriting these results and using ατ=α/τe and βτ=β/τi, we get
(38)E˙=f(E,I,P)=-Eτe+keατc1E-keατc2I+keατP-reατc1E2+reατc2EI-reατEP,I˙=g(E,I,P)=-Iτi+kiβτc3E-kiβτc4I+kiβτP-riβτc3EI+riβτc4I2-riβτIP.
We now further linearize and consider small perturbations δE, δI, δP around the equilibrium indicated by *(39)δE˙=(∂f∂E)*δE+(∂f∂I)*δI+(∂f∂P)*δP,δI˙=(∂g∂E)*δE+(∂g∂I)*δI+(∂g∂P)*δP.
The partial derivatives determined at equilibrium * are:
(40)εE=(∂f∂E)*=(-1τe+keατc1-2reατc1E+reατc2I-reατP)*,εI=(∂f∂I)*=(-keατc2+reατc2E)*,εP=(∂f∂P)*=(keατ-reατE)*,γE=(∂g∂E)*=(kiβτc3-riβτc3I)*,γI=(∂g∂I)*=(-1τi-kiβτc4-riβτc3E+2riβτc4I-riβτP)*,γP=(∂g∂P)*=(kiβτ-riβτI)*.
Combining these expressions and substituting e and i for δE and δI, and making the perturbation δP sinusoidal psinωt (where amplitude p is a small value), we can simplify notation
(41)e˙=εEe+εIi+εPpsinωt,(42)i˙=γEe+γIi+γPpsinωt.
These two 1st order differential equations can also be written as a single 2nd order one. Taking the derivative of (42) and substituting the expression in (41) for i gives
(43)e¨=εEe˙+εI(γEe+γIi+γPpsinωt)+εPωpcosωt.
Using (41), we can substitute i=(e˙-εEe-εPpsinωt)/εI, and we get the standard expression for a 2nd order differential equation with a driving force
(44)e¨-(εE+γI)︸be˙+(εEγI-εIγE)︸ce=εPpdsinωtdt+(εIγP-εPγI)psinωt︸driving"force".
This equation is again of the same form as the equation that governs a parallel RCL circuit as shown in Figure 3. The 2nd-order ODE governing this RCL circuit is (26). This equation is also very similar to the network equation by Jirsa and Haken [45] in which the spatial derivative is set to zero (Section 4.5). Thus the linearized network model is capable of resonance behavior, described in Section 2.1.4. Since we deal with a linearized version of a nonlinear system, unstable solutions of the linearized case may relate to bifurcations in the full nonlinear version of the model (see Section 6.1 for further discussion).


## 4.2.1. Linearization of the Wilson-Cowan Equations

To study resonance phenomena, we will consider the effect onE and I by small sinusoidal perturbations (e.g., evoked by electrical stimulation) close to an equilibrium state. When stimulating electrically in the extracellular space and assuming similar spatial distribution of excitatory and inhibitory cells, it is reasonable to assume that both the excitatory and inhibitory populations receive the same stimulus current. Thus the external inputs P and Q are considered equal. Further simplifying the sigmoid curves Se and Si to lines with slopes α and β, the linearized versions of (35) and (36) become
(37)E˙=-Eτe+(ke-reE)τeα(c1E-c2I+P),I˙=-Iτi+(ki-riI)τiβ(c3E-c4I+P).
Rewriting these results and using ατ=α/τe and βτ=β/τi, we get
(38)E˙=f(E,I,P)=-Eτe+keατc1E-keατc2I+keατP-reατc1E2+reατc2EI-reατEP,I˙=g(E,I,P)=-Iτi+kiβτc3E-kiβτc4I+kiβτP-riβτc3EI+riβτc4I2-riβτIP.
We now further linearize and consider small perturbations δE, δI, δP around the equilibrium indicated by *(39)δE˙=(∂f∂E)*δE+(∂f∂I)*δI+(∂f∂P)*δP,δI˙=(∂g∂E)*δE+(∂g∂I)*δI+(∂g∂P)*δP.
The partial derivatives determined at equilibrium * are:
(40)εE=(∂f∂E)*=(-1τe+keατc1-2reατc1E+reατc2I-reατP)*,εI=(∂f∂I)*=(-keατc2+reατc2E)*,εP=(∂f∂P)*=(keατ-reατE)*,γE=(∂g∂E)*=(kiβτc3-riβτc3I)*,γI=(∂g∂I)*=(-1τi-kiβτc4-riβτc3E+2riβτc4I-riβτP)*,γP=(∂g∂P)*=(kiβτ-riβτI)*.
Combining these expressions and substituting e and i for δE and δI, and making the perturbation δP sinusoidal psinωt (where amplitude p is a small value), we can simplify notation
(41)e˙=εEe+εIi+εPpsinωt,(42)i˙=γEe+γIi+γPpsinωt.
These two 1st order differential equations can also be written as a single 2nd order one. Taking the derivative of (42) and substituting the expression in (41) for i gives
(43)e¨=εEe˙+εI(γEe+γIi+γPpsinωt)+εPωpcosωt.
Using (41), we can substitute i=(e˙-εEe-εPpsinωt)/εI, and we get the standard expression for a 2nd order differential equation with a driving force
(44)e¨-(εE+γI)︸be˙+(εEγI-εIγE)︸ce=εPpdsinωtdt+(εIγP-εPγI)psinωt︸driving"force".
This equation is again of the same form as the equation that governs a parallel RCL circuit as shown in Figure 3. The 2nd-order ODE governing this RCL circuit is (26). This equation is also very similar to the network equation by Jirsa and Haken [45] in which the spatial derivative is set to zero (Section 4.5). Thus the linearized network model is capable of resonance behavior, described in Section 2.1.4. Since we deal with a linearized version of a nonlinear system, unstable solutions of the linearized case may relate to bifurcations in the full nonlinear version of the model (see Section 6.1 for further discussion).


## 4.3. Models Inspired by the Wilson-Cowan Approach

The model described by Lopes da Silva et al. [35] is a model of thalamus and one of the first that  is inspired by the Wilson-Cowan approach. A similar treatment of the mesoscopic model, with altered synaptic function, was given by van Rotterdam et al. [36]. The base model describes a thalamic network; however, the populations are identical to the ones in the cortical Wilson-Cowan model: one excitatory and one inhibitory cell group. Similar to the Wilson-Cowan approach they assume an invertible function f that translates membrane potential, Ve(t) and Vi(t)  (i.e., the input) for the excitatory and inhibitory populations, into a rate of action potentials for each population E(t) and I(t) (i.e., the output). For the excitatory population we then get the following relationships: E(t)=fe[Ve(t)]  and Ve(t)=fe-1[E(t)]. Here fe and fe-1 are the function and its inverse that translates membrane potential into spike rate and vice versa. Both the functions fe and its equivalent for the inhibitory population fi are basically the static nonlinearities in their model. Each nonlinear component is subsequently approached by a Taylor series about the mean value of E denoted as E-: fe-1[E(t)]=ae0E-+ae1(E(t)-E-)+ae2(E(t)-E-)2+⋯. Subsequently Lopes da Silva et al. simplify notation by only considering the deviation from the mean only, e(t)=E(t)-E-. This is done for all variables in the model, and next they linearize about the mean by neglecting the 2nd and higher order terms in the Taylor series so that the original equations can be approximated by ve(t)=fe-1[e(t)]=ae1e(t) and vi(t)=fi-1[i(t)]=ai1i(t). External input P(t) or its demeaned version p(t) perturbs the system (top Left in Figure 11(a)). For the excitatory component in the diagram in the Laplace domain we get Ve(s)=He(s)P(s), with the following Laplace transform pairs: Ve(s)⇔ve(t),  He(s)⇔he(t), and P(s)⇔p(t). Here he and He are the unit impulse response and transfer function of the excitatory network.(a) Diagram of the Model (after [35]). P(t) is the external input; he and hi are the population’s linear unit impulse responses reflecting the synaptic function. Variables Ve and E are the population’s mean membrane potential and firing rate respectively; Vi and I are the same for the inhibitory population. Functions fe and fi are the nonlinear functions that convert population membrane potential into population firing rate. C1 and C2 are coupling constants. (b) Result for the spectrum following (49).
(a)
(b)In a similar fashion we include the effect of the inhibition in the diagram in Figure11(a) and get
(45)Ve(s)=ae1E(s)=He(s)P(s)-C2I(s)Hi(s).
Following the same procedure for the linearized inhibitory component (bottom part in Figure 11(a)) gives
(46)Vi(s)=ai1I(s)=C1E(s)He(s).
If we substitute the expression for I(s) (obtained from (46)) into (45) and use E(s)=Ve(s)/ae1 (from (45)), we obtain the following expression for Ve(s):
(47)Ve(s)=He(s)P(s)1+C1C2Hi(s)He(s)/ai1ae1.
Then the authors use a linear system to model synaptic input with unit impulse responses (UIRs) he and hi. In Lopes da Silva et al. [35], the synaptic UIR is a dual exponential A[e-a1t-e-a2t] and in van Rotterdam et al. [36] it is an alpha-function Ate-at; in both cases, these functions are only considered for t≥0.The associated Laplace and Fourier transforms are, of course, slightly different. Here we follow the paper of Lopes da Silva et al. [35] and use the dual exponential unit impulse response for both the excitatory and inhibitory populations. These functions and their Laplace transforms are
(48)he(t)=A[e-a1t-e-a2t]⟺He(s)=A(a2-a1)(s+a1)(s+a2),hi(t)=B[e-b1t-e-b2t]⟺Hi(s)=B(b2-b1)(s+b1)(s+b2).
We can substitute this into (47) and simplify
(49)Ve(s)=A(a2-a1)(s+b1)(s+b2)(s+a1)(s+a2)(s+b1)(s+b2)+KP(s),
with K=(C1C2/ai1ae1)(a2-a1)(b2-b1)AB.We can substitutes by jω in Ve(s) to obtain an expression for the spectrum of ve. The spectrum can be further quantified by assuming that the thalamic input p(t) is white noise and therefore P(ω) is constant. The numerical values we use for Figure 11(b) are the same as the ones in the example given in Lopes da Silva et al. [35]: A=1.65 mV, B=32 mV, C1=32 cells, C2=3 cells, a1=55 s−1, a2=605 s−1, b1=27.5 s−1, b2=55 s−1, qe1=1/ae1, qi1=1/ai1, and qi1qe1=4.55×106 (K≈3.5e8). It can be seen that the spectrum peaks around 10 Hz, the alpha rhythm frequently encountered in the EEG of subjects with eyes closed.


### 4.3.1. Neural Mass Models

The Wilson-Cowan Equations [29] and Lopes da Silva et al. approach [35] inspired many other models. Instead of considering individual neurons, these models determine the activity of coupled populations of nerve cells, each population representing a neural mass. A very basic model, using two populations representing a coupled excitatory and inhibitory network is described by Dayan and Abbott [69, pages 265–270]. The rate of change for each network is inversely proportional to its activity level, and proportional to the input it receives from itself, the other population, and an external input. This results in a set of two coupled differential equations. Although the model seems to be linear, their model is actually nonlinear because the population activities are limited to positive values only. Therefore their model is capable of generating limit cycles.Jansen’s neural mass model [37] is applied by Grimbert and Faugeras [38] to create a cortical network unit (Figure 12(a)). It is a temporal model of neuronal populations inspired by the Lopes da Silva et al. [35] approach where each population is represented by a Wiener cascade ([70], Chapters 2–5): a dynamic linear component, modeling the synaptic conversion process, and a nonlinear static component that represents the conversion of membrane potential to spike rate (Figure 12(b)). As it can be seen in Figure 12(b), input x1 (input spike rate) is first convolved with the unit impulse response (UIR) h of the linear dynamic component, the result x2=x1⊗h (population membrane potential; ⊗-convolution) is subsequently converted with a static sigmoid to produce output x3=σ(x2) (output spike rate). The structure of the model by Grimbert and Faugeras [38] is similar to the canonical cortical model of Douglas and Martin [71] in the sense that it also includes two excitatory populations E, one inhibitory population I, and an external input P. The development of the equations that govern the activities is straightforward. First we use the UIR to write a convolution expression. Again, the synapse is considered a linear system and in this model, the synaptic impulse response function is defined as an alpha function (as in [36]):
(50)h(t)={αβte-βt,t≥0,0,otherwise,
in which α and β are population-specific constants. The Laplace transform of the unit impulse response is: H(s)=αβ/(s+β)2. In the time/Laplace domain we may link the input x(t)/transformed input X(s) to the output y(t)/transformed output Y(s) of any linear system by the transform pair: Y(s)=H(s)X(s)⇔y(t)=h(t)⊗x(t). Second, we write the convolution relationship as an ODE. Substituting the above expression for H in the Laplace transformed expression and we get (s+β)2Y(s)=αβX(s). This result can be re written as s2Y(s)+2βsY(s)+β2Y(s)=αβX(s). Recall that sY(s) and s2Y(s) denote the Laplace transforms of the 1st and 2nd derivatives of y(t) in the time domain; using this, we can transform the above expression into a 2nd order ODE for the linear synaptic process: y¨(t)=αβx(t)-2βy˙(t)-β2y(t), which can be written as two 1st order ODEs
(51)y˙(t)=z(t),z˙(t)=αβx(t)-2βz(t)-β2y(t).
Now we add the static nonlinearity employed for the membrane potential to pulse rate conversion, a sigmoid curve σ(v) (Figure 12(b)). In this study, it is defined as a Boltzmann function with maximum rate vmax, slope r and a 50% level threshold v0: σ(v)=vmax/(1+e-r(v-v0)). This completes the formalism for one population, a single Wiener cascade in Figure 12(b). The pair of ODEs plus the static nonlinearity can be defined for each of the three populations in the neural mass model (two excitatory and one inhibitory, Figure 12(a)). This means we will have a pair of 1st order differential equations per population, generating a total of six equations (Equation (3) in [38])
(52)y˙0(t)=y3(t),y˙3(t)=Aaσ(y1(t)-y2(t))-2ay3(t)-a2y0(t),y˙1(t)=y4(t),y˙4(t)=Aa{p(t)+C2σ(C1y0)}-2ay4(t)-a2y1(t),y˙2(t)=y5(t),y˙5(t)=BbC4σ(C3y0)-2by5(t)-b2y2(t).
The parameters y0–y2 denote the membrane potentials, and A, a, B, b are the constants in the population specific alpha functions. Note that in these equations, some outputs y and sigmoid functions are weighted by connectivity strengths C1–C4. The focus of the paper is on the variable y, the aggregate membrane potential of the main population, which is thought to be proportional with the local field potential. A simulation using the following parameters, A=3.25, a=100, B=22, b=50, C1=135, C2=0.8*C1, C3=0.25*C1, C4=0.25*C1, vmax=5, r=0.56, and v0=6, shows a signal that is similar to the EEG alpha rhythm (Figure 12(a)).Neural mass models. (a) A diagram of the model of a cortical unit employed by Grimbert and Faugeras [38] and a sample trace showing an oscillation similar to an EEG alpha rhythm. (b) Representation of a Wiener cascade used to model a population of nerve cells. The Wiener cascade consists of a linear dynamical module and a static nonlinear one. (c) Schematic representation of the model described by Freeman [39] and sample traces showing chaotic behavior of the impulse responses (impulse input—arrows).
(a)
(b)
(c)


### 4.3.2. Modeling the Olfactory System

The models created by Freeman in the 1980s used a combined network modeling and experimental approach to study the electrical activity in rabbit olfactory system. Although Freeman’s approach may not have been directly inspired by the Wilson-Cowan equations and/or the models by Lopes da Silva and co-workers, the fundamental building blocks of his model are neural masses. The appeal of Freeman’s approach is that it is strongly based on the anatomy of the olfactory system, and that he relates modeling results to electrophysiological experiments. As compared to the previous models, a unique addition is that it contains latencies L1–L4, modeling the conduction delays occurring in between populations [39]. Differential equations without delays need initial values, whereas equations with delays require initial functions to determine the delayed response of the system. Equations with delays can exhibit behavior that cannot be displayed by equations without delay of the same order. It is given that delays in the nervous system occur because of the finite conduction velocity of nerve fiber activity. In addition to Freeman’s model, several others have explored the role of these delays or included the delays in the formalism (e.g., [65, 66, 72–74]).In his models, Freeman describes different population levels. Within the most basic population type, which is defined as K0, he includes pulse-to-wave conversion, linear spatiotemporal integration, wave-to-pulse conversion, and conduction of the output. This basic population is modeled with a 2nd order linear ODE coupled to a nonlinear static sigmoid function, the Wiener cascade (Figure12(b)). In general, Freeman [39] (and many other papers from Freeman and co-workers) defines the following neuronal population levels (Figure 12(c)).(1) K0: a subset of noninteracting neurons, either all excitatory (K0e) or inhibitory (K0i). The smallest unit in Figure 12(c) signifies a K0 subset. In each K0 subset j, the wave variable vj is governed by a 2nd order ODE: abF(vj)=v¨j+(a+b)v˙j+abvj, with a and b constants. For the output there is a sigmoid function to convert the membrane potential to pulse rate: pj=A-Bexp(C-Cevj), with A, B, and C constants.(2) KI: coupled K0 networks create the KI level. In Freeman’s terminology, twoK0e modules give a KIe module and a pair of K0i gives one KIi (Figure 12(c)). For each population j, the coupling is governed by F(vj), a function representing the incoming pulses originating from the connected populations. As shown above, the dynamics of the system’s basic K0 population is described by a 2nd order differential equation with a driving term (i.e., F(vj)). If we Laplace transform the ODE while assuming that the driving force F(vj) is a unit impulse δ(t), we get ab=[s2+(a+b)s+ab]H(s), with H(s) - transfer function. The inverse transform of the transfer function, the UIR is the dual exponential (also used in other models such as Lopes da Silva et al. [35]).(3) KII: this type of population arises when aKIe and KIi are coupled (Figure 12(c)). The distance between KII populations can be large enough that conduction delays cannot be neglected. These delays are indicated by L in Figure 12(c).(4) KIII: when KII sets are interconnected, a KIII set is created. Note that there is just one such a set in Figure12(c). In Freeman [39] there are three coupled KII sets: the olfactory bulb (OB), the anterior olfactory nucleus (AON), and the prepyriform cortex (PC). The connections in this set are longer range so that lags (latencies) were incorporated. These latencies are included in his expressions of the driving force by a backward counter J.  For example, the effect of a subpopulation from AON to one in OB is included as
(53)(19)∑J=1422PE1(J).
This term represents the weighted sum of the AON population’s pulse values PE1 in the past 14–22 steps. Note that the coefficient is divided by 9, the number of backward steps that are included in the sum.Simulations of Freeman’s model were implemented in Matlab and typical results are depicted in the bottom panel in Figure12(c). To show that the system is initially quiet, we deliver a first stimulus a little after onset of the simulation and to demonstrate that the response of the active intact network is slightly different, we deliver another stimulus at about the halfway point of the simulation (arrows in Figure 12(c)).


## 4.3.1. Neural Mass Models

The Wilson-Cowan Equations [29] and Lopes da Silva et al. approach [35] inspired many other models. Instead of considering individual neurons, these models determine the activity of coupled populations of nerve cells, each population representing a neural mass. A very basic model, using two populations representing a coupled excitatory and inhibitory network is described by Dayan and Abbott [69, pages 265–270]. The rate of change for each network is inversely proportional to its activity level, and proportional to the input it receives from itself, the other population, and an external input. This results in a set of two coupled differential equations. Although the model seems to be linear, their model is actually nonlinear because the population activities are limited to positive values only. Therefore their model is capable of generating limit cycles.Jansen’s neural mass model [37] is applied by Grimbert and Faugeras [38] to create a cortical network unit (Figure 12(a)). It is a temporal model of neuronal populations inspired by the Lopes da Silva et al. [35] approach where each population is represented by a Wiener cascade ([70], Chapters 2–5): a dynamic linear component, modeling the synaptic conversion process, and a nonlinear static component that represents the conversion of membrane potential to spike rate (Figure 12(b)). As it can be seen in Figure 12(b), input x1 (input spike rate) is first convolved with the unit impulse response (UIR) h of the linear dynamic component, the result x2=x1⊗h (population membrane potential; ⊗-convolution) is subsequently converted with a static sigmoid to produce output x3=σ(x2) (output spike rate). The structure of the model by Grimbert and Faugeras [38] is similar to the canonical cortical model of Douglas and Martin [71] in the sense that it also includes two excitatory populations E, one inhibitory population I, and an external input P. The development of the equations that govern the activities is straightforward. First we use the UIR to write a convolution expression. Again, the synapse is considered a linear system and in this model, the synaptic impulse response function is defined as an alpha function (as in [36]):
(50)h(t)={αβte-βt,t≥0,0,otherwise,
in which α and β are population-specific constants. The Laplace transform of the unit impulse response is: H(s)=αβ/(s+β)2. In the time/Laplace domain we may link the input x(t)/transformed input X(s) to the output y(t)/transformed output Y(s) of any linear system by the transform pair: Y(s)=H(s)X(s)⇔y(t)=h(t)⊗x(t). Second, we write the convolution relationship as an ODE. Substituting the above expression for H in the Laplace transformed expression and we get (s+β)2Y(s)=αβX(s). This result can be re written as s2Y(s)+2βsY(s)+β2Y(s)=αβX(s). Recall that sY(s) and s2Y(s) denote the Laplace transforms of the 1st and 2nd derivatives of y(t) in the time domain; using this, we can transform the above expression into a 2nd order ODE for the linear synaptic process: y¨(t)=αβx(t)-2βy˙(t)-β2y(t), which can be written as two 1st order ODEs
(51)y˙(t)=z(t),z˙(t)=αβx(t)-2βz(t)-β2y(t).
Now we add the static nonlinearity employed for the membrane potential to pulse rate conversion, a sigmoid curve σ(v) (Figure 12(b)). In this study, it is defined as a Boltzmann function with maximum rate vmax, slope r and a 50% level threshold v0: σ(v)=vmax/(1+e-r(v-v0)). This completes the formalism for one population, a single Wiener cascade in Figure 12(b). The pair of ODEs plus the static nonlinearity can be defined for each of the three populations in the neural mass model (two excitatory and one inhibitory, Figure 12(a)). This means we will have a pair of 1st order differential equations per population, generating a total of six equations (Equation (3) in [38])
(52)y˙0(t)=y3(t),y˙3(t)=Aaσ(y1(t)-y2(t))-2ay3(t)-a2y0(t),y˙1(t)=y4(t),y˙4(t)=Aa{p(t)+C2σ(C1y0)}-2ay4(t)-a2y1(t),y˙2(t)=y5(t),y˙5(t)=BbC4σ(C3y0)-2by5(t)-b2y2(t).
The parameters y0–y2 denote the membrane potentials, and A, a, B, b are the constants in the population specific alpha functions. Note that in these equations, some outputs y and sigmoid functions are weighted by connectivity strengths C1–C4. The focus of the paper is on the variable y, the aggregate membrane potential of the main population, which is thought to be proportional with the local field potential. A simulation using the following parameters, A=3.25, a=100, B=22, b=50, C1=135, C2=0.8*C1, C3=0.25*C1, C4=0.25*C1, vmax=5, r=0.56, and v0=6, shows a signal that is similar to the EEG alpha rhythm (Figure 12(a)).Neural mass models. (a) A diagram of the model of a cortical unit employed by Grimbert and Faugeras [38] and a sample trace showing an oscillation similar to an EEG alpha rhythm. (b) Representation of a Wiener cascade used to model a population of nerve cells. The Wiener cascade consists of a linear dynamical module and a static nonlinear one. (c) Schematic representation of the model described by Freeman [39] and sample traces showing chaotic behavior of the impulse responses (impulse input—arrows).
(a)
(b)
(c)


## 4.3.2. Modeling the Olfactory System

The models created by Freeman in the 1980s used a combined network modeling and experimental approach to study the electrical activity in rabbit olfactory system. Although Freeman’s approach may not have been directly inspired by the Wilson-Cowan equations and/or the models by Lopes da Silva and co-workers, the fundamental building blocks of his model are neural masses. The appeal of Freeman’s approach is that it is strongly based on the anatomy of the olfactory system, and that he relates modeling results to electrophysiological experiments. As compared to the previous models, a unique addition is that it contains latencies L1–L4, modeling the conduction delays occurring in between populations [39]. Differential equations without delays need initial values, whereas equations with delays require initial functions to determine the delayed response of the system. Equations with delays can exhibit behavior that cannot be displayed by equations without delay of the same order. It is given that delays in the nervous system occur because of the finite conduction velocity of nerve fiber activity. In addition to Freeman’s model, several others have explored the role of these delays or included the delays in the formalism (e.g., [65, 66, 72–74]).In his models, Freeman describes different population levels. Within the most basic population type, which is defined as K0, he includes pulse-to-wave conversion, linear spatiotemporal integration, wave-to-pulse conversion, and conduction of the output. This basic population is modeled with a 2nd order linear ODE coupled to a nonlinear static sigmoid function, the Wiener cascade (Figure12(b)). In general, Freeman [39] (and many other papers from Freeman and co-workers) defines the following neuronal population levels (Figure 12(c)).(1) K0: a subset of noninteracting neurons, either all excitatory (K0e) or inhibitory (K0i). The smallest unit in Figure 12(c) signifies a K0 subset. In each K0 subset j, the wave variable vj is governed by a 2nd order ODE: abF(vj)=v¨j+(a+b)v˙j+abvj, with a and b constants. For the output there is a sigmoid function to convert the membrane potential to pulse rate: pj=A-Bexp(C-Cevj), with A, B, and C constants.(2) KI: coupled K0 networks create the KI level. In Freeman’s terminology, twoK0e modules give a KIe module and a pair of K0i gives one KIi (Figure 12(c)). For each population j, the coupling is governed by F(vj), a function representing the incoming pulses originating from the connected populations. As shown above, the dynamics of the system’s basic K0 population is described by a 2nd order differential equation with a driving term (i.e., F(vj)). If we Laplace transform the ODE while assuming that the driving force F(vj) is a unit impulse δ(t), we get ab=[s2+(a+b)s+ab]H(s), with H(s) - transfer function. The inverse transform of the transfer function, the UIR is the dual exponential (also used in other models such as Lopes da Silva et al. [35]).(3) KII: this type of population arises when aKIe and KIi are coupled (Figure 12(c)). The distance between KII populations can be large enough that conduction delays cannot be neglected. These delays are indicated by L in Figure 12(c).(4) KIII: when KII sets are interconnected, a KIII set is created. Note that there is just one such a set in Figure12(c). In Freeman [39] there are three coupled KII sets: the olfactory bulb (OB), the anterior olfactory nucleus (AON), and the prepyriform cortex (PC). The connections in this set are longer range so that lags (latencies) were incorporated. These latencies are included in his expressions of the driving force by a backward counter J.  For example, the effect of a subpopulation from AON to one in OB is included as
(53)(19)∑J=1422PE1(J).
This term represents the weighted sum of the AON population’s pulse values PE1 in the past 14–22 steps. Note that the coefficient is divided by 9, the number of backward steps that are included in the sum.Simulations of Freeman’s model were implemented in Matlab and typical results are depicted in the bottom panel in Figure12(c). To show that the system is initially quiet, we deliver a first stimulus a little after onset of the simulation and to demonstrate that the response of the active intact network is slightly different, we deliver another stimulus at about the halfway point of the simulation (arrows in Figure 12(c)).


## 4.4. Model for the Electroencephalogram

The electrical signal that can be recorded from the scalp, the electroencephalogram (EEG), reflects the currents generated by the brain’s nerve cells. If we ignore electrical signals from other more remote sources (e.g., muscle, heart) and nonbiological artifacts, we can assume that the EEG is the weighted sum of the underlying neural network in which the weights are determined by the geometry. The EEG signal is used both in research and clinical settings. Because a single EEG signal includes the activity of millions of nerve cells, the relationship between smaller networks of the brain and the EEG signal is not necessarily a trivial one. Therefore this relationship has been the target of many modeling studies; the framework developed by physicist Paul Nunez is an example. Because it is generally thought that the slower synaptic potentials, especially those of the neocortical pyramidal cells, contribute most to the EEG signals on the surface of cortex and scalp, Nunez’ model focused on describing synaptic activity across cortex (e.g., [31, 32]). The spatiotemporal model of cortical activity is based on the schematic of neocortex shown in Figure 13. This diagram shows the effect of action potential firing activity g produced by the left cube on the synaptic activity h of the right cube. The synaptic activity has subscripts e and i for the excitatory and inhibitory synapses, respectively. Both activity variables g and h are a function of location: r→1 (left cube), r→ (right cube), and time t. The remaining variables in Figure 13 are v for conduction velocity and u(r→,t) for external, subcortical input. Consequently the delay for action potentials arriving at the right cube that were generated in the left cube is |r→-r→1|/v; therefore, the action potential firing function from the left cube arriving at the right cube is indicated as g(r→1,t-|r→-r→1|/v). Furthermore, the connectivity within the cortex is governed by functions RE and RI for the excitatory and inhibitory connections, respectively. These functions also depend on the locations r→, r→1, and fiber system (with conduction velocity v), that is, RE=RE(r→,r→1,v). The equation that describes the excitatory synaptic activity as a function of the input to the right cube can now be formulated as
(54)hE(r→,t)=u(r→,t)+∫0∞dv∫SRE(r→,r→1,v)g(r→1,t-|r→-r→1|v)d2r1.
For the inhibitory synaptic activity a simpler expression can be used since there are only local inhibitory fibers all with similar conduction velocity, that is, RI=RI(r→,r→1). Further, because of the local character of inhibition, one can neglect the delays for the inhibitory fibers: that is, now the action potential firing function is simply g(r→1,t). For the inhibitory expression, we get
(55)hI(r→,t)=u0+∫SRI(r→,r→1)g(r→1,t)d2r1.Figure 13
The model presented by Nunez is based on interaction between cortical volume units. The outputg of one volume at position r→1 arrives at its target at position r→ with a delay determined by conduction velocity (v). Every volume unit is characterized by its synaptic activity, both for excitatory (he) and inhibitory synapses (hi). In the mammalian brain, the excitatory cortical synaptic activity is assumed to make the principal contributions to the electroencephalogram (EEG).The external inputu0 represents subcortical input which is assumed constant for a given state. Because the physiological basis for this input is not obvious, Nunez dropped this term in later versions of the model. In the following, it is not very important since the effect of u0 vanishes when linearizing the expression. Note that both in (54) and (55) the connectivity functions RE and RI and/or activity function g must contain a component that translates the action potential rate into a level of synaptic activity. In his development, Nunez assumes this component to be linear by using a simple gain/attenuation factor. The same can be said for external input functions u(r→,t) and u0. They implicitly include a conversion from action potential rate to synaptic activity; the external input in terms of the action potential rate must be corrected to generate synaptic activity and (again) if this relationship is linear, one can use a simple factor. The above expressions in (54) and (55) relate spike input to synaptic activity (pulse-to-wave conversion); of course there is also the effect of the synaptic activity on the spike output (wave-to-pulse conversion). The nonlinear relationship for both synaptic activities hE and hI on action potential rate function g is commonly modeled by a sigmoid function (e.g., [32]). Nunez linearizes this relationship about an assumed fixed state g0. If we consider a small change of the action potential rate δg around state g0 and relate this to small changes in the synaptic activities δhE and δhI, we get δg=(∂g/∂hE)g=g0δhE+(∂g/∂hI)g=g0δhI. Defining HE=δhE and HI=δhI, QE=(∂g/∂hE)g=g0, and QI=-(∂g/∂hI)g=g0, we can simplify notation and get δg=QEHE-QIHI. If we now use (54) to write an expression for δhE, we get(56)δhE(r→,t)︸HE(r→,t)=δu(r→,t)︸U(r→,t)+∫0∞dv×∫SRE(r→,r→1,v)δg(r→1,t-|r→-r→1|v)︸QEHE(r→1,t-|r→-r→1|/v)-QIHI(r→1,t-|r→-r→1|/v)d2r1,
with the changes of notation indicated by the horizontal curly brackets, we obtain (Equation (11.5) in [32]):
(57)HE(r→,t)=U(r→,t)+∫0∞dv∫SRE(r→,r→1,v)[QEHE(r→1,t-|r→-r→1|v)-QIHI(r→1,t-|r→-r→1|v)]d2r1.
Following the same procedure for (55) gives (Nunez’ Equation (11.6))
(58)HI(r→,t)=∫SRI(r→,r→1)[QEHE(r→1,t)-QIHI(r→1,t)]d2r1.
One of the strong aspects of Nunez’ model is that it acknowledges different conduction velocities, the weak part is that it linearizes the relationship between synaptic activity and action potential rate, and vice versa.


## 4.5. A Field Equation for the EEG

An excellent summary of macroscopic modeling is provided in Jirsa and Haken [45]. These authors integrate temporal and spatial aspects, the firing rate models first described by Wilson and Cowan [29, 30], and the synaptic fields described by Nunez (e.g., [31, 32]). Their approach is straightforward. First they develop expressions for the synaptic and action potential activities. Then they simplify by zooming out, determine the expression for the synaptic field, and write that expression in the form of a convolution operation. Finally, they take the Fourier transform of the convolution and transform this result back into the spatiotemporal domain as a partial differential equation for the synaptic field of the excitatory neurons. Just as Nunez (Section 4.4), Jirsa and Haken assume that this field is a representative for the EEG.


### 4.5.1. Spatiotemporal Expression for the Synaptic Field

The excitatory (e) and inhibitory (i) populations each have pulse rate E and I  (E(x,t) and I(x,t)) similar to the approach by Wilson and Cowan [29], and a wave expression ψ(ψe(x,t) and ψi(x,t)) similar to the development by Nunez. In the above, the spatial component x is three dimensional; in the remainder of the text, the spatial variable is reduced to one dimension represented by x. In the following text, we follow the derivation of the macroscopic field equation as it is visualized in Figure 14. In order to find the expression for the synaptic activity, we determine the synaptic component (top, Figure 14), while including the spike rate generator (bottom, Figure 14). The left part of Figure 14 shows the details whereas the right part shows the zoomed out version representing a macroscopic view where some detail is lost.Figure 14
Diagrams of the model described by Jirsa and Haken [45]. They explicitly model the conversion between cellular output and the synaptic input (transformation from pulse-to-wave, upper panels) and the conversion between synaptic activity and spike rate (transformation from wave-to-pulse, bottom panels). Top-left panel shows how the excitatory synaptic field ψe (as a function of time t and place x) is generated by the action potential field E by integration over the 1-dimensional cortical “area” Γ using integration variable X. The outcome (59) depends on the distribution of connectivity fe and the conversion operation He (60). The other panels in the top row show the same picture in various degrees of zooming out. The Bottom-left panel shows the generation of the action potential field E (as a function of time t and place x) by the excitatory synaptic field ψe (61). The outcome depends on the distribution of connectivity fE and the conversion operation HE (62). The other panels in the bottom row show the same in various degrees of zooming out. Note that the distribution functionfE approaches a delta function δ(x-X). The depicted relationships only show the excitatory population; for the inhibitors, similar diagrams could be constructed.(1)Synaptic Component. The synaptic activity ψ  resulting from the incoming spike activity, for the excitatory population is:
(59)ψe(x,t)︸synapticactivityexcitatorypopulationat  locationxand  timet=∫Γ︸SurfaceAreaof  thebraindXfe(x,X)︸Distribution  Functionof  Spatial  Connectivityof  axons(output)  atXtoxHe(x,X,t)︸Output  ofConversionPulse→Wave.
For the inhibitory population a similar expression can be presented. In the equations for synaptic activity, the contribution from the pulse activity over the whole cortex is included: at location x, this contribution is the product of all converted cortical activity H multiplied by the connection distribution f for x.  Synaptic activity ψ  is then computed by integrating (using integration variable X) over cortical surface Γ (Figure 14). The conversion function He determines how incoming action potentials originating from excitatory cells are converted into excitatory synaptic activity (excitatory post synaptic potentials). This is represented by a sigmoid function S that can be linearized over most of the trajectory using the constant slope ae(60)He(x,X,t)≈aeE(X,t-|x-X|ve).
Note that the E(⋯) function describes all action potentials generated in the past that could have traveled to location x originating from location X (Figure 14, top row) in a time interval |x-X|/ve, with ve—conduction velocity. The ones that have really arrived at x depend on the connectivity between x and X; this is described by function fe(x,X).(2)Spike Rate Component. The expression for the spike rate for the excitatory population E is similar to (59)
(61)E(x,t)︸spikedensityexcitatorypopulationat  locationxand  timet=∫Γ︸SurfaceAreaof  thebraindXfE(x,X)︸Distribution  Functionof  Spatial  Connectivityof  dendrites (input)atXtoxHE(x,X,t)︸Output  ofConversionWave→Pulse.
The output conversion HE: wave-to-pulse becomes an expression that reflects the conversion from: (a) the excitatory synapses, (b) the inhibitory synapses, and (c) a term for the extracortical activity. The synaptic activities for the three wave sources are simply added and a (static) nonlinearity Se converts it to pulse rate. Like in the expression in (60), there is a temporal delay term |x-X|/v that corrects for conduction velocity. The expression for the wave-to-pulse conversion becomes
(62)HE(x,X,t)︸ConversionWave→Pulse=Se(ψe(X,t-|x-X|ve)︸Contribution  of  theExcitatory  SynapticActivities  fromlocationXat timest-delay-ψi(X,t-|x-X|vi)︸Contribution  of  theInhibitory  SynapticActivities+pe(X,t-|x-X|v)︸Contribution  of  theExternal  SynapticActivities{ψe(X,t-|x-X|ve)︸Contribution  of  theExcitatory  SynapticActivities  fromlocationXat timest-delay}).The inclusion of a conduction delay is correct when assuming active propagation of synaptic input, for passive propagation (at the speed of light) such a delay is not needed. For the remainder of the derivation by Jirsa and Haken [45] this point is irrelevant because these delays will be neglected. The conversion function HE describes the contributions from all cortical areas X to location x and the distribution function fE in (61) determines the strength of the connection.


### 4.5.2. Zooming out

In the model, cortical connectivity is described by exponential distributions of the form(1/2σ)e-|x|/σ, functions where the area under the curve evaluates to unity: ∫-∞∞(1/2σ)e-|x|/σdx=1. The advantage is that when we zoom out to a global overview of the cortex, these distributions can be replaced by delta functions. In the model, the distributions fi, fI, and fE are represented by a delta function because they are local (see for example fE in Figure 14 bottom row). This procedure for fE is depicted in Figure 14, bottom row: fE(x,X)≈δ(X-x). Distribution fe, describing the excitatory connections over a longer range, remains in the form of the exponential expression (Figure 14, top row).(1)The Expression forψe. Now we use the above to formulate the expression for the synaptic activity because this variable is considered proportional to the EEG. In this approach, we focus on solving for the excitatory component since this is considered the major contributor to the EEG signal. First we substitute the linearized version of HE (60) into ψe (59):
(63)ψe(x,t)=∫ΓdXfe(x,X)aeE(X,t-|x-X|ve).
In this expression, we substitute (61) for E(⋯), while changing integration variable X to X1, and we get ψe(x,t)=∫ΓdXfe(x,X)ae∫ΓdX1fE(X,X1)HE(X,X1,t-|x-X|/ve). Subsequently we substitute the expression for HE from (62) and replace fE with a delta function
(64)ψe(x,t)=∫ΓdXfe(x,X)ae∫ΓdX1δ(X-X1)×Se(ψe(X1,t-|X-X1|ve-|x-X|ve)-ψi(X1,t-|X-X1|ve-|x-X|ve)+pe(X1,t-|X-X1|ve-|x-X|ve)).
Now we apply the sifting property to evaluate the 2nd integral and get
(65)ψe(x,t)=ae∫ΓdXfe(x,X)Se(ψe(X,t-|x-X|ve)-ψi(X,t-|x-X|ve)+pe(X,t-|x-X|ve)).
For the inhibitory population we follow a similar procedure and get
(66)ψi(x,t)=ai∫ΓdXfi(x,X)Si(ψe(X,t-|x-X|vi)-ψi(X,t-|x-X|vi)+pi(X,t-|x-X|vi)).
If we now evaluate the integral in this expression, while replacing the distribution function fi with the delta function δ(X-x), so we can use the sifting property, we obtain the result for the inhibitory synaptic activity: ψi(x,t)=aiSi(ψe(x,t)-ψi(x,t)+pi(x,t)). Subsequently we follow the procedure by Jirsa and Haken [45] and only take the linear aspect of the Si function, using a slope of αi and finally solve for ψi(67)ψi(x,t)=aiαi1+aiαi(ψe(x,t)+pi(x,t)).
Substitution of (67) into (65) gives(68)ψe(x,t)=ae∫ΓdXfe(x,X)Se×{ψe(X,t-|x-X|ve)-aiαi1+aiαiψe(X,t-|x-X|ve)︷ρ~ψe+pe(X,t-|x-X|ve)-aiαi1+aiαipi(X,t-|x-X|ve)︸p(X,t-|x-X|/ve)}.Using the simplifications indicated by the horizontal curly brackets above we have the expression for the excitatory synaptic activity:(69)ψe(x,t)=ae∫ΓdXfe(x,X)Se{ρ~ψe(X,t-|x-X|ve)+p(X,t-|x-X|ve)}.


### 4.5.3. Expression forψe  as a Convolution

The next step is to write the result forψe in the form of a convolution. First, the notation for the integrand in (69) can be simplified
(70)fe(x,X)aeSe{ρ~ψe(X,t-|x-X|ve)+p(X,t-|x-X|ve)}︸ρ(X,t-|x-X|/ve)=fe(x,X)ρ(X,t-|x-X|ve).
Because the excitatory connections cover a larger area its distribution function cannot be replaced by δ. Therefore we use the exponential distribution for fe and get: (1/2σe)e-|x-X|/σeρ(X,t-|x-X|/ve). Although it seems that we make things more complicated, we rewrite this as an integral: (1/2σe)e-|x-X|/σe∫-∞∞ρ(X,T)δ(T-(t-|x-X|/ve))dT. Now we substitute this result and the expression in (70) in (69) and rearrange the order of the terms
(71)ψe(x,t)=∫Γ∫-∞∞δ(t-T-|x-X|ve)12σee-|x-X|/σe︸G(x-X,t-T)×ρ(X,T)dXdT.
The part (1/2σe)e-|x-X|/σeδ((t-|x-X|/ve)-T) can be interpreted as a unit impulse response or Green’s function G.Now we complete this step and write the expression in the form of a convolution with respect to space and time(72)ψe(x,t)=∫Γ∫-∞∞G(x-X,t-T)ρ(X,T)dXdT.


### 4.5.4. From Convolution to the Field Equation

Now we take advantage of the convolution form of (72). With the spatio-temporal Fourier transforms for ψ(x,t), ρ(x,t) and G(x,t), we can write the transform of the convolution as a product in the frequency domain. To summarize we have the following transform pairs:
(73)ψe(x,t)⟺ψe(k,ω),ρ(x,t)⟺ρ(k,ω),G(x,t)⟺ω02+jω0ωve2k2+(ω0+jω)2,withω0=veσe.
Using ψe(x,t)=G(x,t)⊗ρ(x,t)⇔ψe(k,ω)=G(k,ω)×ρ(k,ω), we get:
(74)ψe(k,ω)=ω02+jω0ωve2k2+(ω0+jω)2ρ(k,ω).
Note that in the temporal Fourier transform jω and -ω2 represent the transformed first and second derivatives in the time domain, while the representatives for the first and second spatial derivatives are +jk and -k2. Transforming (74), we complete our last step and get the equation for the spatio-temporal excitatory synaptic activity
(75)(ω02-ve2∂2∂x2)ψe(x,t)+2ω0ψ˙e(x,t)+ψ¨e(x,t)=(ω02+ω0∂∂t)ρ(x,t).
An example of simulated activity is depicted in Figure 15. In this example, we show the model’s response to a half sine wave pulse in the center of a cortical surface. As expected, the evoked activity is damped over both space and time. The result in (75) is also similar to earlier results obtained for linearized cellular and network models. If we ignore the spatial component in (75), we get a second order ODE with a force term, similar to (26) and (44). Although in this case, a nonlinearity is included in the factor ρ of the force term.Simulation of (75) for times 5, 7, 10, and 25 showing the response of a cortical surface to a stimulus pulse in the middle of the surface. The pulse was half a period of a sine wave between times 0 and 10. Parameters used are ω0=0.15; ve=1.52; ρ~=0.9; ae=1.2; αe=1.2.
(a)
t
=
5
(b)
t
=
7
(c)
t
=
10
(d)
t
=
25


## 4.5.1. Spatiotemporal Expression for the Synaptic Field

The excitatory (e) and inhibitory (i) populations each have pulse rate E and I  (E(x,t) and I(x,t)) similar to the approach by Wilson and Cowan [29], and a wave expression ψ(ψe(x,t) and ψi(x,t)) similar to the development by Nunez. In the above, the spatial component x is three dimensional; in the remainder of the text, the spatial variable is reduced to one dimension represented by x. In the following text, we follow the derivation of the macroscopic field equation as it is visualized in Figure 14. In order to find the expression for the synaptic activity, we determine the synaptic component (top, Figure 14), while including the spike rate generator (bottom, Figure 14). The left part of Figure 14 shows the details whereas the right part shows the zoomed out version representing a macroscopic view where some detail is lost.Figure 14
Diagrams of the model described by Jirsa and Haken [45]. They explicitly model the conversion between cellular output and the synaptic input (transformation from pulse-to-wave, upper panels) and the conversion between synaptic activity and spike rate (transformation from wave-to-pulse, bottom panels). Top-left panel shows how the excitatory synaptic field ψe (as a function of time t and place x) is generated by the action potential field E by integration over the 1-dimensional cortical “area” Γ using integration variable X. The outcome (59) depends on the distribution of connectivity fe and the conversion operation He (60). The other panels in the top row show the same picture in various degrees of zooming out. The Bottom-left panel shows the generation of the action potential field E (as a function of time t and place x) by the excitatory synaptic field ψe (61). The outcome depends on the distribution of connectivity fE and the conversion operation HE (62). The other panels in the bottom row show the same in various degrees of zooming out. Note that the distribution functionfE approaches a delta function δ(x-X). The depicted relationships only show the excitatory population; for the inhibitors, similar diagrams could be constructed.(1)Synaptic Component. The synaptic activity ψ  resulting from the incoming spike activity, for the excitatory population is:
(59)ψe(x,t)︸synapticactivityexcitatorypopulationat  locationxand  timet=∫Γ︸SurfaceAreaof  thebraindXfe(x,X)︸Distribution  Functionof  Spatial  Connectivityof  axons(output)  atXtoxHe(x,X,t)︸Output  ofConversionPulse→Wave.
For the inhibitory population a similar expression can be presented. In the equations for synaptic activity, the contribution from the pulse activity over the whole cortex is included: at location x, this contribution is the product of all converted cortical activity H multiplied by the connection distribution f for x.  Synaptic activity ψ  is then computed by integrating (using integration variable X) over cortical surface Γ (Figure 14). The conversion function He determines how incoming action potentials originating from excitatory cells are converted into excitatory synaptic activity (excitatory post synaptic potentials). This is represented by a sigmoid function S that can be linearized over most of the trajectory using the constant slope ae(60)He(x,X,t)≈aeE(X,t-|x-X|ve).
Note that the E(⋯) function describes all action potentials generated in the past that could have traveled to location x originating from location X (Figure 14, top row) in a time interval |x-X|/ve, with ve—conduction velocity. The ones that have really arrived at x depend on the connectivity between x and X; this is described by function fe(x,X).(2)Spike Rate Component. The expression for the spike rate for the excitatory population E is similar to (59)
(61)E(x,t)︸spikedensityexcitatorypopulationat  locationxand  timet=∫Γ︸SurfaceAreaof  thebraindXfE(x,X)︸Distribution  Functionof  Spatial  Connectivityof  dendrites (input)atXtoxHE(x,X,t)︸Output  ofConversionWave→Pulse.
The output conversion HE: wave-to-pulse becomes an expression that reflects the conversion from: (a) the excitatory synapses, (b) the inhibitory synapses, and (c) a term for the extracortical activity. The synaptic activities for the three wave sources are simply added and a (static) nonlinearity Se converts it to pulse rate. Like in the expression in (60), there is a temporal delay term |x-X|/v that corrects for conduction velocity. The expression for the wave-to-pulse conversion becomes
(62)HE(x,X,t)︸ConversionWave→Pulse=Se(ψe(X,t-|x-X|ve)︸Contribution  of  theExcitatory  SynapticActivities  fromlocationXat timest-delay-ψi(X,t-|x-X|vi)︸Contribution  of  theInhibitory  SynapticActivities+pe(X,t-|x-X|v)︸Contribution  of  theExternal  SynapticActivities{ψe(X,t-|x-X|ve)︸Contribution  of  theExcitatory  SynapticActivities  fromlocationXat timest-delay}).The inclusion of a conduction delay is correct when assuming active propagation of synaptic input, for passive propagation (at the speed of light) such a delay is not needed. For the remainder of the derivation by Jirsa and Haken [45] this point is irrelevant because these delays will be neglected. The conversion function HE describes the contributions from all cortical areas X to location x and the distribution function fE in (61) determines the strength of the connection.


## 4.5.2. Zooming out

In the model, cortical connectivity is described by exponential distributions of the form(1/2σ)e-|x|/σ, functions where the area under the curve evaluates to unity: ∫-∞∞(1/2σ)e-|x|/σdx=1. The advantage is that when we zoom out to a global overview of the cortex, these distributions can be replaced by delta functions. In the model, the distributions fi, fI, and fE are represented by a delta function because they are local (see for example fE in Figure 14 bottom row). This procedure for fE is depicted in Figure 14, bottom row: fE(x,X)≈δ(X-x). Distribution fe, describing the excitatory connections over a longer range, remains in the form of the exponential expression (Figure 14, top row).(1)The Expression forψe. Now we use the above to formulate the expression for the synaptic activity because this variable is considered proportional to the EEG. In this approach, we focus on solving for the excitatory component since this is considered the major contributor to the EEG signal. First we substitute the linearized version of HE (60) into ψe (59):
(63)ψe(x,t)=∫ΓdXfe(x,X)aeE(X,t-|x-X|ve).
In this expression, we substitute (61) for E(⋯), while changing integration variable X to X1, and we get ψe(x,t)=∫ΓdXfe(x,X)ae∫ΓdX1fE(X,X1)HE(X,X1,t-|x-X|/ve). Subsequently we substitute the expression for HE from (62) and replace fE with a delta function
(64)ψe(x,t)=∫ΓdXfe(x,X)ae∫ΓdX1δ(X-X1)×Se(ψe(X1,t-|X-X1|ve-|x-X|ve)-ψi(X1,t-|X-X1|ve-|x-X|ve)+pe(X1,t-|X-X1|ve-|x-X|ve)).
Now we apply the sifting property to evaluate the 2nd integral and get
(65)ψe(x,t)=ae∫ΓdXfe(x,X)Se(ψe(X,t-|x-X|ve)-ψi(X,t-|x-X|ve)+pe(X,t-|x-X|ve)).
For the inhibitory population we follow a similar procedure and get
(66)ψi(x,t)=ai∫ΓdXfi(x,X)Si(ψe(X,t-|x-X|vi)-ψi(X,t-|x-X|vi)+pi(X,t-|x-X|vi)).
If we now evaluate the integral in this expression, while replacing the distribution function fi with the delta function δ(X-x), so we can use the sifting property, we obtain the result for the inhibitory synaptic activity: ψi(x,t)=aiSi(ψe(x,t)-ψi(x,t)+pi(x,t)). Subsequently we follow the procedure by Jirsa and Haken [45] and only take the linear aspect of the Si function, using a slope of αi and finally solve for ψi(67)ψi(x,t)=aiαi1+aiαi(ψe(x,t)+pi(x,t)).
Substitution of (67) into (65) gives(68)ψe(x,t)=ae∫ΓdXfe(x,X)Se×{ψe(X,t-|x-X|ve)-aiαi1+aiαiψe(X,t-|x-X|ve)︷ρ~ψe+pe(X,t-|x-X|ve)-aiαi1+aiαipi(X,t-|x-X|ve)︸p(X,t-|x-X|/ve)}.Using the simplifications indicated by the horizontal curly brackets above we have the expression for the excitatory synaptic activity:(69)ψe(x,t)=ae∫ΓdXfe(x,X)Se{ρ~ψe(X,t-|x-X|ve)+p(X,t-|x-X|ve)}.


## 4.5.3. Expression forψe  as a Convolution

The next step is to write the result forψe in the form of a convolution. First, the notation for the integrand in (69) can be simplified
(70)fe(x,X)aeSe{ρ~ψe(X,t-|x-X|ve)+p(X,t-|x-X|ve)}︸ρ(X,t-|x-X|/ve)=fe(x,X)ρ(X,t-|x-X|ve).
Because the excitatory connections cover a larger area its distribution function cannot be replaced by δ. Therefore we use the exponential distribution for fe and get: (1/2σe)e-|x-X|/σeρ(X,t-|x-X|/ve). Although it seems that we make things more complicated, we rewrite this as an integral: (1/2σe)e-|x-X|/σe∫-∞∞ρ(X,T)δ(T-(t-|x-X|/ve))dT. Now we substitute this result and the expression in (70) in (69) and rearrange the order of the terms
(71)ψe(x,t)=∫Γ∫-∞∞δ(t-T-|x-X|ve)12σee-|x-X|/σe︸G(x-X,t-T)×ρ(X,T)dXdT.
The part (1/2σe)e-|x-X|/σeδ((t-|x-X|/ve)-T) can be interpreted as a unit impulse response or Green’s function G.Now we complete this step and write the expression in the form of a convolution with respect to space and time(72)ψe(x,t)=∫Γ∫-∞∞G(x-X,t-T)ρ(X,T)dXdT.


## 4.5.4. From Convolution to the Field Equation

Now we take advantage of the convolution form of (72). With the spatio-temporal Fourier transforms for ψ(x,t), ρ(x,t) and G(x,t), we can write the transform of the convolution as a product in the frequency domain. To summarize we have the following transform pairs:
(73)ψe(x,t)⟺ψe(k,ω),ρ(x,t)⟺ρ(k,ω),G(x,t)⟺ω02+jω0ωve2k2+(ω0+jω)2,withω0=veσe.
Using ψe(x,t)=G(x,t)⊗ρ(x,t)⇔ψe(k,ω)=G(k,ω)×ρ(k,ω), we get:
(74)ψe(k,ω)=ω02+jω0ωve2k2+(ω0+jω)2ρ(k,ω).
Note that in the temporal Fourier transform jω and -ω2 represent the transformed first and second derivatives in the time domain, while the representatives for the first and second spatial derivatives are +jk and -k2. Transforming (74), we complete our last step and get the equation for the spatio-temporal excitatory synaptic activity
(75)(ω02-ve2∂2∂x2)ψe(x,t)+2ω0ψ˙e(x,t)+ψ¨e(x,t)=(ω02+ω0∂∂t)ρ(x,t).
An example of simulated activity is depicted in Figure 15. In this example, we show the model’s response to a half sine wave pulse in the center of a cortical surface. As expected, the evoked activity is damped over both space and time. The result in (75) is also similar to earlier results obtained for linearized cellular and network models. If we ignore the spatial component in (75), we get a second order ODE with a force term, similar to (26) and (44). Although in this case, a nonlinearity is included in the factor ρ of the force term.Simulation of (75) for times 5, 7, 10, and 25 showing the response of a cortical surface to a stimulus pulse in the middle of the surface. The pulse was half a period of a sine wave between times 0 and 10. Parameters used are ω0=0.15; ve=1.52; ρ~=0.9; ae=1.2; αe=1.2.
(a)
t
=
5
(b)
t
=
7
(c)
t
=
10
(d)
t
=
25


## 5. Models with a Stochastic Component

With the exception of the spin model and Hopfield’s approach (Section3.1), all models we discussed are deterministic. Most of the deterministic approaches can easily be modified to include a stochastic component [48]. Such a component may represent intrinsic noise in the membrane and/or the synaptic input originating from sources external to the modeled unit. In the 1960s and 1970s, membrane noise was investigated to quantify electrical channel properties (e.g., [75]). Especially when modeling a network embedded in the brain, its input from the rest of the brain is generated by enormous networks and a stochastic approach to this input is appropriate and also significant because it affects the network function (e.g., [76]). A critical shortcoming of the deterministic models (e.g., the mean field models discussed in Section 4), that can be mitigated with a stochastic approach, is the absence of fluctuations and correlations of the neural activity, especially spike trains (e.g., [77, 78]). Moreover, a statistical approach is applied to spike trains to relate to information coding in networks in general [47, 79]. Mathematically, stochastic modeling often results in a diffusion equation approach, where the models are frequently approximated by either the Fokker-Planck equation and/or Langevin equation (e.g., Appendix in [48]). Further discussion of these approaches is beyond the scope of this paper.Recently, Benayoun et al. [19] and Wallace et al. [20] used a network of switches in which the update rule was asynchronous and stochastic known as the Gillespie’s exact stochastic simulation algorithm [80, 81]. To show this approach, I describe a simplified version
(76)Af~↑↓αQ
The active neurons A are created from an infinite pool of inactive ones Q at a rate f~ and they become inactive again with a rate α. The rate function f~ depends on the activity level in the network, however, for the sake of the example, we will not worry about this “detail” here. Rate α is the intrinsic rate at which active cells become inactive (this also includes the refractory period). The deterministic kinetic rate equation for this process is
(77)dAdt=f~Q-αA=f-αA.
We assume here that Q is available in large excess (i.e., an infinite pool) so that it can be considered constant. This equation can be solved both analytically and numerically.The next step is to consider the so-called master equation for this process. In contrast to the deterministic approach above, we now use a stochastic approach and describe the probability that there aren active neurons (i.e., n cells in our infinite pool are in state A). We may visualize this by a line on which we see all the possible states for the number of active neurons, that is, 0,1,2,3,4,5,…,n-1,n,n+1,n+2,n+3,…,∞. If we update the system fast enough, we can safely assume that a state change only involves a single step.Under this assumption, the probabilities describing the dynamics around staten are easily seen in the following diagram that depicts part of the line of possible states
(78)n+1⋯Pn+1f↑↓αn⋯Pnf↑↓αn-1⋯Pn-1
This leads to the master equation describing the dynamics of the probability Pn for the state with n active neurons
(79)dPndt=α(n+1)Pn+1-αnPn+fPn-1-fPn.
To simplify this further, we can apply shift operators E+ and E- that operate on function f(n) and increase or decrease the number n by one, respectively [82]. This allows one to rewrite (79) as a function of Pn(80)dPndt=[αE+n-αn+fE--f]Pn.
The problem is that in most cases the master equation is difficult to solve analytically or cannot be solved at all. Gillespie [80, 81] formulated a remarkable and exact algorithm to simulate chemical systems in a stochastic fashion. Within this approach, one starts from the system in a particular state at time t and determines the probability density function (PDF) for the next reaction. Briefly, if we consider 1,2,…,μ,…,M reactions; if one simulates such a system over time, we only need to know (a)
when the next reaction (any of theM) will occur, and(b)
what kind of reaction (μ) it will be.In order to resolve these questions, Gillespie follows the system of reactions from any timet to t+τ+dτ.  Accordingly, he formulates the joint probability P(τ,μ)dτ as the probability that none of the reactions occurs in the interval between t and t+τ,  and that reaction μ  does occur in the interval between t+τ  and t+τ+dτ (Figure 16(a)). It is straightforward to show that P(τ,μ)dτ is characterized by the product of both probabilities. Let us first determine the first component P0 that no reaction occurs in the interval between t and +τ. Consider a time axis τ′ set at 0 immediately after one of the reactions occurred. Consequently, the probability that none of the reactions occurs in dτ′ is
(81)1-∑μ=1Maμdτ′.
Thus one can define the probability of no reaction at τ′+dτ′ is the probability of no reaction at τ′ multiplied by the expression in (81)
(82)P0(τ′+dτ′)=P0(τ′)[1-∑μ=1Maμdτ′]
which can be written as
(83)P0(τ′+dτ′)-P0(τ′)P0(τ′)=-∑μ=1Maμdτ′=-a0dτ′,
where a0 is the sum over all aμ. Integration of (83) from 0 to τ results in
(84)P0(τ)=e-a0τ.
Component (b) that reactions μ  do occur in the interval between t+τ  and t+τ+dτ is given by aμdτ. Thus the probability  P(τ,μ)dτ is the product of this result and (84), therefore
(85)P(τ,μ)dτ=aμe-a0τdτ→P(τ,μ)=aμe-a0τ.
For all values of μ (the reaction number) and τ (time steps >0). Now we go back to the initial questions (a) and (b) above, and we state that
(86)P(τ,μ)=P1(τ)︸When  is  the  next  reaction?P2(μ∣τ)︸Givenτwhat  reaction  is  it?.
The first factor indicates the probability for the next reaction and it does not matter which one, therefore
(87)P1(τ)=∑μ=1MP(τ,μ)=∑μ=1Maμe-a0τ=a0e-a0τ.
From substitution of (87) in (85), we have
(88)P2(μ∣τ)=P(τ,μ)∑μ=1MP(τ,μ)=aμa0.
So now we have the probabilities for both τ and μ and the only question that remains is how we use this knowledge to implement a simulation. Usually the particular probability functions in (87) and (88) are not readily available in a simulation environment. But most computers do have a random number generator based on a uniform distribution between 0 and 1. In a simulation, the uniform distribution can be used to generate values taken from any arbitrary PDF by using the uniform distribution to pick random points on the cumulative distribution of the target PDF. The procedures for τ and μ are shown in Figures 16(c) and 16(d). Because τ is a continuous variable, we can compute the cumulative (Cum[⋯]) function by integration of (87)
(89)Cum[P1(τ)]=1-e-a0τ.
In order to determine both the values for τ and μ, we generate two random numbers r(1) and r(2) respectively. Now we can connect the random number r(1) to the cumulative distribution of P1 (Figure 16(c)). However, since the cumulative function ranges between 0–1, we can also connect r(1) directly to exp(-a0τ) (Figure 16(c); note that 0 and 1 are inverted in this case). If we now solve for τ, we get
(90)r(1)=e-a0τ→τ=1a0log1r(1).
Subsequently we create the cumulative distribution for the reactions and use r(2) to update the system (Figure 16(d)). This procedure of picking a time step and update the system can be repeated for the simulation epoch. Examples of two simulations of our simplified system are depicted in Figure 16(b). For the sake of this example (see Matlab script in Appendix B), we update individual switches (neurons) at each time step. In homogenous all-to-all networks, the simulation algorithm of the number of neurons active in each population may be simplified along the lines of Gillespie's original presentation for a well-mixed chemical system, since the downward and upwards transition rates are identical for all neurons. This simplification will result in a faster computation since the time steps will increase (due to a smaller a0) as compared to our approach in the example.Gillespie’s approach for Monte Carlo simulation of the master equation. (a) The joint PDF that determines the reactionsμ as function of time τ. Note that the variable τ is continuous and μ discrete.  (b) shows two examples of deterministic (black lines) and associated stochastic (blue lines) simulations superimposed. The two cases start from a low and high initial value. The inset is a detail of  how the stochastic solutions fluctuate around the deterministic ones. (c) The procedure of picking the time steps. (d) depicts how to determine the cell to update at each time step.
(a)
(b)
(c)
(d)Because the master equation can be used to recover the deterministic rate equation, the stochastic approach is more general than the mean field models we discussed in Section4. We recover the mean field equation by determining the mean number of active neurons n from the master equation via its definition
(91)〈n〉=∑n=0∞nPn.
In our case we assumed an infinite pool of nerve cells so there is no finite limitation to the number of active neurons. Therefore the sum in (91) is taken from 0→∞. Other models have assumed a limited number of inactive and active nerve cells.We now plug in (91) into (79) (i.e., we multiply each term by n and sum over all n), and we get
(92)ddt∑n=0∞nPn=α∑n=0∞n(n+1)Pn+1-α∑n=0∞‍n2Pn+f∑n=0∞nPn-1-f∑n=0∞nPn.
This generates an expression for the time derivative of the mean (the term left of the equal sign). We now substitute n′ and n′′ for n-1 and n+1, respectively and we have
(93)ddt∑n=0∞nPn=α∑n′′=1∞n′′(n′′-1)Pn′′-α∑n=0∞‍n2Pn+f∑n′=-1∞(n′+1)Pn′-f∑n=0∞nPn.
Note that we changed the range over which we sum in the 1st and 3rd term. We now sum again from 0→∞. Note that we can make this change in the summation interval since P-1=0 in the 3rd term, and for n′′=0, n′′(n′′-1)Pn′′=0 in the 1st term. Furthermore, we expand the first and third terms and get
(94)ddt∑n=0∞nPn=α∑n′′=0∞n′′2Pn′′-α∑n′′=0∞n′′Pn′′-α∑n=0∞‍n2Pn+f∑n′=0∞n′Pn′+f∑n′=0∞‍Pn′-f∑n=0∞nPn.
Terms 1 and 3, 4 and 6 cancel and terms 2 and 5 remain, that is
(95)ddt∑n=0∞nPn=-α∑n′′=0∞n′′Pn′′+f∑n′=0∞‍Pn′.
Using (91), and knowing that all probabilities must add up to one, we obtain an expression similar to the kinetic rate equation (77)
(96)ddt〈n〉=f-α〈n〉,
where the mean of n is the variable A. Thus the mean number of active neurons in the master equation of the process produces the kinetic rate equation. Accordingly, the work by Benayoun et al. [19] and Wallace et al. [20] showed that the mean activity generated by the stochastic model in which f depended on the network activity, corresponded to the mean field Wilson-Cowan model (Section 4.2). In addition, the simulations of their stochastic model made predictions about the fluctuations of the neural activity. These fluctuations followed avalanche statistics as observed in real neuronal networks or even in the scalp or intracranial EEG [83–85].A recent, novel approach to solve the master equation for a neuronal network model used the equation of motion of the generator function in operator notation(97)|G〉=∑n=0∞Pn|n〉.
Here, we use the Dirac bra-ket notation (for a didactic overview of the notation see [86]). The ket of G represents the probability Pn for each state n represented by the ket of n in an infinite dimensional (Hilbert) space. The advantage of this approach is that the extensive mathematics developed for the operator notation can now be used to analyze neuronal networks (e.g., [49, 87]). Using (79), multiplying with ket of n and summing over n, we get
(98)ddt∑n=0∞‍Pn|n〉=α∑n=0∞‍(n+1)Pn+1|n〉-α∑n=0∞‍nPn|n〉+f∑n=0∞‍Pn-1|n〉-f∑n=0∞‍Pn|n〉.
Now we use the equality in (97) and  n′ and n′′ for n-1 and n+1, respectively and rewrite the expression in (98)
(99)ddt|G〉=α∑n′′=1∞‍(n′′)Pn′′|n′′-1〉-α∑n=0∞‍nPn|n〉+f∑n′=-1∞‍Pn′|n′+1〉-f|G〉.
For the same reasons as we did above, we summate n′ and n′′ from 0 to ∞ and, in addition, we use the raising and lowering operators a+ and a-(100)a+|n〉=|n+1〉,a-|n〉=n|n-1〉.
Note that the eigenvalues are 1 and n instead of the conventional √(n+1) and √n.Now we get(101)ddt|G〉=α∑n′′=0∞‍Pn′′a-|n′′〉-α∑n=0∞‍Pna+a-|n〉+f∑n′=0∞‍Pn′a+|n′〉-f|G〉.
Using (97), we can simplify this to
(102)ddt|G〉=αa-|G〉-αa+a-|G〉+fa+|G〉-f|G〉.
With a bit of algebra, we may rewrite this, and recognize it as a linear equation of motion
(103)ddt|G〉=-(a+-1)(αa--f)|G〉→|G˙〉=-L|G〉.
If we use the eigenfunctions |j〉 of L, to expand |G〉(104)|G〉=∑jGj|j〉.
For each eigenfunction there is a corresponding eigenvalue j of L, and therefore
(105)G˙j=-jGj.
The solution of this differential equation is
(106)Gj(t)=Kje-jt.
Now we can retrieve the PDF by projecting G on bra-n, and using the expressions in (104)–(106)
(107)Pn(t)=〈n∣G(t)〉=〈n|∑jGj(t)|j〉=∑jKje-jt〈n∣j〉.
The Kj coefficients are determined from the initial distribution Pn(0).


## 6. Epileptic Seizures

Epilepsy is a serious neurological disorder characterized by spontaneous recurrent seizures. In the electroencephalogram (EEG) of patients with epilepsy, one may observe seizures (ictal events) and interictal events: for example, epileptic spikes or spike-waves. The real problem is that about 30% of the current population of 60 million patients with epilepsy do not respond to any treatment [88]. This problem is mostly due to an incomplete understanding of the mechanisms that underlie this pathology (e.g., [89]). As is the case for many other neurological diseases, this directly relates to a poor understanding of network function in general. Because of lack of experimental tools for studying network behavior at sufficient scale with the associated detail (e.g. [70, 90]), there is a significant role for modeling in this field (e.g., [3, 5, 6, 27, 41, 61, 91–102]).


### 6.1. Modeling Epileptic Activity

Although there are many clinical types of seizure [104], the seizure activity in the EEG is often characterized by rhythmic activity usually lasting for minutes. Examples of typical seizure onset activity are depicted in Figure 17. In contrast, interictal spikes and waves are sudden, short bursts of electrical activity. Although there are exceptions, in between the epileptiform events (i.e., seizures, spikes, waves), the EEG of patients with epilepsy cannot be distinguished from normal EEG patterns. Because epilepsy is characterized by its EEG patterns representing large networks, it is impossible to define epileptiform activity at the level of single neurons or even small networks. Although there seem to be cellular activity patterns frequently found in epileptic brain tissue, such as the paroxysmal depolarization shift (PDS), it is by no means clear how these PDSs relate to the typical network phenomena in epilepsy. The oscillatory activity of cells and small networks is of general interest to the neuroscientist but equally of interest in epilepsy research. However, at this level, the distinction between normal and pathologic oscillatory activity is undefined. At the scale of large networks, as reflected in the EEG signal, the distinction between normal and abnormal is not without problems but definitely easier to make. This makes the study of networks at mesoscopic and macroscopic scales of particular interest in epilepsy research. Based on the activity patterns observed in patients with epilepsy, an effective network model should be able to explain(1)
normal and epileptiform states,(2)
prototypical oscillations and bursts, and(3)
transitions between the states.Can seizure onset be modeled with a Hopf bifurcation? (a) A partial complex seizure spreads and grows in amplitude, suggesting one might model it with a supercritical Hopf bifurcation. (b) A generalized seizure occurs suddenly and at full amplitude, similar to a subcritical Hopf bifurcation. (c) Note that the same seizure as in (b) can also be explained by a perturbation in a bistable regime associated with the subcritical Hopf bifurcation. In this case the transition into the seizure state is not a bifurcation but due to a change of basin of attraction (e.g., [103]).
(a)
(b)
(c)From these properties, a few general conclusions about the nature of the model can be determined. First, a model of network activity must be at least two dimensional to explain oscillations.  Second, transitions between normal and epileptiform activity can be modeled by a bifurcation in a nonlinear dynamical system. Of course, there are many bifurcations that could explain sudden onset and offset of seizure-like oscillations. Examples using the simplest candidate: co-dimension-1 bifurcations transitioning a system between steady state and oscillatory behavior, the sub- and supercritical Hopf bifurcations [105], are shown in Figure 17. In Figure 17, the two top panels show the bifurcation superimposed on seizure activity, the bottom panel indicates that in case of the sub-critical scenario, the transition can also occur if the system exceeds the boundary between the basins of attractors (one attractor being the stable node and the other the stable limit cycle) [103].Some of the models employed in epilepsy research go beyond the minimalistic approach of two dimensions. The detailed simulation models with compartmental neuron models with realistic ion channels may have hundreds of parameters. Independent of the level of complexity of the network models, it is important to know if the neuron model itself is capable of oscillatory activity so that the individual network node can be a pacemaker for oscillation.


### 6.2. Population Models and Epileptiform Activity

Here we examine the population models we described in Section4 for their capability of generating seizure-like behavior. The Wilson and Cowan [29] model is both nonlinear and two-dimensional and is therefore capable of showing a range of behaviors that include limit cycles (Figure 10) and transitions from stable behavior to an oscillatory one. The Wilson and Cowan model satisfies the three criteria of an effective network model for studying epilepsy described in the previous section. In addition, it is also a minimal model because the two dimensions (E and I, (35), (36)) are minimally required to include oscillation and the nonlinearity (sigmoids Se and Si) are essential for the existence of sudden transitions (bifurcations) between nonoscillatory and oscillatory states. An example of how a gradually changing excitatory-coupling (WE, represented by constants c1 and c3 in Figure 9) can cause sudden transitions of the population activity is depicted in Figure 18. Interestingly, it can be seen that, depending on the initial state, oscillations may occur with both increasing and decreasing strength of the coupling (arrows, Figure 18).Figure 18
Bifurcations, caused by a gradually varying excitatory coupling (blue, top trace), creating a seizure-like epoch in population activity governed by the Wilson and Cowan Equations (Excitatory population–red and inhibitory population–green).The modeling approach by Lopes da Silva et al. [35] was used to formulate the Jansen neural mass model, which was later examined for epileptic behavior by Grimbert and Faugeras [38]. Several versions of neural mass models were recently investigated [41, 91, 93, 98, 102]. Because these neural mass models are at  least two dimensional and nonlinear (see for example (52)), they are capable of displaying transitions between activities associated with a stable node and oscillations or even chaotic behavior if they have three or more dimensions. The neural mass models by Freeman [39] are capable of displaying these behaviors as well since they are both relatively high-dimensional and nonlinear. However, due to Freeman’s parameter choices, the isolated K0 populations (Section 4.3.2) show stable behavior only. The nondriven K0 system obeys the ODE v¨j-τv˙j+Δvj=0. In this notation, τ is the sum and Δ the product of the eigenvalues. In Freeman’s simulations the eigenvalues are −0.72 and −0.22, thus their sum τ=-0.94 and product Δ=0.158. Since τ is negative, Δ positive, and τ2-4Δ=0.25>0, the K0 units have the behavior of a stable node (e.g., [105, Chapter 5]). Therefore, any behavior in the model other than that associated with a stable node, must be generated externally by the driving force F, that is, the external input and interconnections between the populations. Using Freeman’s choice of parameters, the model displays chaotic behavior similar to ongoing discharges during a seizure (Figure 12(c)).Linearized models do not support bifurcations and therefore fall short in explaining a critical aspect of epilepsy (seizure onset and offset). The linearized model of Nunez [31, 32] has this shortcoming. Nonetheless, it is useful to check if it supports oscillatory activity. First, the oscillatory activity can model ongoing seizure activity. Second, the point where an unstable oscillation arises in a linearized system is an indicator for a Hopf bifurcation in the full nonlinear model. Models in neuroscience are often represented by electronic equivalent circuits (e.g., Figures 1, 2, 3, and 5). In contrast, Nunez compares the spatiotemporal behavior of his cortical model with that of the mechanics of a string (e.g., [32]). If we consider a string in vacuum attached by its ends and subject to a tension T (but unattached to any other structure, Figure 19(a)), we only need to consider forces acting in the string itself. At rest the string is straight but now we perturb the string a little in the vertical direction y. Let us analyze the forces by considering a small part of the string of length Δx. In the detail of Figure 19(a), we can see that the forces acting on each side of this small piece are equal to the tension T in the string. Because of the curve adopted by the string, the two forces at each side of Δx are not exactly opposed and therefore they do not cancel. First we consider the horizontal component of the force: FH=Tcosα2-Tcosα1. This horizontal component can be neglected since, for small values of the angles α1 and α2, we have cosα1≈1 and cosα2≈1. Therefore we find that FH≈0. The upward force is Tsinα2 and the downward force is Tsinα1. The vertical force FV acting on Δx is
(108)FV=Tsinα2-Tsinα1.
We divide this equation by T and (recall that cosα1≈1 and cosα2≈1) and we get
(109)FVT=Tsinα2Tcosα2-Tsinα1Tcosα1=tanα2-tanα1.
According to Newton’s second law, the net force must equal the product of the mass of Δx and its net acceleration ∂2ϕ/∂t2. If we define the mass per unit length of the string by m, the mass of our piece is mΔx. Using the above for FV, using (109), and rearranging terms, we get
(110)1Δx(tanα2-tanα1)=mT∂2ϕ∂t2.
Because the tangents can be written as the slope of the curve ϕ adopted by the string, we can rewrite the expression left of the equal sign as
(111)1Δx(tanα2-tanα1)=1Δx(∂φ∂x|x2-∂φ∂x|x1),
and for Δx→0, this expression becomes
(112)∂2ϕ∂x2.
Using this result and substituting propagation speed v=T/m in (110), we obtain the wave equation
(113)∂2ϕ∂t2=v2∂2ϕ∂x2.
Now we add a local spring connected to our little piece of string, an extra external force and we embed our mass of the string in a damping environment (e.g., molasses); this situation is depicted in Figure 19(b). Now we have to include these forces in Newton’s law for FV in (113)
(114)m∂2ϕ∂t2=T∂2ϕ∂x2+FLocal+FDamping+FExternal.
The local force due to the spring is proportional to the position of the spring ϕ; that is FLocal=-Kϕ, K is a spring constant and there is a minus sign because a positive deflection of the piece of string causes a force in the opposite direction. The damping force is proportional to the vertical velocity of the piece of string ∂ϕ/∂t; also in this case there is a minus sign because a damping force works in a direction opposite of the velocity: FDamping=-B(∂ϕ/∂t). Thus both the local spring and damping forces are restoring forces. Putting this all into (114) yields.
(115)(Km-v2∂2∂x2)ϕ+Bm∂ϕ∂t+∂2ϕ∂t2=1mFExternal.
It can be seen that this result and the one from Jirsa and Haken [45], (75), are very similar. Equation (75) explicitly includes a nonlinear component in the force term, whereas FExternal here may be either linear or nonlinear.  The result in (115) is an expression that can explain multiple types of behaviors. If B (damping) is very large the equation starts to resemble the diffusion equation ∂ϕ/∂t=D(∂2ϕ/∂x2) whereas for very small B the equation resembles more the wave equation ∂2ϕ/∂t2=v2(∂2ϕ/∂x2) [32]. If on the other hand, the spatial derivative and the external force term dominate, the expression becomes similar to Poisson’s Equation ∇2ϕ=f.Diagram of string models used for EEG. (a) A string attached at its ends but unattached to local structures; the associated detail shows the forces acting on a small pieceΔx between x1 and x2. Forces acting on each side of this small piece of string are equal to the tension T in the string. The curve ϕ adopted by the perturbed string is associated with angles α1 and α2 at x1 and x2respectively. (b) shows a detail of the same string but now there is an external force, a local effect (a local spring), and the string is embedded in a damping medium (molasses).
(a)
(b)Finally we examine the result (75) from Jirsa and Haken [45] for a special case, a generalized seizure (e.g., Figures 17(b) and 17(c)). Since we have generalized activity we assume that there is spatial homogeneity, and since the seizure sustains itself we assume there is no external input; due to these assumptions, the partial derivative term ∂2ψe/∂x2 vanishes and so does the external input part in the ρ variable (see (70) for the definition of ρ). The Taylor series approximation of this variable is ρ(x,t)≈aeαeρ~ψe(x,t). Substituting this approximation into (75) and using the above considerations, we can simplify
(116)ψ¨e(t)+(2ω0-ω0aeαeρ~)︸bψ˙e(t)+(ω02-ω02aeαeρ~)︸cψe(t)=0.
The general solution of this ODE is in the form A1e-λ1t+A2e-λ2t, with λ1,2=-b/2±b2-4c/2. Therefore, we may expect intrinsic global oscillations only if λ1,2 are imaginary, that is if the condition b2-4c<0 is satisfied. Now note that in the simplified equation above c=ω0(b-ω0), if we substitute this into the condition above: b2-4c=b2-4bω0+4ω02=(b-2ω0)2. Since the outcome is a quadratic expression, the condition b2-4c<0 will never be satisfied and therefore, without oscillatory input (e.g., an external network or a set of pacemaker neurons in the network), oscillations cannot occur in this generalized and linearized version of the model.


## 6.1. Modeling Epileptic Activity

Although there are many clinical types of seizure [104], the seizure activity in the EEG is often characterized by rhythmic activity usually lasting for minutes. Examples of typical seizure onset activity are depicted in Figure 17. In contrast, interictal spikes and waves are sudden, short bursts of electrical activity. Although there are exceptions, in between the epileptiform events (i.e., seizures, spikes, waves), the EEG of patients with epilepsy cannot be distinguished from normal EEG patterns. Because epilepsy is characterized by its EEG patterns representing large networks, it is impossible to define epileptiform activity at the level of single neurons or even small networks. Although there seem to be cellular activity patterns frequently found in epileptic brain tissue, such as the paroxysmal depolarization shift (PDS), it is by no means clear how these PDSs relate to the typical network phenomena in epilepsy. The oscillatory activity of cells and small networks is of general interest to the neuroscientist but equally of interest in epilepsy research. However, at this level, the distinction between normal and pathologic oscillatory activity is undefined. At the scale of large networks, as reflected in the EEG signal, the distinction between normal and abnormal is not without problems but definitely easier to make. This makes the study of networks at mesoscopic and macroscopic scales of particular interest in epilepsy research. Based on the activity patterns observed in patients with epilepsy, an effective network model should be able to explain(1)
normal and epileptiform states,(2)
prototypical oscillations and bursts, and(3)
transitions between the states.Can seizure onset be modeled with a Hopf bifurcation? (a) A partial complex seizure spreads and grows in amplitude, suggesting one might model it with a supercritical Hopf bifurcation. (b) A generalized seizure occurs suddenly and at full amplitude, similar to a subcritical Hopf bifurcation. (c) Note that the same seizure as in (b) can also be explained by a perturbation in a bistable regime associated with the subcritical Hopf bifurcation. In this case the transition into the seizure state is not a bifurcation but due to a change of basin of attraction (e.g., [103]).
(a)
(b)
(c)From these properties, a few general conclusions about the nature of the model can be determined. First, a model of network activity must be at least two dimensional to explain oscillations.  Second, transitions between normal and epileptiform activity can be modeled by a bifurcation in a nonlinear dynamical system. Of course, there are many bifurcations that could explain sudden onset and offset of seizure-like oscillations. Examples using the simplest candidate: co-dimension-1 bifurcations transitioning a system between steady state and oscillatory behavior, the sub- and supercritical Hopf bifurcations [105], are shown in Figure 17. In Figure 17, the two top panels show the bifurcation superimposed on seizure activity, the bottom panel indicates that in case of the sub-critical scenario, the transition can also occur if the system exceeds the boundary between the basins of attractors (one attractor being the stable node and the other the stable limit cycle) [103].Some of the models employed in epilepsy research go beyond the minimalistic approach of two dimensions. The detailed simulation models with compartmental neuron models with realistic ion channels may have hundreds of parameters. Independent of the level of complexity of the network models, it is important to know if the neuron model itself is capable of oscillatory activity so that the individual network node can be a pacemaker for oscillation.


## 6.2. Population Models and Epileptiform Activity

Here we examine the population models we described in Section4 for their capability of generating seizure-like behavior. The Wilson and Cowan [29] model is both nonlinear and two-dimensional and is therefore capable of showing a range of behaviors that include limit cycles (Figure 10) and transitions from stable behavior to an oscillatory one. The Wilson and Cowan model satisfies the three criteria of an effective network model for studying epilepsy described in the previous section. In addition, it is also a minimal model because the two dimensions (E and I, (35), (36)) are minimally required to include oscillation and the nonlinearity (sigmoids Se and Si) are essential for the existence of sudden transitions (bifurcations) between nonoscillatory and oscillatory states. An example of how a gradually changing excitatory-coupling (WE, represented by constants c1 and c3 in Figure 9) can cause sudden transitions of the population activity is depicted in Figure 18. Interestingly, it can be seen that, depending on the initial state, oscillations may occur with both increasing and decreasing strength of the coupling (arrows, Figure 18).Figure 18
Bifurcations, caused by a gradually varying excitatory coupling (blue, top trace), creating a seizure-like epoch in population activity governed by the Wilson and Cowan Equations (Excitatory population–red and inhibitory population–green).The modeling approach by Lopes da Silva et al. [35] was used to formulate the Jansen neural mass model, which was later examined for epileptic behavior by Grimbert and Faugeras [38]. Several versions of neural mass models were recently investigated [41, 91, 93, 98, 102]. Because these neural mass models are at  least two dimensional and nonlinear (see for example (52)), they are capable of displaying transitions between activities associated with a stable node and oscillations or even chaotic behavior if they have three or more dimensions. The neural mass models by Freeman [39] are capable of displaying these behaviors as well since they are both relatively high-dimensional and nonlinear. However, due to Freeman’s parameter choices, the isolated K0 populations (Section 4.3.2) show stable behavior only. The nondriven K0 system obeys the ODE v¨j-τv˙j+Δvj=0. In this notation, τ is the sum and Δ the product of the eigenvalues. In Freeman’s simulations the eigenvalues are −0.72 and −0.22, thus their sum τ=-0.94 and product Δ=0.158. Since τ is negative, Δ positive, and τ2-4Δ=0.25>0, the K0 units have the behavior of a stable node (e.g., [105, Chapter 5]). Therefore, any behavior in the model other than that associated with a stable node, must be generated externally by the driving force F, that is, the external input and interconnections between the populations. Using Freeman’s choice of parameters, the model displays chaotic behavior similar to ongoing discharges during a seizure (Figure 12(c)).Linearized models do not support bifurcations and therefore fall short in explaining a critical aspect of epilepsy (seizure onset and offset). The linearized model of Nunez [31, 32] has this shortcoming. Nonetheless, it is useful to check if it supports oscillatory activity. First, the oscillatory activity can model ongoing seizure activity. Second, the point where an unstable oscillation arises in a linearized system is an indicator for a Hopf bifurcation in the full nonlinear model. Models in neuroscience are often represented by electronic equivalent circuits (e.g., Figures 1, 2, 3, and 5). In contrast, Nunez compares the spatiotemporal behavior of his cortical model with that of the mechanics of a string (e.g., [32]). If we consider a string in vacuum attached by its ends and subject to a tension T (but unattached to any other structure, Figure 19(a)), we only need to consider forces acting in the string itself. At rest the string is straight but now we perturb the string a little in the vertical direction y. Let us analyze the forces by considering a small part of the string of length Δx. In the detail of Figure 19(a), we can see that the forces acting on each side of this small piece are equal to the tension T in the string. Because of the curve adopted by the string, the two forces at each side of Δx are not exactly opposed and therefore they do not cancel. First we consider the horizontal component of the force: FH=Tcosα2-Tcosα1. This horizontal component can be neglected since, for small values of the angles α1 and α2, we have cosα1≈1 and cosα2≈1. Therefore we find that FH≈0. The upward force is Tsinα2 and the downward force is Tsinα1. The vertical force FV acting on Δx is
(108)FV=Tsinα2-Tsinα1.
We divide this equation by T and (recall that cosα1≈1 and cosα2≈1) and we get
(109)FVT=Tsinα2Tcosα2-Tsinα1Tcosα1=tanα2-tanα1.
According to Newton’s second law, the net force must equal the product of the mass of Δx and its net acceleration ∂2ϕ/∂t2. If we define the mass per unit length of the string by m, the mass of our piece is mΔx. Using the above for FV, using (109), and rearranging terms, we get
(110)1Δx(tanα2-tanα1)=mT∂2ϕ∂t2.
Because the tangents can be written as the slope of the curve ϕ adopted by the string, we can rewrite the expression left of the equal sign as
(111)1Δx(tanα2-tanα1)=1Δx(∂φ∂x|x2-∂φ∂x|x1),
and for Δx→0, this expression becomes
(112)∂2ϕ∂x2.
Using this result and substituting propagation speed v=T/m in (110), we obtain the wave equation
(113)∂2ϕ∂t2=v2∂2ϕ∂x2.
Now we add a local spring connected to our little piece of string, an extra external force and we embed our mass of the string in a damping environment (e.g., molasses); this situation is depicted in Figure 19(b). Now we have to include these forces in Newton’s law for FV in (113)
(114)m∂2ϕ∂t2=T∂2ϕ∂x2+FLocal+FDamping+FExternal.
The local force due to the spring is proportional to the position of the spring ϕ; that is FLocal=-Kϕ, K is a spring constant and there is a minus sign because a positive deflection of the piece of string causes a force in the opposite direction. The damping force is proportional to the vertical velocity of the piece of string ∂ϕ/∂t; also in this case there is a minus sign because a damping force works in a direction opposite of the velocity: FDamping=-B(∂ϕ/∂t). Thus both the local spring and damping forces are restoring forces. Putting this all into (114) yields.
(115)(Km-v2∂2∂x2)ϕ+Bm∂ϕ∂t+∂2ϕ∂t2=1mFExternal.
It can be seen that this result and the one from Jirsa and Haken [45], (75), are very similar. Equation (75) explicitly includes a nonlinear component in the force term, whereas FExternal here may be either linear or nonlinear.  The result in (115) is an expression that can explain multiple types of behaviors. If B (damping) is very large the equation starts to resemble the diffusion equation ∂ϕ/∂t=D(∂2ϕ/∂x2) whereas for very small B the equation resembles more the wave equation ∂2ϕ/∂t2=v2(∂2ϕ/∂x2) [32]. If on the other hand, the spatial derivative and the external force term dominate, the expression becomes similar to Poisson’s Equation ∇2ϕ=f.Diagram of string models used for EEG. (a) A string attached at its ends but unattached to local structures; the associated detail shows the forces acting on a small pieceΔx between x1 and x2. Forces acting on each side of this small piece of string are equal to the tension T in the string. The curve ϕ adopted by the perturbed string is associated with angles α1 and α2 at x1 and x2respectively. (b) shows a detail of the same string but now there is an external force, a local effect (a local spring), and the string is embedded in a damping medium (molasses).
(a)
(b)Finally we examine the result (75) from Jirsa and Haken [45] for a special case, a generalized seizure (e.g., Figures 17(b) and 17(c)). Since we have generalized activity we assume that there is spatial homogeneity, and since the seizure sustains itself we assume there is no external input; due to these assumptions, the partial derivative term ∂2ψe/∂x2 vanishes and so does the external input part in the ρ variable (see (70) for the definition of ρ). The Taylor series approximation of this variable is ρ(x,t)≈aeαeρ~ψe(x,t). Substituting this approximation into (75) and using the above considerations, we can simplify
(116)ψ¨e(t)+(2ω0-ω0aeαeρ~)︸bψ˙e(t)+(ω02-ω02aeαeρ~)︸cψe(t)=0.
The general solution of this ODE is in the form A1e-λ1t+A2e-λ2t, with λ1,2=-b/2±b2-4c/2. Therefore, we may expect intrinsic global oscillations only if λ1,2 are imaginary, that is if the condition b2-4c<0 is satisfied. Now note that in the simplified equation above c=ω0(b-ω0), if we substitute this into the condition above: b2-4c=b2-4bω0+4ω02=(b-2ω0)2. Since the outcome is a quadratic expression, the condition b2-4c<0 will never be satisfied and therefore, without oscillatory input (e.g., an external network or a set of pacemaker neurons in the network), oscillations cannot occur in this generalized and linearized version of the model.


## 7. Concluding Remarks

In the previous sections, I presented a variety of modeling approaches used in neuroscience. As shown in Section6, most models are capable of explaining aspects of epileptiform behavior. Interestingly, some of these explanations are counterintuitive (e.g., [96, 106]). Considering the rich repertoire of model behaviors, it is not surprising that a single approach to treatment, such as adjusting the excitation-inhibition balance, is likely to fail. A common theme in the different models is that oscillatory behavior arises from the force term in the differential equations; in the neuronal networks, these force terms may represent membrane currents of the cells that create an oscillatory input and/or the properties of the network that embeds them. A logical next step would be to create a systematic summary of hypothetical mechanisms leading to epileptiform activity and to examine them in experimental models and clinical cases.In general, the interest in and acceptance of modeling in neuroscience is growing (e.g., [46, 50]). However, an important roadblock for the application of modeling is that, in part due to the complexity of the nervous system, there is little consensus amongst neuroscientists on what properties are to be considered critical in a theoretical approach.  For this reason, traditionally, the neuroscientist prefers “realistic models” in which many experimental details are represented.  Of course, in this sense, the skeptical neuroscientist is right that all models are wrong, that is, no model fully captures reality, both by design and because some details are simply not known. To cite Rosenblueth and Wiener [107], “the best material model for a cat is another, or preferably the same cat.” It is obvious that such a “best material model” is not the most helpful tool for a theoretical study of specific aspects of cat behavior because it clearly contains too much detail. Similarly, if one models traffic flow, one probably would not include some details (e.g., color of the cars) while including others (e.g., number of traffic lanes between suburbs and industrial zones). Just like in physics, when modeling in biology, it is good to recall Einstein's dictum: “models should be as simple as possible, but not more so.” (e.g., [108]). Due to simplifications, models do not include all possible parameters, but they can nonetheless be very useful because they allow us to make predictions and they may generate uncluttered insight into the critical components of the underlying mechanisms of the modeled process. Of course, the above general statements are valid for modeling in neuroscience, whether we deal with models that describe the microscopic details of cellular function or ones that deal with the macroscopic networks of billions of nerve cells in the human brain.Because of the necessary simplifications required in modeling, there are limitations in the models reviewed here; they neglect many aspects that play important roles. A critical part of the dynamics of neuronal networks is determined by synaptic plasticity, that is, the change in network connection strength. These changes play a crucial role in the process of learning but, when pathologic, may also contribute to diseases such as epilepsy. It was more than a half century ago that Hebb [109] recognized this and postulated his law. Hebb’s law describes how synaptic coupling strength depends on local activity. It is often presented in a simplified form as “neurons that fire together wire together.” Much later, an experimental approach has confirmed the existence of activity dependent plasticity at the level of the coupling between nerve cells, that is, the synapse (e.g., [110, 111]). Here the plastic changes appeared to depend on the timing of the action potential (spike) activity of the pre- and postsynaptic neurons and is therefore also called spike-time-dependent plasticity (STDP). STDP has been used in simulations to find appropriate synaptic strength (e.g., [25]), effectively as a parameter estimation technique. It is fairly obvious that if connections between excitatory neurons that fire together also wire together (i.e., get strengthened) represents a positive feedback mechanism, which may lead to instability and possibly seizure-like hyperactivity [112]. Therefore, under normal physiological conditions, there must be additional mechanisms counteracting this positive feedback. Recently, an interesting attempt was made to analyze effects of plasticity of inhibitory synapses at the population level and successfully relate it to the network’s balance between excitation and inhibition [113]. Since this model describes strengthening of the inhibitory connection, the instability issue does not arise here. Although plasticity is an important aspect of network behavior, the models reviewed here did not include it as one of the properties. Plasticity and many other aspects involved in nervous system function, such as the role of glia (e.g., [114]), remain to be investigated.Finally, one important conclusion is that there is not “a best approach” when modeling neural function. For example, nonlinear models can help us to understand specifics such as sudden transitions, but linear ones are more suitable to analyze subthreshold properties such as oscillation. The linearized versions of the models also make it easier to see similarities between the models across different levels of organization. Another example: complex models can be studied with simulations, they are more complete since they can deal with more parameters than the more abstract models used for mathematical analysis. However, both have a place and they can be complementary. Analysis of simplified mathematical models can generate fundamental insight in spite of the fact that they lack detail. Simulation can be too complex to directly generate such fundamental insight into the interaction of the ongoing processes, but they have the advantage to be closer to reality. This property appeals to the experimenter since it may produce data with a temporal-spatial resolution that is not (yet) feasible experimentally. In this context, the ultimate goal in understanding neural networks in physiological and pathological states, is to create a framework of experimental models, detailed simulations, and mathematical formalisms that allow us to understand and to predict dynamics of network activities including state transitions; that is, results in one model can be used to inspire work in another type of model. In the case of analyzing a heterogeneous network disease such as epilepsy, modeling can provide an overview of hypothetical network mechanisms involved in ictal activity that can be employed as a guide for experimental and clinical investigation.



---

*Source: 871472-2013-03-07.xml*