# Metabifurcation analysis of a mean field model of the cortex

## Abstract

Mean field models (MFMs) of cortical tissue incorporate salient, average features of neural masses in order to model activity at the population level, thereby linking microscopic physiology to macroscopic observations, e.g., with the electroencephalogram (EEG). One of the common aspects of MFM descriptions is the presence of a high dimensional parameter space capturing neurobiological attributes deemed relevant to the brain dynamics of interest.

We study the physiological parameter space of a MFM of electrocortical activity and discover robust correlations between physiological attributes of the model cortex and its dynamical features. These correlations are revealed by the study of bifurcation plots, which show that the model responses to changes in inhibition belong to two archetypal categories or “families”. After investigating and characterizing them in depth, we discuss their essential differences in terms of four important aspects: power responses with respect to the modeled action of anesthetics, reaction to exogenous stimuli such as thalamic input, distribution of model parameters and oscillatory repertoires when inhibition is enhanced. Furthermore, while the complexity of sustained periodic orbits differs significantly between families, we are able to show how metamorphoses between the families can be brought about by exogenous stimuli.

We here unveil links between measurable physiological attributes of the brain and dynamical patterns that are not accessible by linear methods. They instead emerge when the nonlinear structure of parameter space is partitioned according to bifurcation responses. This partitioning cannot be achieved by the investigation of only a small number of parameter sets and is instead the result of an automated bifurcation analysis of a representative sample of 73,454 physiologically admissible parameter sets. Our approach generalizes straightforwardly and is well suited to probing the dynamics of other models with large and complex parameter spaces.

###### keywords:

bifurcation analysis, brain dynamics, mean field model, electroencephalogram (EEG), anesthesia, thalamus[cor1]Corresponding author

## 1 Introduction

Understanding brain function is one of the longstanding and unresolved quests of science. Despite the efforts and contributions of a wide range of scientists from different disciplines, we are still very far from a comprehensive account of what the brain does, and how and why it does it. The formidable hardness of this problem is embodied in the complexity of the structure and organization of this organ, which supports interactions over a wide range of temporal and spatial scales. One of the scientific approaches that aims to shed some light on the complexity of brain activity is the use of realistic mathematical models and computer simulations: constrained by known physiology and anatomy, these theories try to recreate and understand the patterns emerging in the electrical activity of the cortex. In particular, so-called mean field models (MFMs) ref:JH96; ref:RRW97; ref:LCD02 try to capture the results of neuronal interactions through spatial and/or temporal averaging, which then describe the population activity of regions of cortex in a parsimonious manner. Accompanied by a sufficiently flexible parametrization, MFMs take into account existing uncertainties in cortical anatomy and physiology in an effective way.

Despite the apparent biological simplicity of MFMs, such theories produce physiologically and behaviorally plausible dynamics. Since their birth forty years ago, they have been employed successfully to investigate a large number of cognitively and clinically relevant phenomena, such as the action pharmacological agents have on the activity of cortical tissue ref:WC72; ref:WC73; ref:LHSZ73; ref:Fre75; ref:LvRB+76; ref:ZKM78, the origin of evoked potentials ref:JR95; ref:RRW02, the sleep-wake cycle and circadian rhythm ref:SSS+05; ref:PR07, the emergence of gamma band oscillations and its relationship to cognition ref:Wri97; ref:RWR00; ref:BL07 and the onset and characteristics of epileptic seizures ref:WBBC00; ref:RRR02; ref:LBK+03; ref:KKS05; ref:LB05; c.f. Ref. ref:DJR+08 for a recent review. We use here Liley’s MFM ref:LCD02, which is capable of reproducing the main spectral features of spontaneous (i.e., not stimulus-locked) EEG, in particular the ubiquitous alpha rhythm ref:LAWA99. Firstly, this model is biologically constituted: all state variables and parameters can be constrained on the basis of existing anatomical and physiological measurements in the literature. Secondly, it supports a rich repertoire of behaviors both physiologically relevant and dynamically interesting. For example, parametrically widespread, robust chaotic activity of various origins has been found ref:DLC01; ref:vVL06; ref:FvVB+08, and multistability, i.e., the presence of various coexisting dynamical regimes, has been demonstrated. Multistability has been speculated to correspond neurobiologically to the formation of memories ref:DFCL09.

In the spirit of the dynamical approach ref:Fre01, in this paper we link nonlinear electrical activity and neurobiologically significant attributes of cortex. To this end, we consider a representative sample of parameter sets that have previously been found to generate physiologically plausible behavior ref:BL05. This sample contains 73,454 sets and for each of these we computed the bifurcation plots when varying two parameters related to inhibition, i.e., the qualitative changes of recurrent activity patterns, when neuronal inhibition is altered. It turns out that we can sort the sets into two distinct “families” of dynamical behavior. These families are found to correlate with EEG signal power and responses to anesthetics, whereas family membership is determined by specific neurobiological parameters.

The paper is divided into three main sections. First, we briefly introduce Liley’s ordinary differential Equations (ODEs). We next discuss the theoretical and numerical tools employed in some detail, in particular how one arrives at a systematic bifurcation analysis procedure to show the qualitative changes of the model solutions when inhibition is varied. Then, we present the main distinctive features of the families, such as their responses to the simulated induction of anesthetics and their correlations with model parameters of interest. We are also able to show how exogenous agents, i.e., input from the thalamus to the cortex, can induce dramatic changes in those patterns and stimulate transitions from one type of family to the other. That, in particular, provides a compelling example for the modulation thalamus is thought to exert on the cortex ref:SG06. All these relations cannot be discovered with standard linear or nonlinear analyses of the physiological parameter space, and represent the main result of this work. A discussion of open problems concludes the paper.

