Small changes in synaptic gain lead to seizure-like activity in neuronal network at criticality

Epilepsy is a neurological disorder characterised by spontaneous recurrent seizures. The mechanisms by which multiple molecular and cellular changes lead to seizures is not well understood. Here, we study cortical seizure generation by simulating the activity of neuron groups in a network using the laminar cortex model. We identified a clear boundary between low-amplitude, asynchronous activity and high-amplitude, rhythmic activity, around which small changes in excitatory synaptic gain led to strong oscillatory activity. Neuron groups only responded significantly to stimulation around the boundary. The consequences of biophysical changes induced by epilepsy-related SCN1A mutations were also examined. Marked reduction in neuronal inhibition, as caused by mutations underlying Dravet syndrome, invariably led to strong neuronal firing, whereas small reductions in inhibition could cause significant changes when the network was poised close to the boundary. The study highlights the critical role of network dynamics in seizure genesis.


Results
Architecture of the LCM. For the reader's convenience, we introduce the structure of the LCM (see Fig. 1).
The LCM represents a small region in the visual cortex of the cat as a sheet of cortical columns, within which neuron groups are situated (see Methods and Fig. 1B). Neuronal dynamics are modelled using a 'mean field approximation' approach, which treats the same type of neurons within a column as a group acting as a single entity within the network. The dynamics and connections of single neurons in a neuron group are averaged using the mean field approximation 14,16 . We expanded the original LCM to incorporate a thalamocortical loop. Hence, the model now incorporates eleven neuron types in the cortex and three neuron types in the thalamus (see Fig. 1C,D and Table S1). Interactions between neuron groups are controlled by synaptic connections within the cortex and the thalamus (see Fig. 1C and Table S5). The LCM models the stimulus input to the thalamus as point spike sources, with spike rates generated using a Monte Carlo method (see Methods). For each run of the LCM, the membrane potentials of all neuron groups are first initialised at their resting values (−65 V), and then updated every 0.5 msec through five processes: action potential (AP) generation, AP conduction and distribution, synaptic transmission, postsynaptic potential (PSP) aggregation, and membrane potential (MP) formation (see Fig. 1D). Each of these processes is modelled using a set of equations elaborated upon in the Methods. The results reported below were generated by an instantiation of the LCM configured to simulate a cortical area of 1.12 × 1.12 mm 2 , which contains 20 × 20 cortical columns with the size of 56 μm 17 . The results were produced from a 5.12-second recording of neuronal firing rates after 10 seconds of running of the model.
Temporal features of different states of neuronal activity. Parameters used to model the behaviour of neuron groups were derived from published experimental results (See Tables S1-S5). We found that the excitatory and inhibitory synaptic gains (g E and g I ), which respectively control the efficiency of excitatory and inhibitory synaptic transmission, dramatically influence neuronal firing rates. Mathematically, synaptic gain can be expressed as (also see Equation (7)), where PSP qp is the amplitude of PSPs at an afferent spike rate φ p , N qp is the number of synapses, and f(V q , φ p ) is a scaling factor for spike adaptation and membrane potential dependency (refer to Equation (7) for more details). From a network perspective, excitatory and inhibitory synaptic gain respectively reflect the excitatory and inhibitory connection strength between neurons in the network with changes in gain corresponding to the effects of short-term synaptic plasticity. Figure 2 displays the neuronal firing rates for a range of synaptic gains. Three types of neuronal activity with markedly different firing rates and oscillation states were noted. Low amplitude asynchronous activity occurred at low excitatory synaptic gain [<0.38 uV/(mV · Hz)] or high inhibitory synaptic gain (see Fig. 2). Here, the mean firing rate was less than one AP/sec. Neuronal activity was of low frequency with no specific frequency being dominant. High amplitude asynchronous activity was produced when the excitatory synaptic gain was high and the inhibitory gain was small. The mean firing rate exceeded 60 AP/sec, and in the frequency spectrum of neuronal firing rates, low-frequency oscillations were enhanced. When both excitatory and inhibitory synaptic gain were large, strong rhythmic activity occurred. In this state, the neurons were highly activated (mean firing rates >10 AP/sec), and the frequency spectra showed obvious peaks. Figure 3 displays four quantities calculated from the neuronal firing rates obtained using different synaptic gains: temporal mean firing rates (MFR), the mean power spectrum densities (PSD) of the low (2-15 Hz, PSD L ) and the high (16-50 Hz, PSD H ) frequency bands, the mean magnitude-squared coherence (MSC) of the low (MSC L ) and the high (MSC H ) frequency band. The first three measurements quantify temporal properties of neuronal firing and the MSCs measure coherence between neuron firing in columns. In Fig. 4, parameters are plotted against excitatory and inhibitory synaptic gains. Figures 3 and 4 reveal a complex relationship between neuronal firing rate and synaptic gains. The firing rate generally increases as excitatory synaptic gain (g E ) increases or inhibitory synaptic gain (g I ) decreases in magnitude. However, the change in firing rate is not continuous; a series of critical synaptic gain combinations are observed, around which small changes in excitatory or inhibitory gain cause dramatic jumps in the firing rate. For example, with an inhibitory gain of 1 uV/(mV · Hz), the firing rate jumps from about 0.3 AP/sec to >60 AP/sec as excitatory gain increases from 0.394 to 0.397 uV/(mV · Hz) (refer to Fig. 3D). On the synaptic gain maps shown in Fig. 3, the critical combinations of excitatory and inhibitory synaptic gains fall along a line: + . = . + . g g g g 0 018 0 376, or 5 5 6 20 89, where g E,C and g I,C are the critical excitatory and inhibitory gains, respectively. Secondly, this line divides neuronal firing rates into two regions: Region 1: where g E < g E,C and g I > g I,C , in which the neurons fire at low rate (mean firing rate <0.4 AP/sec), and Region 2: where g E > g E,C and g I < g I,C , in which the neurons fire at medium to high rates (mean firing rate >10 AP/sec, refer to Fig. 4A,B). Neuronal firing rates in the two regions exhibit different dependencies on synaptic gains. In Region 1, the firing rate is positively correlated with both the excitatory (Pearson's correlation coefficient r = 0.65) and the inhibitory gain (r = 0.38). In Region 2, firing rate is weakly negatively correlated with excitatory gain (r = −0.31), but strongly negatively correlated with inhibitory gain (r = −0.93). Neuronal firing rate increases exponentially as the magnitude of inhibitory gain decreases. The firing rate increases from about 20 AP/sec at g I = 6 uV/(mV · Hz) to 40 AP/sec at g I = 2 uV/(mV · Hz) and about 100 AP/ sec at g I = 0.2 uV/(mV · Hz). Excitatory gain does not affect the firing rate significantly (see Fig. 4B), but the range of firing rates increases slightly with a decrease in the lower limit (refer to Fig. 4A). In Region 1 (g E < g E,C ), oscillations in both low and high frequency bands increased with excitatory gain (r = 0.78 for low frequency band; and r = 0.87 for high frequency band; refer to Fig. 4C and E), and were not significantly affected by changes in inhibitory gain (r = −0.1 for low frequency band; and r = 0.11 for high frequency band; refer to Fig. 4D and F). In Region 2 (g E > g E,C ), the PSDs remained low [<0.1 (AP/sec) 2 /Hz] for excitatory gain <0.43 uV/(mV · Hz). When excitatory gain exceeded 0.43 uV/(mV · Hz), the PSDs varied between 0-2 (AP/ sec) 2 /Hz for the low-frequency band and 0-8 (AP/sec) 2 /Hz for the high-frequency band. As with neuronal firing rate, the PSDs were more strongly dependent on inhibitory gain than on excitatory gain (refer to Fig. 4C-F). PSDs in both low and high-frequency bands increased as inhibitory gain increased, implying that strong neuronal oscillations require high inhibitory activity. The PSDs of high but not low-frequency bands appeared to saturate at g I = 3 uV/(mV · Hz) (refer to Fig. 4F). The MSCs in both low and high-frequency bands did not vary significantly with either excitatory or inhibitory gains in Region 1 (g E < g E,C ). In Region 2 (g E > g E,C ), the correlations between the MSCs and excitatory gains were moderate in the low-frequency band (r = 0.59) and low in the high-frequency band (r = 0.27). The MSCs had a strong, non-linear relationship with the inhibitory gains in the region. Highly coherent neuronal firing occurred at high excitatory gain [g E ≥ 0.43 uV/(mV · Hz)] and moderately low inhibitory gain (between 1.3 and 3.3 uV/ (mV · Hz); refer to Fig. 4G,H).
The LCM incorporates spike inputs from two external sources: cortico-cortical projections to cortical neurons and sensory inputs to thalamic neurons. Cortico-cortical input was modelled as low-magnitude (1 AP/sec), low-frequency (frequency cut off at 20 Hz) white noise resembling temporal features of alpha activity generated in the absence of a stimulus. Afferent thalamic input was modelled for two conditions: (1) the non-stimulated state, modelled as low magnitude (1 AP/sec) low frequency (frequency cut off at 20 Hz) white noise; and (2) the stimulated state, modelled as high magnitude (50 AP/sec) high frequency (frequency cut off at 50 Hz) white noise. Cortico-cortical connections to each column were treated as being independent (i.e., asynchronous input) whereas thalamic inputs were the same for all columns (i.e., synchronous input). The neuronal firing rates and their response to stimulation are shown in Figs 2 and 3. Stimulation did not change mean firing rates significantly ,C but firing rate increased by as much as 80 AP/sec when g E was close to g E,C (see   Response to transient changes in synaptic gain. The results shown above were simulated using fixed synaptic gains. However, the efficiency of synaptic transmission may be influenced by many factors. We therefore examined the changes in neuronal firing rates induced by slight, transient changes in synaptic gain (see Fig. 5). For most combinations of synaptic gains, a change in firing rate occurred within milliseconds and neuronal activity returned to its previous state as soon as gains were restored to their initial values. The transition was more complex, however, for certain combinations of synaptic gain. For example, for the combination of g E = 0.4 uV/(mV · Hz) and g I = 2 uV/(mV · Hz) (see Fig. 5), the effects of changes in excitatory gain persisted for hundreds of milliseconds. For the combination of g E = 0.403 uV/(mV · Hz) and g I = 2 uV/(mV · Hz), a 0.015 uV/(mV · Hz) increase in excitatory gain resulted in significant changes in neuronal activity only after a 1-second delay with a strong neuronal oscillation that lasted almost 0.5 second before steady state was reached. Neuronal activity did not return to its original state after synaptic gains were restored to their initial values. Similar behaviours were also observed with temporary changes in inhibitory synaptic gain. It should be noted that the complex transitions in neuronal activity occurred when g E increased from below g E,C to above g E,C or g I decreased from above g I,C to below g I,C . Response to changes in the reversal potential of inhibitory synapses. We also examined neuronal network response to changes in the reversal potentials of inhibitory synapses. Figure 6 displays neuronal firing rates obtained using different combinations of inhibitory synaptic gains and reversal potentials of inhibitory synapses. We observed a boundary in the map, around which a small positive shift in the reversal potential led to a jump in neuronal firing rates (Fig. 6A). The critical values of the reversal potentials moved positively when inhibitory gain increased. In addition, neuronal oscillations did not change significantly with the reversal potential under weak inhibition but the oscillation amplitudes in both low-and high-frequency bands increased significantly when reversal potentials shifted positively (Fig. 6B,C).
Neuronal activity with changes associated with SCN1A mutations. The LCM can be used to analyse the behaviour of the network with changes in the biophysical properties or synaptic connections of neurons. We examined the effects of biophysical changes associated with SCN1A gene mutations. This gene encodes the alpha 1 subunit (Nav1.1) of the neuronal voltage-gated sodium channel. Mutations are associated with a wide range of epilepsy syndromes, ranging from relatively mild simple febrile seizure (FS) and Generalised Epilepsy with Febrile Seizures Plus (GEFS+) to more severe Intractable Childhood Epilepsy with Generalized Tonic-Clonic seizures (ICEGTC) and Dravet syndrome (DS) [18][19][20][21][22] . A variety of biophysical changes have been associated with different mutations, including impaired firing capability and positive shifts in the firing thresholds of inhibitory neurons [23][24][25] . Impaired firing capability was simulated by gradually reducing maximum firing rates F I [max] (see Equation (5)) from 200 to 40 AP/sec, in steps of 10 AP/sec, in all inhibitory neuron groups. Figure 7 displays the neuronal firing rates with reduced F I [max] and different synaptic gains. As expected, the neuronal firing rate increased as F I [max] decreased for all synaptic gains tested, but the firing rate did not change continuously. It jumped from a low rate to a high rate around certain values of F I [max] , and the closer g E was to g E,C , the smaller the reduction in F I [max] required for an abrupt transition (refer to Fig. 7E). PSDs at both high and low frequency increased as F I [max] was reduced and there was an abrupt transition at the same F I [max] as that at which the abrupt transition in firing rates was observed. PSDs at both high and low frequency decreased when F I [max] was extremely low (F I [max] < 70 AP/ sec).
A positive shift in the firing thresholds of inhibitory neurons was simulated by gradually increasing the voltage at half-maximum-firing (V I [HMF] ; see Equation (5)) for all inhibitory neuron groups from −45 mV to −19 mV in steps of 2 mV. Figure 8 displays the neuronal firing rates with changes in V I [HMF] and different synaptic gain combinations. Positive shifts in firing thresholds caused the neuronal firing rate to increase discontinuously. Even a small, 2-4 mV, shift in neuronal firing threshold caused a significant increase in the firing threshold (refer to Fig. 8E). The PSDs at both high and low frequencies increased after a small positive shift in firing threshold, but they returned to minimum values for V I [HMF] > −29 mV.

Discussion
In this paper, we study the neuronal activity in the cortex using computer simulations based on the LCM. The LCM incorporates many realistic features of neuronal networks allowing us to gain insights into neuronal dynamics under various conditions that are relevant to the generation of seizures.
Neuronal network dynamics. For the normal function of a brain network, outputs (neuronal firing rates) associated with different input conditions (stimuli) should be distinguishable from each other. Our results suggest that this is the case when synaptic gains are close to their critical values (i.e., g E close to g E,C ). In the cortex, thalamic projections account for less than 5% of the total synapses on a neuron (refer to Table S5 in the Supplementary Information) and are greatly outnumbered by distant cortico-cortical projections and local connections. Spike inputs from the thalamus are unlikely to evoke significant responses in the firing states of the cortical neurons unless their membrane potentials are already poised close to the firing threshold. Therefore, the neuronal network must maintain some level of activity even in the absence of the stimulus. This may help to ). explain previous observations that noise may help to improve the function of neuronal networks (see, for example, paper of Hansel and Vreeswijk 26 ). The neuronal network exhibited significantly different dynamics in the two synaptic gain regions, indicating that transitions of neuronal network phases occurred around the boundary between the two regions. The network in Region 1 (g E < g E,C ) is in a disordered phase, because the activity of neuron groups is more strongly influenced by noise than by neighbouring neurons. The network in Region 2 (g E > g E,C ) operates in an ordered phase because the activity of neuron groups is highly correlated. Our results suggest that the neuronal network of the brain operates in a balanced critical state between the two phases. This finding aligns with the proposal that neuronal networks in the brain operate in a self-organised critical state 27,28 . The criticality hypothesis has gained support from both theoretical considerations and experimental observations. Theoretically, a neuronal network in a critical state can maintain some order to ensure consistent responses to specific stimuli but allows a degree of disorder to enable adaptation to variations in external conditions 29 . Experimentally, neuronal avalanches with an inverse power-law distribution for avalanche size and duration, a characteristic of systems in a self-organised critical state, have been observed in the spontaneous activity of organotypic neuronal cultures from rats 28,30 and humans 31 . In diseased tissue, such as that from patients with epilepsy, neuronal avalanches deviate from the power-law distribution 32 . Our simulations provide further compelling support for the criticality hypothesis insofar as neuronal networks only show a significant response to stimuli when they are situated close to the critical state (refer to Fig. 3).
Many interesting behaviours were observed when the network was close to the critical state. Firstly, the neuronal network in a critical state displayed some degree of adaptation to changes in neuronal excitation and inhibition (refer to Fig. 4). The network first reacted to increased neuronal excitation with strong oscillatory activity. Then, it gradually absorbed the changes by increasing firing, and oscillation amplitude fell significantly. The observed adaptation may help to clarify a long-standing conundrum: how does a neuronal network at criticality avoid slipping into either a completely disordered or a completely ordered state in a noisy environment. Our results suggest that the neuronal network itself is capable of absorbing temporal changes in neuronal excitation and inhibition. This, together with adaptation in synaptic transmission and the refractory period of neuron firing, may confer tolerance in the network to significant fluctuations in background activity. Secondly, the network showed hysteresis -network responses were delayed (by about one second) after the change in synaptic gain (refer to Fig. 4). The hysteresis phenomenon is unlikely to be caused by AP propagation along axons or conduction of PSPs along dendrites, because the former only lasts several milliseconds and the latter less than 100 milliseconds. A more likely cause for the delay is the cascade of neuronal excitability within the network (see below). This suggests that the network may be able to accumulate neuronal excitation for a period much longer than the duration of an AP or PSP, which may have important implications for the generation of seizures.
Implications for seizure generation. Seizure initiation in the brain is a temporal dynamical process resulting from changes in neurons, synapses, and ion concentrations. Simulation of different states of the neuronal network predicts the changes in activity that are associated with changes in state. This can help to inform our understanding of how cellular and molecular changes can evoke or inhibit seizure generation.
Although it is possible for the neuronal network to achieve a critical state over a range of inhibitory activity, seizure vulnerability is greatest at a low level of inhibition. At a small inhibitory gain (<1 uV/(mV · Hz)), the rate of neuronal firing is close to maximal when the excitatory gain exceeds the critical value (see Fig. 3). The high firing rates persist even after excitatory gain falls below the critical value (further discussed below). Furthermore, stimulation significantly increased synchronisation when inhibitory gain was small. Seizure generation in neuronal networks with weak inhibition may result from the amplification of synchronised inputs in noisy environments.
The simulations revealed distinctive roles of neuronal excitation and inhibition in regulating network dynamics. Excitation imposed a threshold on neuron firing with firing rates being lowest when neuronal excitation was below threshold and high firing rates when excitation exceeded the threshold. Above the threshold, inhibition plays an important role in controlling neuronal firing. The results can be explained by the interactions between three factors that determine neuronal membrane potential: leaking currents through the membrane, depolarisation caused by excitatory PSPs, and repolarisation caused by inhibitory PSPs (see Fig. 9A). Leaking currents reduce the membrane potential according to τ ⋅ − u t exp( / ) m , where u is the membrane potential above the resting potential, and τ m is the membrane time constant (see, Equation (3) in Methods) 6 . The amplitudes of PSPs are influenced by neuron membrane potentials through two mechanisms. Firstly, the efficiency of excitatory synapses decreases and that of inhibitory synaptic transmission increases with the membrane potential of postsynaptic neurons. The relationships were modelled as linear functions in the LCM (see Equation (7) and Fig. 9E,F). Secondly, the afferent spike rates to synapses increase with the membrane potentials of presynaptic neurons (see Fig. 9C,D). The membrane potentials of pre-and postsynaptic neurons are highly correlated in the network due to the strong recurrent synaptic connections between neurons. The three factors play different roles in regulating neuron firing at different levels of excitatory gain. When the excitatory gain is low, PSP amplitudes do not change significantly because the afferent spike rates to synapses are low. However, leaking currents, which increase linearly with the membrane potential (see Fig. 9B), strongly restrain the growth of the membrane potential. As the excitatory gain increases, the amplitude of excitatory PSPs may outweigh the repolarisation caused by leaking currents, and membrane potentials of neurons begin to grow. When they exceed the firing thresholds, neuronal firing rates increase explosively, which leads to further depolarisation. The cascade of neuronal excitability occurs as long as the excitatory PSP amplitudes are higher than the membrane potential reduction due to leaking currents, and even a small net neuronal depolarisation eventually leads to explosive neuron firing. The cascade stops when the amplitudes of inhibitory PSPs match that of excitatory PSPs because increases in membrane potential reduce the efficiency of excitatory synapses and boost the efficiency of inhibitory synapses. The larger the inhibitory gain, the sooner the cascade stops. Therefore, the final membrane potentials depend more strongly on the inhibitory gain than on the excitatory gain.
Epilepsies associated with SCN1A mutations exhibit extraordinary clinical heterogeneity. Seizure manifestations vary between members of the same pedigree bearing identical mutations 33,34 . Previous studies have ascribed the phenotypic heterogeneity to the diversity of functional changes in the voltage-gated sodium channel caused by individual mutations 35,36 . Our simulations, however, suggest the effects of SCN1A mutations on the neuronal network depend on the levels of neuronal excitation and inhibition. Though severe loss of neuronal inhibition, as caused by the SCN1A mutation associated with DS, leads inevitably to high neuron firing, a small reduction in neuronal inhibition, as caused by the SCN1A mutation associated with the GEFS+ syndrome, only causes significant changes in neuronal activity when the network is poised close to the boundary. The mutation also causes other neuronal changes, such as changes in GABA A chloride reversal potential, which may contribute to seizure generation 37 .
, where g p is the synaptic gain, V q is the membrane potential of the postsynaptic neuron,

Methods
In this section, we outline the structure of the LCM model, which enables cortical neuronal dynamics to be modelled under various conditions using model parameters that represent different neuronal behaviours. Additional details about the simulation methods and parameter values used in the paper are provided in the Supplementary Information.
Laminar cortex model. The LCM was developed to simulate local field potentials (LFP) in the visual cortex of the cat. The LCM simulates neuronal activity of a flat sheet of cortical columns 14,15 . A cortical column may contain thousands of neurons, however, neurons in a column display similar responses to specific stimuli. Therefore, the LCM treats the same type of neurons within a column as a group which acts as a single entity in a network.
A neuron group has similar features to a single neuron but its dynamics and connections are averaged using the mean-field approximation 14,16 . The approach allows us to model the neuronal activity of a large cortical region without knowing the detailed physiology of individual neurons.
In the current LCM, the evolution of neuronal membrane potentials is modelled using the following differential equation 38 is the membrane potential relative to the resting state, and V and V [0] are the membrane potential and resting membrane potential of neuron groups, respectively; C m is the total capacitance of the membrane; G m is total conductance of the membrane; the term −G m u represents the currents leaking through the membrane; and I E is the total current flowing into the neuron through synaptic transmission, i.e., postsynaptic currents (PSC). The LCM calculates membrane potentials at discrete time points. For a small time interval, a simple iteration equation can be derived from Equation (3): where Δt is the time interval; τ m = C m /G m is the membrane time constant characterising the time for the membrane to repolarize after depolarisation 6 is the membrane potential change caused by the injected currents during the period. The LCM uses equation (4) to calculate the membrane potentials of neuron groups iteratively. The membrane potential change ΔV is assumed to be the aggregate postsynaptic potential (PSP) at the soma, and is determined by three processes: neuronal firing, synaptic transmission and postsynaptic potential (PSP) aggregation. Each of the three processes is modelled using a set of equations, which are described in details below.
The LCM uses a set of equations to model neuronal processes. For notational convenience, the following conventions are adopted for the equations: (1) the LCM treats the mean (somatic) membrane potentials of neuron groups (denoted by V) as state variables; (2) the LCM uses spike rates (denoted by φ or F) to measure the interaction between neuron groups. The spike rate is defined as the average number of spikes a neuron of one group receives from another group per unit time; (3) neuron groups can be differentiated by neuron type (as listed in Table S1 in the Supplementary Information) or column number. Subscripts q or p may be used to denote a neuron group in a column. For synaptic transmission, p and q respectively indicate the presynaptic and postsynaptic neuron group; (4) the LCM simulates PSPs (denoted by PSP) generated through effects on different neurotransmitter receptors. The superscript [rcpt] denotes the receptor type; (5) Subscripts and superscripts may be omitted when the equation applies to all relevant objects in the same class (such as neuron groups, receptors, or synaptic connections) equally but parameter values may differ.
The LCM determines the mean spike rate generated by a neuron group using a sigmoid function 16,38 : [max] [HMF] where F(t) is the mean spike rate, k is the firing gain, F [max] is the maximum firing rate of the neuron groups, and V [HMF] is the membrane potential at half maximum firing rate. Equation (5) is similar in behaviour to a power-law function at low voltage 39 and converges to F [max] at high voltage. The LCM also incorporates time delays of spike propagation. The time delay for a spike to reach a postsynaptic neuron group is calculated with: [AP] [axon] [AP] where s [axon] is the length of the axon connecting the two neuron groups, and v [AP] is the spike propagation speed, which is set to 1.0 m/sec 6 . In the LCM, axon length is assumed to be the distance between the two groups (i.e. direct synaptic connections between neuron groups), and spike rates remain constant during propagation. Afferent spikes change the membrane potentials of the target neuron groups. The amplitude of the postsynaptic potential (PSP) depends on the afferent spike rate and the membrane potential of the target neuron group. PSP amplitude is determined by 16 : where the subscript 0 denotes the PSP amplitude at the synapse (see below), φ p is the afferent spike rate (see below), g p is the synaptic gain at the resting membrane potential (referred to as synaptic gain), N qp is the synapse number between the two neuron groups, V p [rev] is the reversal membrane potential (0 for excitatory spikes and −70 mV for inhibitory spikes) and λ [rcpt] is the spike adaptation parameter. Three kinds of neurotransmitter receptors were considered: AMPA, NMDA, and GABA receptors. The excitatory synaptic gain (g E ) controls the PSPs of both AMPA and NMDA receptors, and inhibitory synaptic gain (g I ) controls PSP of the GABA receptors. The time course of the PSP is modelled with modified bi-exponential functions: [rcpt] 2 is the normalising constant, and H(x) is a Heaviside step function, with a value of 1 when x ≥ 0 and 0 elsewhere. The LCM also considers the time delay for PSPs to reach the neuron body using an approach similar to Equation (6). It also assumes that PSP amplitudes decay exponentially during propagation, i.e.: where PSP qp,0 is the PSP amplitude at the synapse, PSP qp is the corresponding membrane change in neuron body, s [dend] is the distance between the synapse and the neuron body, and λ [PSP] is the PSP decay factor. The change in membrane potential in the target neuron group is the aggregate of all PSPs 16 : where ⊗ denotes convolution in time, and aggregation is performed over all receptors and presynaptic neuron groups in all columns. The model incorporates the laminar architecture of the cortex represented as five laminae (layers I to VI with layers II and III combined). The cortex is divided horizontally into columns containing eleven neuron groups spanning the five cortical layers (see Fig. 1). The model uses a quantitative synaptic connection map to define the synaptic connections between neuron groups (see Table S5 in Supplementary Information). This map was derived from the literature 40 . Each entry in the map represents the average number of synapses on a neuron that project from a specific type of neuron. A neuron may receive thousands of synaptic connections from many surrounding cortical columns. The synaptic connection map only shows the types of presynaptic neurons, and does not contain information on columnar origin. The LCM also assumes that the density of presynaptic neurons connected to a specific neuron is normally distributed according to the following equation: where r(s qp ) is the proportion of synapses from neurons located at a distance s qp from the target neuron; σ [synp] is the standard deviation of the spatial distribution of the presynaptic neurons, set to 80 um for pyramidal neurons 41 and 40 um for spiny stellate and inhibitory neurons. X is a Gaussian random number with mean of 1 and standard deviation of 0.2, to avoid unrealistic synaptic connections that are totally symmetric.
The LCM also incorporates the thalamocortical loop. The thalamus is modelled as an additional layer to the previously described cortical model, allowing the propagation delay between the thalamus and cortex to be simulated. Three groups of thalamic neurons are considered: interneurons in the reticular nucleus of the thalamus (IRTN), relay neurons in the lateral geniculate nucleus (RLGN), and interneurons in the lateral geniculate nucleus (ILGN). Synaptic connections formed between thalamic neuron groups and cortical neuron groups are added. The source code of the LCM may be downloaded from the GitHub website (http://github.com/jiaxin-du). simulation parameters. A large number of parameters is used to describe the behaviour of neuron groups and their connections. These are divided into four categories: (1) parameters that describe the physiological properties of neuron groups (see Table S3 in the Supplementary Information); (2) parameters that describe the properties of synapses and receptors (see Table S4 in the Supplementary Information); (3) parameters that describe cortical structure (see Table S2 in the Supplementary Information), and (4) parameters that describe synaptic connections between neuron groups (see Table S5 in the Supplementary Information). Most parameter values are based on published experimental data. We note that the physiological parameters are not independent and that changes in one parameter can usually be compensated by changes in other parameters (refer to our previous paper 14 ).

Data analysis.
To determine the excitation level of a cortical column, we calculate the mean firing rates of excitatory neuron groups in the column, from which the MFRs, PSDs and MSCs of neuronal firing are determined. The MFRs and PSDs were estimated from the firing rates in the central column, and the MSCs were estimated from these in the central 100 columns. The MFRs were the temporal mean of firing rates weighted by a tapered cosine window function, which was used to minimise the cut-off effects of oscillatory firing rates. The PSDs and MSCs of firing rates were estimated using Welch's method 42 . The firing rates were first split into ten segments of equal length with 50% overlap between consecutive segments. Each segment was then multiplied by a Hamming windowing function. PSDs and MSCs were determined from the windowed segments using a fast where PSD(F i ) computes the PSDs of firing rates in column i (F i ); CPSD(F i , F j ) computes the cross power spectrum densities between F i and F j ; and N = 100 is the number of columns used for MSC calculation. The final PSDs and MSCs were the mean PSDs and MSCs of the segments, respectively.