## 2 Neuronal mean field equations

Liley’s MFM aims to provide a mathematically and physiologically parsimonious description of average neuronal activity in a human cortex, with spatially coarse-grained but temporally precise dynamics. One excitatory and one inhibitory neuronal population, respectively, is considered per macrocolumn, which is a barrel-shaped region of approximately diameter comprising the whole thickness of cortex (thus deep). Cortical activity is locally described by the mean soma membrane potentials of the excitatory () and inhibitory () neuronal populations, along with four mean synaptic inputs , , , and . These inputs convey the reciprocal interaction between neuronal populations, where double subscripts indicate first source then target (each either excitatory or inhibitory ). The connection with measurements is through , which is linearly related to the EEG signal ref:Nun81. Lumped neuron populations are modeled as passive compartments, into which all synaptically induced ionic currents terminate. According to population types , synaptic activity drives the mean soma membrane potentials from their resting values. The equations for and are given by

(1) | |||||

(2) |

where and are mean resting potentials, and and are the membrane time constants of the respective neuronal populations. The reversal potentials of the transmembrane ionic fluxes mediating excitation and inhibition are given by and , respectively. Note that the synaptic inputs are weighted with (excitatory ) and (inhibitory ) at the resting potential of the respective excitatory or inhibitory neuronal population, and that these weights then vary linearly with voltage.

The mean synaptic inputs describe the postsynaptic activation of ionotropic neurotransmitter receptors by presynaptic action potentials, arising from the collective activity of neurons both nearby and distant. The time course of such activity, based on well-established experimental data ref:Tuc88, is modeled by a critically damped oscillator driven by the mean rate of incoming excitatory or inhibitory axonal pulses. We thus have, for :

(3) | |||||

(4) |

where the terms in curly brackets correspond to sources of the axonal pulses from three origins: local, i.e., in the same macrocolumn of the cortex , arriving through long-range, excitatory cortico-cortical connections from other macrocolumns , and extracortical, i.e., primarily of thalamic origin . For subsequent simplicity we assume the absence of any extracortical inhibitory input, i.e. . quantifies the strength of anatomical population connectivity. The maximal postsynaptic potential (PSP) amplitude occurs in the target population at time after the arrival of the presynaptic spike from the source population . A schematic illustration of the architecture of interactions in the Liley model can be found in Fig. 1.

Local mean soma potentials are nonlinearly transformed to mean neuronal population firing rates with a sigmoidal function

(5) |

where and indicate the firing thresholds and their associated standard deviations for the respective neural population. The chosen form is a computationally convenient approximation to the error function, which results when one assumes a Gaussian distribution of firing thresholds.

The final part of the model concerns the propagation of excitatory axonal activity in long-range fibers. In the full formulation, the following inhomogeneous, two-dimensional damped wave equation is used:

(6) |

Two core assumptions form the basis of this equation: first, action potentials traverse much larger distances in the cortico-cortical fibres than the local (intra-columnar) axonal fibres, and thus their conduction delays can no longer be considered negligible. Second, cortical areas further away share fewer long range connections. Given the average number of long range excitatory connections onto a local neuron , experiments suggest that the fraction thereof that stems from distant neurons decays approximately exponentially with distance, with a scale constant . Activity is then modeled as spreading omnidirectionally with constant conduction speed . We have also made the additional assumption that long-distance propagation dynamics is not differentiated according to excitatory and inhibitory targets, i.e., and .

In the following, we force the term and thus consider only bulk dynamics for our model cortex. Furthermore, we neglect the spatially inhomogeneous effects of long-range fibre systems. Thus, Eq. (6) now becomes

(7) |

with . We see that Eq. (7) has turned into an inhomogeneous critically damped oscillator. This is in contrast to simulating an unconnected local population only, which would involve dropping Eq. (6) altogether. However, it is clear that the steady state solution to Eq. (7), , simply changes in Eq. (3). Thus, the steady state of bulk oscillations corresponds to local activity with increased connectivity.

The fourteen ODEs of Eqs. (2)-(4) and (7) are studied in the present paper. We note for the interested reader that a thorough discussion of realistic axonal velocity distributions in MFMs has been recently published by Bojak and Liley ref:BL10 and that Liley et al. ref:LCD02 contains more in-depth discussions of the physiological and mathematical arguments that lead to Eqs. (2)-(6). Also, parameters used in the equations and their physiological ranges are shown in Tab. 1.

Parameter | Definition | Minimum | Maximum | Units |
---|---|---|---|---|

resting membrane potential | mV | |||

passive membrane decay time | ms | |||

excitatory reversal potential | mV | |||

inhibitory reversal potential | mV | |||

EPSP peak amplitude | mV | |||

IPSP peak amplitude | mV | |||

EPSP characteristic rate constant | ||||

IPSP characteristic rate constant | ||||

no. of cortico-cortical synapses, target excitatory | – | |||

no. of cortico-cortical synapses, target inhibitory | – | |||

no. of excitatory intracortical synapses | – | |||

no. of inhibitory intracortical synapses | – | |||

axonal conduction velocity | ||||

decay scale of cortico-cortical connectivity | cm | |||

maximum firing rate | ||||

firing threshold | mV | |||

standard deviation of firing threshold | mV | |||

extracortical synaptic input rate |

## 3 Methods

### Parameter sets and bifurcation parameters

A large number of parameter sets is considered in this study, representing a significant sample of the physiologically admissible parameter space of Liley’s MFM. To generate such realistic parameter sets, the properties of the model have been studied around fixed points of linearized version of Eqs. (2)-(6). Biologically plausible selection criteria were then invoked concerning the properties of the simulated EEG and the firing rates of their associated neuronal populations. Limits on the admissible ranges of parameter values were also applied, so that they are consistent with well-established neurophysiological features of mammalian cortex. Details on the methods and the constraints applied in this search can be found in Bojak and Liley ref:BL05. As a result, a collection of 73,454 physiologically meaningful parameter sets is available and considered in this study as representative of physiologically realistic parametric instantiations of model cortical activity in Liley’s MFM. Bifurcation theory is an essential tool for discussing the qualitative changes that solutions of dynamical systems undergo under variations of their parameters. Loosely speaking, this type of analysis provides a picture similar to a phase diagram for physical systems, so that changes in the properties of the solutions (i.e., bifurcations) may be considered as a generalization of phase transitions ref:SS10. This allows us to explore parametric boundaries between qualitatively different activities.

One of the a priori strongest motivations for the use of this approach is represented by the failure of standard methods of statistical linear analysis, which are unable to supply relevant information about the model’s 29-dimensional parameter space. An example of this is illustrated in Fig. 2. As a result of a Principal Component Analysis (PCA) over the parameters of the 73,454 sets, the percentage of total variance explained by the first ten PCs is plotted, individually and cumulatively. Essentially, if the parameter space presented strongly linear correlations among the variables, it would be possible to identify a limited number of components capturing the vast majority of the variance, indicating that the significant degrees of freedom in the space were in fact much less than its dimension. This is clearly not the case, since the first ten components cumulatively account for only about 42% of the total variance, with low and very similar singular contributions. It is not possible to find a linear combination of parameters that expresses the “bulk” properties of the parameter space. If we want to extract valuable statistical information, a different, non standard method of analysis is needed. On this basis, the aim of this paper is to highlight correlates between dynamical “snapshots” of the model, i.e. two parameters bifurcation plots, and spectral and parametric properties of Liley’s MFM.

In choosing the parameters to vary for our bifurcation analysis, we are guided by the existence of literature pointing to inhibition being a sensitive locus of control of brain dynamical responses. For instance, it is well-known that pharmacological agents alter cortical activity primarily through acting on inhibition ref:RA04. Factors that modify inhibitory PSP amplitudes and anatomical coupling strengths of the inhibitory neural population are thus selected as continuation parameters. This is in line with previous work ref:BL07 and implies changes in Equation 4 for the inhibitory PSP amplitudes

(8) |

and local inhibitory-inhibitory connectivity

(9) |

where and are chosen as our bifurcation parameters.

This choice is not exclusive and other parameters or factors could have reasonably been varied, possibly highlighting other poignant features that have not emerged through continuations in and . Nonetheless, these two parameters allow biologically plausible control: changes the amplitude of the inhibitory PSPs for both the excitatory and inhibitory populations, which for example general anesthetics and sedative agents can affect to a large degree. At the same time inhibitory-inhibitory connection strengths are altered via , which changes the impact on the inhibitory feedback loop largely responsible for generating specific oscillations ref:LCD02. Note that since in Eq. 4, the effect of inhibition on the excitatory population scales precisely with , but on the inhibitory population with . Thus, for example, setting would simulate an exclusive impact on the excitatory population. Essentially then, we can control overall inhibitory “potency” with and “specificity” with .

### One and two parameters bifurcation analysis procedure

There are two types of bifurcation plots that will be discussed: one parameter (1par) plots in and two parameters (2par) graphs in and . Let us first explain how the diagrams are obtained: we integrate Eqs. (2)-(4) and (7) for a chosen parameter set up to an equilibrium solution for and , i.e., in the absence of changes to inhibition. Initial conditions are given at by , while all the other state variables are set to zero. Once an equilibrium solution is found, it is continued in and a 1par plot is obtained. In all the -plots there are always two types of codimension one (cod1) bifurcations present: saddle-node and Hopf points, respectively. The presence of Hopf bifurcations gives us some information about the character of sustained oscillations produced by the action of inhibition, which we will discuss below.

The saddle-node and Hopf bifurcations are subsequently continued in parameters and . In our diagrams, lines describing the loci of saddle-nodes are labeled sn, and similarly for the Hopf points hb.

Physiological limits for minimal and maximal parameters intervals are set at , and . Higher or lower values of are considered to alter inhibitory PSP amplitudes and inhibitory-inhibitory connectivity strengths either to unphysiological values or to unusually large or small figures, which, for example, do not capture the commonly observed effects of pharmacological agents on the cortex. However, as discussed below, continuations are performed till particular features of the bifurcation topology become apparent. Typically, this occurs far outside of the physiological range. Bifurcation analysis is performed here with the widely used continuation software packages AUTO-07P and MATCONT ref:Doe81; ref:DGK03.

## 4 Results

### Two parameters bifurcation patterns

To establish our findings, we proceeded with a two-step analysis. First, 2par bifurcation plots for a reduced sample of sets were obtained and analyzed, visually and numerically. Analysis of plots in for 405 randomly selected sets out of the 73,454 in the generated sample revealed the existence of two unique, and clearly distinct, categories of patterns which recurrently appear. Our primary aim in this subsection is to describe these archetypal patterns and discuss their shared and distinct features. In the next subsection, a third parameter () is varied and this shows that the two families are “transformable”: it is possible to metamorphose one family into the other, by varying exogenous stimuli. Once the two families are characterized, generalization through the whole batch is straightforward and global features of the physiological parameter space naturally emerge.

The general character and topology of the 2par plots for Liley’s MFM is reproduced in Fig. 3. These represent two archetypal -plots for so-labelled Family 1 (F1) and Family 2 (F2) types of sets. The major difference among the two families is represented by how the sn lines are organized and divide the parameter space in areas with different dynamical properties. In general, because of the form of the MFM’s equations, a region bound by sn lines has always three equilibria, whereas the region outside has only one. For F1, see Fig. 3A, sn lines travel almost parallel for the whole range of the space, and do not form cusp points (c), dividing the space essentially in three major regions: one inside sn lines and two outside. In contrast, F2, see Fig. 3B, is characterized by the presence of two cusp points and , so that the region containing three equilibria is the union of two separated “wedge-shaped” areas, with cusps as their vertices. The way hb and sn lines interact is, for the majority of cases, by so-called fold-Hopf (fh) points, which appear as tangencies between such lines. Other cases we saw but are not present in Fig. 3 happen when hb lines connect two Bogdanov-Takens (bt) points. These instead correspond in the diagram to points where a hb line terminates on a sn line. In general, these bt points occur on the left, lower wedge of the sn lines, close to fh for F1 and F2. We will give examples of these bifurcations in the next subsection.

For both families, there are also other cases that slightly differ from Figure 3, although the structure of the sn lines never changes. Variations are minimal and local, i.e. one of the fh points can be missing or a small difference in the shape of the hb lines can be present, but the essential structure of the plot is not altered. What instead shows high variation, independently of the family, is the extension of the plots, i.e., the values of and at which bifurcation points appear. It is important to realize that the physiological region is typically a tiny area of the plots, as indicated by the rectangles in Fig. 3. There are plots, as those shown in Fig. 3, where the whole pattern is contained within orders of magnitude of the bifurcation parameters, but there are also a large number of patterns for which hb and sn lines extend up to immense values, like . This is an interesting, important point: scales vary widely, but the nature and structure of the patterns are invariant, and this is appreciated only if plots are considered in their complete extension in -space. Moreover, it appears that similar bifurcations maintain the same type of normal forms among equivalent patterns, no matter what scale is involved. This implies that the dynamics in the proximity of bifurcations is topologically equivalent, i.e. qualitatively very similar, among plots in the same family. Thus, parameter sets in the same family share not only the topology of the saddle-node bifurcation lines, but also yield qualitatively similar dynamics close to the codimension two points.

### Thalamic input and pattern family metamorphosis

The magnitude of extracortical (thalamic) input to inhibitory neurons is chosen as a third continuation parameter for investigating the transition between families. Our choice of has an intuitive meaning: unlike and , whose variation represent endogenously driven effects on the cortex, changes in the thalamic input capture the action of exogenous causes. Furthermore, excites the inhibitory population and hence can be expected to affect significantly the dynamical repertoire of families based on variations of inhibitory PSPs.

Let us discuss a prototypical transition from F1 to F2. To follow our description the reader will need to refer to Fig. 4 . We start at a very large and unphysiological value of , indicated in Fig. 4A, where a representative of F1 is present. Compared to the (standard) diagram depicted in Fig. 3, two fh points on the left wedge of line are not yet formed. When the input is decreased to in Fig. 4B, thus approaching the physiologically plausible range of , the hb branch extends towards lower values of , creating the point, as found in the standard F1 type in Fig. 3. Analogously, the point is created at larger values of and . At the same time, on the right wedge , two cusps and appear in Fig. 4B. This happens in a so-called swallow tail bifurcation (see, e.g., Ref. ref:arnold). All equilibria belong to the same branch, meaning that they are connected at the saddle-node bifurcations, at which the branch is folded. At this swallow tail bifurcation, a pair of saddle-node bifurcations appears on the unique branch of equilibria. Notice in the inset that the point is now on the right wedge of the little triangle formed by the two cusps and . Moreover, the single branch hb has lost its previous elliptical shape, and has produced intersections. However, these are not bifurcation points since only one unique branch of hb lines is present, and the diagram still represents a type of plot belonging to F1. The emergence of the two cusps is not considered a major change, since the topological characteristic that we choose in discriminating F1 from F2 is the global shape of the sn lines. When separate and do not intersect but keep almost parallel, as in Fig. 4A and B, a diagram of F1 type is still present.

As the extracortical input is further decreased to in Fig. 4C, the triangle of line with the two cusps and expands and moves to the left of the diagram, intersecting the opposite branch (see also the inset in Fig. 4C). Now is overlapping the line of saddle-nodes , and the cusp has crossed and is above it. Notice, in particular, the vicinity of the blue and green lines of sn in the inset: this represents a stage very close to the transition of F1 into F2. In fact, at smaller values, the cusp further lowers and separation between the two branches of sn further diminishes, so that increasing portions of and become nearly tangent. This condition persists up to a value of where sn lines finally exchange, so that at the two cusps no longer belong to the same branch, as Fig. 4D and its inset show. The cusp is now part of and belongs to : a second swallow tail bifurcation has taken place, signaling the appearance of an F2 type of diagram. In other words, two regions of the -space having triangular shapes have formed and overlapped, and now tends to move towards the right as is decreased, separating further from . A further decrease of the thalamic input will eventually result in the same plot as in Fig. 3B, i.e., the prototypical F2 with two disjoint V and -shaped curves of saddle-nodes.

A physiologically relevant observation can be made at this point. One can investigate the effect of varying the excitatory thalamic input to either inhibitory () or excitatory () cortical neurons. It turns out that there is a relatively robust family-specific effect of thalamic excitation, which is opposed for inhibitory and excitatory cortical targets: F1 (F2) is the preferred family at high (low) values of , whereas F1 (F2) is the preferred family at low (high) values of . The incidence of family metamorphosis for variations within the physiological interval of thalamic input is quite strong: about 56% of the 73,454 sets undergo a transition when is varied and approximately 42% when is. Once again, inhibition seems to be somewhat more effective in controlling the global dynamical properties of the cortex, in this case as driven by excitatory exogenous inputs.

### Changes within families due to thalamic input

To complete our analysis of 2par plots, let us make some remarks on the other instantiations of F1 and F2 that have not been illustrated yet, but do occur in the 405 analyzed sets. These diagrams are only slight, local variations of the ones shown so far, and still exhibit the global common characteristics typical of the respective families. As an example, patterns belonging to F1 very similar to Fig. 4A are present, but they miss the point. This type of diagram is still considered as part of F1, because the structure of its sn lines does not change. If we start looking at the transition between F1 and F2 from this particular plot, rather than from Fig. 4A, the fate of the pattern as is decreased is the same, i.e., it metamorphoses into an F2 diagram.

There is one final example of possible variations, regarding the appearance of Bogdanov-Takens (bt) points, which we anticipated in the previous subsection. Remarkably, these changes can be accessed as well by varying the same parameter : exogenous agents do not only produce transitions between families, but can also trigger modifications within the same family. We study here the annihilation of a couple of bt points, by increasing for the same parameter set employed in Fig. 4. Once again the reader is asked to look closely at Fig. 5 and its insets to follow our description. If we start at a low, physiological , two Bogdanov-Takens points, and , are present just below and we want to demonstrate how they disappear for larger . These bts occur on the left wedge of the sn line culminating in the cusp point . Notice also that the branch is responsible for the point close to , in analogy with the F2 in Fig. 3. The similarity between the two plots is evident. As the thalamic input is increased to in Fig. 5B, a couple of generalized Hopf (gh) points appear on . These points signal a change in the hb line from subcritical to supercritical, i.e. the Hopf bifurcation associated to the branch switches from a hard to a soft type ref:Kuz04.

If is made even larger, the two hb lines come closer and closer, up to a value of the input where they glue together and separate, as shown in Fig. 5C for . After this, acquires a large portion of what previously was , with the latter now extending over a much smaller distance. The split occurs in between and , each of which now belongs to a separate Hopf branch: the tiny reaches from to and the long bends in proximity of and continues to larger and values. Notice in the inset of Fig. 5C that also (in green) bends close to and heads back to , where it stops. Then, when is further incremented, see Fig. 5D with , annihilates and the smaller branch further shrinks, with the bt points getting closer and closer as grows. In the end, these two bt points coalesce, and we are left with only a single branch, as in Fig. 3B. We note that this line of hb delimited by and only exists at positive values of and . For the interested reader, we point out that the plot in Fig. 5A is responsible for chaos in an unphysiological range of and , and its route through a so-called homoclinic doubling cascade is analyzed in detail in ref:FvVB+08.

Let us conclude with some observations. First, the tiny branch exists for a very limited interval of and this is reflected in its very low occurrence in the 405 sets: it was found only twice. Second, an identical mechanism for F1 is also present, where the creation or annihilation of bt points takes place on the left lines in Fig. 3A. Finally, the diagrams that have been presented for the annihilation of bt points in F2, the similar plots for F1 (not shown), and the plots in Figs. 3 and 4 with their small, local variations, are the only ones that have emerged in our analysis of the reduced batch of 405 sets. Within the previously discussed assumptions, we conclude that these diagrams constitute all the possible modeled cortical responses to inhibition that are physiologically meaningful in Liley’s MFM.

### Partition of parameter space: families reaction to anesthetics and distributions of parameters

Having shown the properties of F1 and F2 types of diagrams, and concluded that any of the 405 plots of the preliminary sets falls into one of the two, we process the whole 73,454 sets in search of global, statistically relevant, correlations between family membership and parameter distributions. The approximation we introduce to cope with the large number of sets in an automated fashion is simple: F1 and F2 are those sets respectively without and with separate lines of sn, which are divided by the exchange of sn branches illustrated in Figs. 4C-D. Hence, Figs. 3A and 4A-C are considered as F1, Figs. 3B and 4D as F2. The second swallow tail bifurcation mentioned previously is the boundary between the so-defined F1 and F2. The algorithm we employed searches for the position and number of cusps on sn lines, and assigns the family type consequently. Equilibration, continuation and attribution of a single set to its family takes on average seconds on an everyday desktop computer.

The first, global difference between families we wish to show lies in their response to modeled anesthetic action. According to the methods explained in Bojak and Liley ref:BL05, it is possible to simulate the increase or decrease in EEG power when anesthetics are induced, and extract the ratio between the power of the anesthetized cortex and the cortex at rest. Herein power is predicted for the system at a stable equilibrium subject to fluctuations induced by noise added to the thalamocortical input . In line with clinical practice, they adopted the minimum alveolar concentration (MAC) of anesthetic agent at 1 atm pressure as the reference measure for describing the behavior of the anesthetized cortex. One MAC is the inspired anesthetic concentration needed to prevent movement in 50% of people to a noxious (surgical) stimulus. For example, assuming that isoflurane, a common anesthetics, is employed during surgery, a patient would be maintained usually at to isoflurane in an oxygen- nitrous oxide mixture, or at without the nitrous oxide. According to previous studies, it is assumed that 1 MAC is equal to a for isoflurane, which corresponds to aqueous concentration ref:Map76; ref:FL96.

The distributions of family responses to anesthetics according to power ratios are depicted in Fig. 6, with F2 types being in the majority in the whole batch, since they turn out to be about one and half times as many as F1. The histograms look similar in shape but are shifted with respect to one another. They are characterized by similar averages, with a relative difference of only about , but are strikingly different at the respective tails. At low ratios, F1 sets have a tendency of appearing more frequently than F2 sets, and vice versa for high ratios. For example, for power ratio values smaller than , the cumulative probability for F1 is almost three times that of F2, whereas for ratios bigger than the situation is reversed, with F2 having a cumulative probability more than three times that of F1. This shows that global, dynamical patterns have an interesting degree of correlation with total EEG power when responding to modeled isoflurane, whose mode of action is representative of most anesthetics.

Intuitively, anesthetized cortices should show a decrease in total EEG power when anesthetics are induced, as behavior is impaired as the cortex “slows down”, with power ratios at 1 MAC expected to be smaller than one. This is not always true, and it is not uncommon to rather observe a transient surge in EEG power that is clinically known as the biphasic response ref:BGC+97. For the generated 73,454 sets, a power increase larger than 1.4 at 1 MAC, has been taken as indicative of a biphasic response in ref:BL05, i.e., the simulated total power integrated over all frequencies was at least 1.4 times larger with than without the presence of 1 MAC isoflurane. Only 86 sets out of the whole batch show a biphasic response, of which around 90% are of F2 type. This is interesting, but the statistics are too poor to draw a general conclusion about the relation of F2 to the appearance of a biphasic response. It is possible that the predominance of F2 for high ratios is related to the cortical mechanism eliciting a surge in power when anesthesia is present.

Besides responses to anesthetics, our family partition also shows that some frequency distributions of parameters exhibit significant differences for F1 and F2. In order to objectively assess these differences, we use the following procedure: first, we construct parameter frequency histograms with bins within the limits provided for each parameter by Tab. 1, separately for the two families. Two particular cases arise in this procedure. On one hand for and , which have flexible upper limits, we choose the maximum possible one () as bin limit. On the other hand, occupies almost exclusively the lower part of the allowed parameter range. For the sake of better accuracy with the same number of bins, we set as the upper bin limit for . Parameter sets with larger values, i.e. 12 out of 29,131 for F1 and 24 out of 44,323 for F2, are counted in the highest bin. Second, we now compute the square root of the information radius (also known as the Jensen-Shannon divergence) between the histograms for each parameter – for F1 and for F2 – with bins

(10) |

which is a proper metric for the similarity of discrete probability distributions ref:ES03. Note that for empty bins, one defines . Identical distributions have , whereas maximal dissimilarity is indicated by . An example for maximal dissimilarity would be histograms and with Kronecker deltas and .

The eight largest dissimilarities between the frequency distributions of the families are found for the parameters , , , , , , and , where the number in brackets is the corresponding value. These distributions are shown in Figs. 7 and 8. In the case of , see Fig. 7A, F1 reaches a maximum around , with a slow decline for higher values, while F2 has a steadily increasing trend, perhaps saturating above . This is also the case for F2 sets in , see Fig. 7B, whereas F1 gives rise to a broad unimodal distribution with a maximum around mV. For both and , the difference in frequency for F1 and F2 at the edges of the intervals is striking. This strongly suggests that and play a relevant role in selecting the global dynamics of the model, with a clear bias for one type of family over the other.

Although differences in the distributions of other parameters are still significant, they do have a minor effect in predicting family membership, since the information radius distance of and is by far the largest of all the distributions. For example, in Fig. 7C and in Fig. 7D retain separated maxima for the occurrence of F1 and F2, showing a tendency for larger likelihood of F1 at low physiological values and a faster decrease in frequency as the parameters increase. Notice how in Fig. 8A shows bell-shaped distributions concentrated in the subinterval , in a similar configuration to that of total power ratios discussed above in Fig. 6. In contrast, for PSP amplitudes , , and connection strengths , see Fig. 8B-D, we observe larger changes for F2, while F1 is closer to constant throughout the allowed parameter range. In these three cases at the edges of the physiological parameter interval one finds clear bias for one family over the other, but overall the difference is not as strong as for and . The distributions of the remaining parameters do not show marked differences between F1 and F2, with values between 0.012 and 0.10.

### Different oscillatory responses of families to altered inhibition

Belonging to F1 or F2 only weakly predicts the presence of a surge or decrease in power when isoflurane action is simulated, since only the tails in Fig. 6 have a predominant presence of one family over the other. Nevertheless, the variabilities of patterns belonging to F1 and F2 embody different responses to changes in inhibition. Thus we focus now on the oscillatory activity associated with F1 and F2, when changes in inhibition are modeled solely via alteration of the parameter , to see how alterations in inhibitory PSP modify the oscillatory “landscape” families can produce. To this end we inspected recurrent 1par plots at fixed , within a physiological range of , and characterize the oscillatory activity.

Examples of this are depicted in Figs. 9 and 10, obtained from selected sets in the initial batch of 405. We want to stress that the variations within each family are extensive: the plots proposed are only schematic indications and are not meant to exhaustively describe the behaviors of F1 or F2 types. Nonetheless, some general conclusions can be drawn, based on the most commonly appearing dynamical traits in the continuations in . Firstly, F1 seems to be less prone to stable oscillatory activity in the physiological range than F2. Changes of R, and thus of and , tend to trigger unstable periodic orbits with a limited parameter span for F1, whereas F2 shows stable cycles in a larger parameter range. In particular, plots where the size of the parameter interval of stable orbits is either very limited (Fig. 9A) or zero (Fig. 9B) are frequently observed in F1 types, whereas plots with significant oscillatory regimes (Fig. 9C and D) are infrequent.

Secondly, the presence of two or more stable regimes is rarer in sets belonging to F1, as compared to F2, and furthermore they typically extend over more limited ranges within the physiological interval. Hence the occurrence of separate periodic orbits at the same value of is much rarer in F1 than in F2. In other words, multistability is distinctive of F2. For instance, the plots in Fig. 9C and D have a very limited extension of stable orbits for physiological values of compared with Figure 10B, C and D. Since a relation between multistability and memory has been hypothesized, this could suggest that F2 types possess dynamical structure that facilitates the formation of memory. In general and on average, it appears that F1 is associated with a more restricted dynamical repertoire than that of F2, and that F2 responds to the changes in inhibition with more activity than F1.

The distribution of the Hopf bifurcations in different families gives some insight into this. In F1, the curve of Hopf bifurcations can simply be closed, with no interaction with the saddle-node bifurcation. This is responsible for the dynamics in Fig. 9A and B. In this case, there exists a time-periodic solution but it is stable only on a small domain. A slightly more complicated situation arises if the Hopf bifurcation intersects with the saddle-node bifurcation at the fh points, as in Fig. 3A. This leads to a more folded branch of periodic solutions, as shown in Fig. 9C and D. Again, stable oscillatory behavior is observed only in a small domain. In F2 instead the hb curve is stretched out and interacts with both branches of sn curves and generally with both sides of the -shaped branch, as in Fig. 3B. This last interaction is, in particular, closer to the physiological range than the fh points usually observed for F1. Hence, the configuration in F2 does have a greater potential for sustained oscillations and the coexistence of stable periodic solutions, as shown in Fig. 10B.

## 5 Discussion

We have illustrated how bifurcation diagrams of Liley’s MFM in two parameters, and , account for core effects of inhibition over the cortex and can be classified into two topologically distinct families. These parameters modify inhibitory PSPs and inhibitory self-coupling, respectively, and the two archetypal families are identifiable by invariant dynamical features of model cortex. Across a large number of parameter sets, we have listed and identified all the possible diagrams. The relatively small number of bifurcations and combinations between branches we have found indicates that Liley’s MFM, when it supports “realistic alpha activity” at nominal parameter values ref:BL05, results in a relatively limited dynamical repertoire for cortical activity in response to parametric variations within physiologically admissible space. These dynamical structures have been classified into families according to their shared global topology and equivalent local properties of the solutions around their bifurcations. Family membership weakly correlates with typical changes in total spectral power when the cortex is acted upon by anesthetic agents, but strongly with unusually large changes. At the same time, the likelihood for family membership is sensitive to changes in specific neurophysiological attributes, in particular to the standard deviation of the excitatory neural population firing threshold (i.e. the steepness of the neuronal firing rate function ) and the excitatory membrane decay time constant . As far as we are aware, these connections between dynamical properties of a MFM, physiological attributes of cortex and changes in EEG power spectra under anesthesia have never been found before.

It is important to appreciate the generality and relative simplicity of the methods we have used here. We have chosen two meaningful parameters for the bifurcation analysis, and , and have then classified bifurcation diagrams into two families, F1 and F2, according to the absence or presence of separate sn lines. Family membership turned out to be separated by a swallow tail bifurcation. More sophisticated criteria for labeling resulting 2par plots could be envisaged. For example, families might have been categorized by the number, location and types of bifurcation they presented, or by the values of the coefficients in their normal forms. However, in our case the simple approach proved sufficient. The main methodological result we wish to underline is that by grouping the global dynamics for a representative sample of 73,454 parameter sets of an EEG model according to agreed criteria, strong, distinctive and robust correlations with EEG power spectra on one hand and some of the physiological parameters on the other hand have been derived. This is novel per se, because it shows that systematic partitioning of complex parameter spaces according to global dynamical patterns can lead to significant results which are not otherwise accessible. This method seems to be one the few viable strategy for unveiling interesting relations in high dimensional nonlinear spaces and exploring typical responses in complex models.

In this paper, we have focused on the action of anesthetics and ratios of total power, but many other developments are possible. For example, we could select those sets which give rise to epileptic activity or have a high power in the EEG gamma band, categorize their dynamical patterns and look for parameter correlations. Our approach is hence promising for other theoretical studies of the EEG, or neural activity in general.

We want to stress that the proposed method is not dependent on Liley’s MFM: in general, any other set of coupled, nonlinear ODEs with a complicated but meaningful parameter space can be treated likewise. The difficulties in extracting relevant information from equations with many parameters and involved feedback and feedforward mechanisms should not be underestimated, and our dynamical approach can offer new guidance. Here a PCA on the 73,453 analyzed parameter sets turns out to be completely useless. This is also true if we employ a PCA after dividing the sets into F1 and F2: the first ten principal components then account for about 55% of the total variance, up from 42%. Again, no valuable information on parameter relations can be gleaned from these types of statistical methods. If we compare Fig. 2 with the clear separation between distributions in power ratio in Fig. 6 or in Fig. 7 for parameters and , we can appreciate the strength of our method to tease out hidden correlations.

What is the biological meaning of family differences with respect to parameters such as and ? First, a small excitatory rate constant enhances the speed with which adapts to changes in PSP inputs, as is clear from Eq. 2. For , the mean excitatory soma membrane potential simply becomes a weighted sum of the PSPs. Thus, low enslave more directly to synaptic input and results in a depression of self-sustained oscillations, which is characteristic for F1. parametrizes the local excitatory firing rate response, see Eq. 3, and hence does not directly act upon dynamics, as does. However, a smaller means that the excitatory firing rate changes more rapidly for variations of , although over a more limited range. In the limit of , one obtains a step function which instantaneously switches the population from zero to maximal firing at the threshold value . As mentioned above, in a steady state situation excitatory self-feedback is essentially proportional to , and a low increases this self-feedback. Thus, as for , the reaction of to changing synaptic input will be more rapid for low , limiting the presence of self-sustained oscillations and hence favoring membership in F1. One can then view the somewhat lesser differentiation of with regards to the families due to its boosting excitatory reactivity more indirectly. Overall, enhancements in the reactivity of excitatory populations result in a diminished dynamical repertoire.

These considerations on the primary role of and in influencing family membership are not weakened by the presence of complex feedbacks between neuronal populations in Liley’s model, c.f. Fig. 1. Sets of F1 types appearing at low or show homogeneous distributions for all the other parameters, with no formation of clusters or any other form of bias. This suggests that the model cortex generally can be pushed towards F1 by appropriately reducing the excitatory membrane decay constant or the excitatory standard deviation of the firing threshold, and towards F2 by increasing them. Thus, we can also identify excitatory reactivity as an endogenous counterpart to the exogenous control. In fact, it has been shown in detail here also how exogenous activity, namely the excitatory thalamic input to inhibitory (excitatory) cortical populations described by (), can cause changes across and within families. There exists a kind of mirror symmetry between and : at the increase of , or decrease of , we have observed a transition from F2 to F1, see Fig. 4, corresponding to a transition from complex to simple dynamics. This indicates that the exogenous driving of inhibition (or a reduction in the driving of excitation) moves the cortex to dynamics which, on average, appear simpler.

Thalamic input appears in Liley’s MFM as an active source of control on the complexity of brain dynamics, having the role of modulator in cortical interactions. This is quite different to the classical view of thalamus as being the gateway through which peripherally derived sensory information reaches cortexref:W38. For example, retinal axons project to the lateral geniculate nucleus of the thalamus, which then projects to visual cortex. In this way the thalamus is seen as the conduit through which all information about lower levels of the nervous system reaches the cortex for “processing”. Viewed from the perspective of our model, this would imply that thalamic input simply drives the cortex. However, based on our metabifurcation analysis we predict that such input also configures (i.e. “modulates”) local and global dynamical responses from the cortex.

It is now known that all areas of cortex, not just the primary sensory cortices (somatosensory, visual, auditory) receive thalamic input, and that most areas of cortex, in some reciprocal manner, innervate thalamusref:SG06. As a consequence, contemporary approaches to understanding the role of thalamus are now aimed towards better articulating the functional inter-relationships that exist between cortex and thalamus. One popular approach has been the modeling of integrated dynamics of thalamo-cortical activity. On this basis, it has been speculated that the human alpha rhythm emerges as a consequence of reverberant activity between cortex and thalamus. While we have not here considered such thalamo-cortical feedback, doing so in the context of our results leads to a number of interesting speculations regarding the auto-regulation of cortical dynamics: cortical feedback through the thalamus could initiate a sequence of transitions between bifurcation pattern families, providing a way of reconfiguring cortical dynamics “on the fly” and thus on a time scale quite different to that associated with activity dependent synaptic plasticity. Unfortunately, not enough is known experimentally about thalamocortical feedback at present to turn such speculations into strong model predictions.

On the other hand, thalamus also appears to stimulate internal variations within families, since a lower value of (higher value of ) corresponds to more complex bifurcations as depicted in Fig. 5. As said, this type of 2par plot with two bt points has been shown to produce chaotic activity, meaning that “cortical state” F2 is capable of producing highly nonlinear events, and an exogenous enhancement of cortical inhibition tends to reduce this ability. This and the fact that F2 types have shown to occur with the largest likelihood in parameter space (i.e., their total number is one and a half time more than F1), stimulates some speculations. F2 may be considered as the cortical default state of Liley’s MFM, associated with the richest dynamical responses to inhibition and acting as a sort of cognitive readiness state. Both endogenous (cortical) enhancements of excitatory reactivity and exogenous, subcortical increases (decreases) of inhibitory (excitatory) activity reduce the cortical state to much simpler behavior.

Finally, we have implemented a partition of the parameter space according to the bifurcation plots. As mentioned, types of normal forms appear to be constant within each family and continuation lines of the diagrams are similar, whilst the extensions of these patterns in the -space vary wildly. This suggest the hypothesis that a single, high dimensional master bifurcation diagram exists and governs all sets, and that changes in parameters only change its orientation, position and scale in the (hyper)volume. By continuing in and , we intersect the master diagram with a plane. The position of this plane determines whether we see F1 or F2 diagrams, and the two can be continuously deformed into each other.

We would like to conclude with some future directions. First, an inventory of all the possible types of oscillatory dynamics associated with F1 and F2 is surely attractive. How do many of them support complex events such as chaos and how is this reflected in parameter space? Do chaos and complex dynamics occur for limited, selected intervals in the parameters? The great hurdle is represented by the computational demand: it is still unfeasible to process all sets one by one and see what types of orbits are present. Second, a systematic study of the unfolding of the bifurcations in 2par plots that characterize families should be attempted. This has the potential to increase our knowledge about the type of dynamics these parameter sets can produce. Finally, there is an interesting theoretical question concerning the optimization of the method we presented in this paper. What is the best choice of continuation parameters for achieving the highest degree of separation among bifurcation plots? We wonder if it would be possible to choose parameters or combinations of parameters such that their degree of clustering with respect to families is maximal. Some canonical discriminant analysis on the sets has been preliminarily performed, but results are not particularly illuminating. Achieving optimal partitioning has the potential to uncover biologically relevant dependencies among parameters that are still unknown, and to shed more light on clinically or pharmacologically relevant questions.

## 6 Acknowledgements

FF was supported by Australian Research Council Discovery Grant DP 0879137 and the Swinburne University of Technology Researcher Development Scheme 2010. FF would like to acknowledge the hospitality of the Department of Mathematics at Concordia University, Canada, where the idea for this study was conceived. LvV was supported by NSERC Discovery Grant 355849-2008. IB would like to thank Prof. Rolf Kötter for the time granted to complete this project.