Extreme value theory of evolving phenomena in complex dynamical systems: firing cascades in a model of neural network

TL;DR

This paper extends extreme value theory to evolving phenomena in complex systems, exemplified by neural network cascades, providing new statistical laws for analyzing such events.

Contribution

It introduces a novel framework for applying extreme value theory to non-instantaneous, evolving phenomena in complex dynamical systems, with a focus on neural cascades.

Findings

Derived extreme value laws for neural cascades.

Proposed a new definition of neuronal cascade.

Analyzed cascade statistics using exceedances and block maxima.

Abstract

We extend the scope of the dynamical theory of extreme values to cover phenomena that do not happen instantaneously, but evolve over a finite, albeit unknown at the onset, time interval. We consider complex dynamical systems, composed of many individual subsystems linked by a network of interactions. As a specific example of the general theory, a model of neural network, introduced to describe the electrical activity of the cerebral cortex, is analyzed in detail: on the basis of this analysis we propose a novel definition of neuronal cascade, a physiological phenomenon of primary importance. We derive extreme value laws for the statistics of these cascades, both from the point of view of exceedances (that satisfy critical scaling theory) and of block maxima.

Equations

$$H(z_n,z_{n+1},\dots ,z_{n+m})=H(z_n,{\phi }^1(z_n),\dots ,{\phi }^{m-1}(z_n))=\mathcal{H}(z_n).$$

$$\left\{\begin{array}[]{ll}x^{i}_{n+1}$ = $F(x^{i}_{n},y^{i}_{n},I^{i}_{n},x^{i}_{n-1})\\ y^{i}_{n+1}$ = $y^{i}_{n}-\mu(x^{i}_{n}+1)+\mu\sigma+\mu I^{i}_{n}.\end{array}\right.$$

$$F(x_{n},y_{n},I_{n},x_{n-1})=\left\{\begin{array}[]{ll}\frac{\alpha}{1-x_{n}}+y_{n}+\beta I_{n}&$ if $x_{n}\leq 0\\ \alpha+y_{n}+\beta I_{n}&$ if $x_{n}>0$ and $\{x_{n}<\alpha+y_{n}+\beta I_{n}$ and $x_{n-1}\leq 0\}\\ -1&$ if $x_{n}>0$ and $\{x_{n}\geq\alpha+y_{n}+\beta I_{n}$ or $x_{n-1}>0\}\end{array}\right.$$

$$\xi_{n+1}=\left\{\begin{array}[]{ll}1$ if $x_{n}>0$ and $\{x_{n}\geq\alpha+y_{n}+\beta I_{n}$ or $x_{n-1}>0\}\\ 0$ otherwise$.\end{array}\right.$$

$$\omega_{ij}=\left\{\begin{array}[]{ll}0&\mbox{ if neuron }j\mbox{ does not act on neuron }i,\\ 1&\mbox{ if neuron }j\mbox{ has excitatory action on neuron }i,\\ \psi&\mbox{ if neuron }j\mbox{ has inhibitory action on neuron }i.\end{array}\right.$$

$$I_{n+1}^i-\eta I_n^i=W\sum _{j=1}^N{\omega }_{ij}({\chi }_{ij}-x_n^i){\xi }_n^j.$$

$$F_H(h)=\underset{L\to \infty }{\lim }\frac{1}{L}\#\{H({\phi }^n(z_0))\le h,0\le n<L\}.$$

$$E_H(h)=\mu (\{z\in Z\text{mbox}s.t.H(z)>h\})$$

$$\begin{array}[]{ll}(n^{\prime},i^{\prime})\rightarrow(n,i)\Longleftrightarrow&\xi_{n^{\prime}}^{i^{\prime}}=1,\;\xi_{n}^{i}=1,\\ &\exists\;n^{*},\;n^{\prime}<n^{*}<n,\mbox{ such that the following holds:}\\ &x_{m}^{i}<x_{th}\mbox{ for }n^{\prime}\leq m<n^{*},\\ &x_{m}^{i}\geq x_{th}\mbox{ for }n^{*}\leq m<n,\\ &\xi_{m}^{j}=0\mbox{ for }n^{\prime}<m<n^{*},\;\forall j\mbox{ such that }\omega_{i,j}=1.\end{array}$$

$${\rho }_S(S)\sim S^{{\beta }_S},{\beta }_S=\frac{{\beta }_G+\gamma -1}{\gamma },$$

$$I(z,n)=H(z_n)=H({\phi }^n(z))$$

$$\begin{array}{rl}{H}_L(z)& =\max \{I(z,0),I(z,1),\dots ,I(z,L-1))\}\\ & =\max \{H(z),H(\phi (z)),\dots ,H({\phi }^{L-1}(z))\}.\end{array}$$

$$F_{H_L}(h)\sim (F_H(h){)}^{\lambda L}.$$

$$\begin{array}{rl}\log (F_{H_L}(h))& \sim \lambda L\log (F_H(h))\\ & \sim \lambda L(F_H(h)-1)\\ & =-\lambda L{E}_H(h).\end{array}$$

$$F_{H_L}(h)\simeq \exp \{-\lambda L{A}_H{h}^{-{\beta }_H}\}.$$

$$F_{H_L}(h)\simeq \exp \{-\lambda L\exp \{A_H{h}^{B_H}\}\}.$$

$$F_{H_L}(h)\sim \exp \{-\lambda L{E}_H(h)\}.$$

$$\lambda L{E}_H(h_L(t))\to t,\text{mbox}asL\to \infty .$$

$$\mu (\{z\in Z\text{mbox}s.t.H_L(z)\le h_L(t)\})=F_{H_L}(h_L(t))\to F(t)=e^{-t},\text{as }L\to \infty .$$

$$h_L(t)\sim E_H^{-1}(\frac{t}{\lambda L}).$$

$$h_L(t)\sim {\left(\log (\frac{\lambda L}{t})\right)}^{\frac{1}{B_H}}.$$

$$h_L(t)\sim (\frac{\lambda L}{t}{)}^{\frac{1}{{\beta }_H}}$$

$$\mu (\{z\in Z\text{mbox}s.t.(\lambda L{)}^{-1/{\beta }_H}{H}_L(z)\le u\})\to e^{-u^{{\beta }_H}},$$

$$\left\{\begin{array}[]{lcl}x&\rightarrow&\frac{\alpha}{1-x}+y\\ y&\rightarrow&y-\mu(x+1)+\mu\sigma.\end{array}\right.$$

$$\left\{\begin{array}[]{lcl}\bar{x}&=&\sigma-1\\ \bar{y}&=&\sigma-1-\frac{\alpha}{2-\sigma}\end{array}\right.$$

$$J=\left[\begin{array}[]{lcl}\frac{\alpha}{(2-\sigma)^{2}}&&1\\ -\mu&&1\end{array}\right]$$

$$\text{mbox}Det(J)=A(\alpha ,\sigma )+\mu .$$

Full text

*keywords: dynamical extreme value theory, complex systems, neural networks, criticality, neuronal avalanches

1 Introduction

The dynamical theory of extreme values, well described in a series of works (see e.g. the review book [1] and references therein), has so far mainly considered low dimensional systems in a very simple setting that we will describe momentarily. To the contrary, its statistical ancestor, the so–called extreme value theory (EVT) [2], has been fruitfully applied to phenomena like earthquakes, floods, epidemics [3, 4], that all involve complex, many dimensional interactions that operate for a certain period of time and are otherwise silent. It is the purpose of this paper to develop a theoretical framework to extend the dynamical treatment of extreme events to this kind of systems. We shall do this by focusing on a specific example, a model of neural network that has been introduced to describe various phenomena occurring in the electrical activity of the cerebral cortex [5, 6, 7, 8].

This model draws its origin, object and motivation from a vast body of existing research involving mathematics, physics and physiology. Nonetheless, the goal of this work is not to reveal new phenomena in the last field—although we will propose a new conceptual approach to a much studied problem, the statistical description of firing cascades, better known as neuronal avalanches [9, 10, 11, 12, 13]. Rather, the model is particularly suited to introduce a formal advancement in the dynamical theory of extreme values, which is also relevant to the neuronal case, among many others. In fact, at difference with standard theory, we will consider extreme events that occur in complex, many dimensional systems and that are not restricted to a specific instant in time, but which evolve during a finite time interval.

1.1 Dynamical Extreme Value Theory

The dynamical theory of extreme values can be summarized by considering a phase space $Z$, in which motions are determined by the repeated action of a deterministic map $\varphi:Z\rightarrow Z$ (that can also be probabilistic, but we will restrict ourselves to the first case). In addition, an intensity function $H:Z\rightarrow\mathbb{R}$ gauges the value, i.e. the intensity $H(z)$ of a phenomenon, an event, which occurs at the precise instant of time when the motion visits the point $z\in Z$. Extreme value theory is the study of the statistics of particularly intense events, i.e. large values of $H$, accumulated in time along trajectories of the motion. Without loss of generality we shall assume that time is discrete and denote it by the letter $n$. Trajectories of the system will naturally be written as $z_{n}=\varphi^{n}(z_{0})$ and $I(n,z_{0})$ will denote the intensity of the event occurring at time $n$ after the motion started at the point $z_{0}$: $I(n,z_{0})=H(\varphi^{n}(z_{0}))$.

Originally [1], very simple intensity functions $H$ have been investigated, which depend only on the Euclidean distance of $z$ from a single point $z^{*}$, and increase as this distance diminish. In this simplified setting extreme events are close approaches of the motion to such location. This situation was widely generalized by one of the present authors in [14], by introducing intensity functions that depend on the distance from an uncountable set of points populating a fractal set $K\subset Z$. This analysis has recently been confirmed in full mathematical rigor [15].

In addition, recall that phase space $Z$ and map $\varphi$ are complemented by an invariant measure $\mu$ to fully define a dynamical system. The measure $\mu$ determines the frequency by which the dynamics visit different regions of phase space. The combined action of this measure and of the geometry of $K$ in determining extreme value statistics was also investigated in [14], revealing the rôle of fractal quantities like Minkowski dimension and Minkowski content, which were extended to the case of general (i.e. not necessarily Lebesgue and possibly singular) invariant measures.

Suppose now that the intensity of a phenomenon cannot be determined instantly at time $n$, when the dynamics visits the phase–space point $z_{n}$, because this phenomenon develops over a certain time span, so that its true intensity can only be assessed when it terminates. This is the case of neuronal avalanches, when the electrical activity of the brain is sustained for a few milliseconds. But it is perhaps easier, to focus ideas, to consider a geophysical analogy [16], to which the theory exposed herein can be equally applied.

An earthquake can be triggered by a single seismic event, a rupture, which propagates and multiplies in space and time. The total energy released by the quake is the intensity of the event, which can be assessed only after the shock terminates. Notice at this point that even if this event takes place continuously over a finite time–span, it is usually ascribed in geophysical catalogs to the moment of the first, triggering rupture. To put this observation into the theoretical framework of dynamical extreme values, the intensity function must take the form $H:Z^{m}\rightarrow\mathbb{R}$, in the sense that the event at time $n$ has an intensity $I(n,z_{n})$ of the form $I(n,z_{n})=H(z_{n},z_{n+1},\ldots,z_{n+m-1})$, where $m$ is the duration of the phenomenon. This generalization of the conventional theory will be studied in this paper.

A first observation must be made at this point. Since we are in the presence of a deterministic dynamics, $z_{n+j}=\varphi^{j}(z_{n})$, the intensity function takes on the special form

[TABLE]

Intensity can therefore be thought of as depending only on the initial point $z_{n}$ via a global function ${\cal H}:Z\rightarrow\mathbb{R}$. This might lead one to think that the proposed approach is void of any generalization. As a matter of facts, because of the dependence of $H$ on future points of the evolution, the function $\cal H$ is presumably very complex and its phase–space portrait, i.e. its level sets, might feature a fractal, hierarchical organization, thereby providing a concrete dynamical example of the theory of [14] mentioned before. Moreover, when the dynamical system is chaotic, the function $\cal H$ in formula (1) cannot be given a more efficient computational scheme other than what expressed by the central expression in Eq. (1): the best one can do is to follow the dynamics as it evolves and record its data. Think again of the earthquake example: while it develops, it is practically impossible to forecast its intensity. On a more theoretical vein, recall the aphorism by Joe Ford: a chaotic system is its own fastest computer and most concise description [17].

There is a further consideration that renders the generalization (1) far from trivial: the time–span of the phenomenon, $m$, is not known a priori, nor it is constant over different events: the integer $m$ in equation (1) depends on $z_{n}$. Because of the reasons just outlined there is no significantly better way to compute $m$ than to follow the dynamics until $m$ reveals itself: that is, the phenomenon under investigation—the quake—terminates.

1.2 Neuronal avalanches

It is now possible to return from the geophysical analogy to the paradigmatic example described in this paper: neuronal avalanches. It is observed that electrical stimuli originate (in the presence of, or without a sensorial stimulus) and propagate in the cerebral cortex to involve the firing of thousands of neurons. In fact, a neuron may either be silent, or fire, which means, polarize and send an electric stimulus via synaptic connections to other neurons. There is evidence that this activity composes a fundamental part of the computation performed by the brain [18, 19, 20]. It can be recorded by suitable devices [11] and it appears as bursts of firings, at times well separated between each other, but more often adjacent. In certain instances this activity may be due to synchronized oscillations of a large number of neurons [21], but more complicated patterns are frequently observed, somehow intermediate between the former and homogeneous, uncorrelated firing of individual neurons. A substantial amount of literature (see e.g. [9, 10, 11, 12, 13]) has shown that this regime can be described as critical, characterized by power laws and scaling relations. We focus our attention on both the critical and the subcritical regime, but we introduce a seemingly new operational definition of avalanche.

This definition is akin to what Rangan and Young [22, 23] have termed multiple firing events, as the maximal sequence of firings that are caused by a single neuron. By simulating each neuron in a network (typically, several thousands of these are considered) as characterized by microscopic variables (voltages and conductances) evolving according to coupled ODEs, and by the technical expedient of setting certain synaptic time constants to zero, Rangan and Young have been capable of determining exactly which neuron causes the firing (also called the spike) of a post-synaptic neuron.

The situation in the model that we will study is different and more general, in that the combined electrical activity of many different pre–synaptic neurons can cause the spike of a neuron. Nonetheless, by studying the detailed dynamics of a single, unperturbed neuron—in particular, the basin of attraction of the spiking event—we will be able to introduce a rule according to which the firing of a particular neuron can be identified as the precipitating cause of the spike of another. We shall then define a firing cascade as the tree that mathematically encodes such chain of excitations. We believe that this approach is general enough to be applicable to a wide set of complex systems, composed of many interacting components.

Typical characteristics of these trees/cascades, like their intensity (number of nodes/number of firings), duration, and degree (number of levels in the tree) provide intensity functions $H$ to which extreme value theory can be applied, both from the point of view of exceedances and of block maxima (concepts to be made precise below), thereby providing a concrete example of the abstract theory described in the previous section.

This analysis will produce scaling relations that also permit to compute the probability that all cascades in a certain time interval have intensity less than a given threshold. This kind of extreme value laws for the size of neuronal cascades has not been investigated before, to the best of our knowledge—see however [24] for a discussion of the relevance of extreme value theory to sensorial perception.

1.3 Contents of the paper

The plan of this paper is the following. We first describe in Sect. 2 the dynamical model of neural network that we adopt, borrowed from [8]. It consists of a collection of two–dimensional maps, each representing an individual neuron, initially proposed by Rulkov [5, 6]. In sections 2.1 and 2.3 the theory follows closely, with minor modifications, the original references [5, 6, 7, 8]. This is a necessary introduction, not only to set notation and define concepts, but also to investigate a few details of dynamical import: in particular, the basin of attraction of the spiking event, described in Section 2.2, which is instrumental to define causality interactions. Neurons do not work in isolation, but are linked in a sparse network, where different neurons interact only when the spiking state is reached. In Sect. 2.3 we describe the model of network connections, following [7, 8], and then in Sect. 3 the collective dynamics that takes place on it.

We start in Sect. 3.1 from the problem of the existence of an invariant measure that we heuristically replace by Birkhoff time averages, having nonetheless in mind that this procedure might be dubious due to possible long transients in the dynamics [25]. Section 3.2 contains the fundamental definition of firing cascades as trees of causal excitations, as well as its practical implementation, which is based on the concept of basin of attraction of the spiking event. This permits to define and measure important quantities like the intensity $S$ of a cascade, its degree/generation level $G$ and its time span $T$. In Sect. 4 the statistics of these quantities is analyzed, showing a transition to critical behavior as the intensity of coupling among different neurons is increased; the exponents of the power–laws describing the above quantities obey critical scaling theory and agree with previous experimental findings and theoretical research. From the point of view of extreme value theory, this analysis can be classified as the study of exceedances.

Finally, in Sect. 5 we further widen the scope of our investigation: we show the emergence of Extreme Value Laws in the statistics of extrema in time–intervals of finite length. That is, we will derive scaling relations for the distribution of the largest event in a certain time interval, both in the critical and the sub–critical regime. Numerical experiment will confirm the formal theory. In the considered example, this offers a new statistical characterization of neuronal cascades, which might perhaps be of physiological relevance [24].

The conclusions briefly summarize the work, and an Appendix describes in an elementary way the stability regions, in parameter space, of the fixed point of Rulkov map.

2 The dynamical model: a sparse network of Rulkov maps

A variety of models have been proposed to study numerically neuron dynamics. Many of these are random. Purely dynamical models, on the other hand, can be continuous time models, involving systems of differential equations [26]. They are able to reproduce a whole variety of observed biological patterns [22, 23, 27], but are expensive from a computational point of view. On the contrary, Rulkov 2d–map [5, 6] is a good compromise between capacity of well reproducing real biological behaviors, such as irregular spiking and bursts of spikes [6, 7] and fast computations, two features that are crucial for analysis. In this section, we briefly describe the dynamics of a single neuron, which is effected by Rulkov map. We follow the implementation presented in [7, 8] that also contains illustrative pictures. We need only to mention that, at difference with these references, we do not allow a random external input: our dynamical system is fully deterministic.

Rulkov map describes the evolution of a fast variable $x_{n}$ ($n$ is integer time) that models the membrane potential of a neuron, and a slow control variable $y_{n}$ that, although it is void of physiological meaning, is required to pilot the evolution of the former, so to reproduce some patterns observed in real experiments. This pair of variables describes the state of every neuron in the system. We label neurons by the variable $i=1,\ldots,N$, so that the phase space of the system is $Z=\{(x^{i},y^{i}),i=1,\ldots,N\}$. Recall that $n$ labels time; Rulkov map is then defined as follows:

[TABLE]

The variable $I^{i}_{n}$ is the synaptic input, a current, acting on neuron $i$ at time $n$, defined below in Eq. (6). Different neurons are coupled only via this term. It is convenient to first discuss the dynamics when this external input is set to a constant (that may be taken to be zero, for sake of definiteness) and each neuron evolves independently of the others.

2.1 Unperturbed dynamics of the single neuron

In this subsection, accounting for the fact that neurons do not interact, we drop the neuron index $i$ for clarity. The function $F$ in Eq. (2) is defined as follows:

[TABLE]

The parameters $\alpha,\beta,\mu,\sigma$ appearing in equations (2),(3) are all positive.

The parameter $\beta$ modulates the relative influence of the synaptic input $I_{n}$ with respect to $y_{n}$ in Eq. (3). We choose $\beta=0.133$. Notice that a memory effect is present in the above definition, via the value of the variable $x$ at time $n-1$.

The parameter $\mu$ is taken to be small ($\mu=10^{-3}$), so that the variable $y$ evolves slowly with respect to time, following the second equation in (2). Also, its values cover a limited range. The combination of $\mu$ with $\alpha$ and $\sigma$ determines the structure of motions, as discussed in detail in [5, 6]. In Appendix we sketch a simple analysis of the stability of the fixed point $(\sigma-1,\sigma-1+\frac{\alpha}{\sigma-2})$. Stoop et al. [8] choose $\alpha=3.6$, so that only $\sigma$ is left free to vary. As shown in Appendix, under these conditions the map admits a bifurcation at $\sigma=\sigma_{cr}=2-\sqrt{\frac{\alpha}{1-\mu}}$. Parameter space is accordingly parted in four regions, I to IV, plotted in Fig. 13 in the Appendix. Choosing $\sigma=0.09$ puts the parameter-space point in region II, implying an attracting fixed point, while $\sigma=0.103$ leads to region III, where the fixed point is unstable. This permits to define two kinds of neurons: Intrinsically Non-Spiking Neurons (INSN) (parameters in region II, motion in the left frame in Fig. 1) and Intrinsically Spiking Neurons (ISN) (region III, right frame). These latter are the initiators of the dynamics of the network. The different behaviors are apparent in Fig. 1 where we plot motions in the single–neuron phase space (as in Fig. 1 in [8], the $y$ axis is in abscissa and the $x$ axis in ordinate).

A neuron is said to fire when its membrane potential $x_{n}$ reaches a local positive maximum (depolarization) and then suddenly drops to the hyperpolarized state $x=-1$: in physiology one talks of increasing polarization when the potential becomes more negative. At this moment a signal is transmitted to the other neurons, as described in the next two sections. It is instructive to follow the motion in the right panel of Fig. 1. Three cycles are plotted (two are almost coincident and require attention to be distinguished), each of which contains a spike. Start from the rightmost part of a cycle and follow points rotating counter-clockwise, so that $y_{n}$ initially decreases and $x_{n}$ increases. This happens because of the first alternative in Eq. (3). The map (2) is defined in such a way that when $x$ becomes positive at time $n$ ($x_{n-1}$ being negative: draw a horizontal line at ordinate $x=0$ in the figure for visualization), either $x_{n}$ is larger than $\alpha+y_{n}+\beta I_{n}$, so that at the next step the membrane potential drops to the negative value $x_{n+1}=-1$ (third alternative), or else it further increases to the value $x_{n+1}=\alpha+y_{n}+\beta I_{n}$ (second alternative), before achieving the value $x_{n+2}=-1$ at the following iteration, again in force of the third alternative. This motivates the introduction of the spike variable $\xi_{n}$ that is always null, except at those times when the neuron potential reaches the hyperpolarized state, in which case it takes the value one. It can be formally defined as follows (see eq. (3)):

[TABLE]

The spike variable $\xi_{n}$ will be used in Sect. 2.3 to define the dynamical action of synaptic coupling.

2.2 Perturbed dynamics of the single neuron: basin of attraction of the firing event

The synaptic input $I_{n}$ received by a neuron can either be stochastic or deterministic, as in our case. How this latter is generated in our model is inessential for the present discussion and will be explained in the next section; suffices here to observe that $I_{n}$ enters equation (3) linearly. In particular, a positive current $I_{n}$ increases the potential $x_{n}$ (hence, it depolarizes the neuron) and can therefore be associated with an excitatory action, while the reverse is true for a negative current $I_{n}$. Notice however that this rule of thumb has non–trivial exceptions that we will describe momentarily.

Figure 2 shows how the dynamics of ISN and INSN observed in Fig. 1 is perturbed by $I_{n}$. Typical trajectories are shown, in which sharp angles appear when $I_{n}$ is non–null. Notice however that time being discrete, lines in the figure join distinct points–albeit typically very close, so to simulate a continuous trajectory that does not exist. In the left frame one can observe how INSN, silent when isolated, can be triggered to spike by the input $I_{n}$. The dynamics of ISN (right panel) can also be perturbed, altering their regular spiking regime. In fact, consider Fig. 3, left panel, where the potential $x_{n}$ of an ISN (red curve) is plotted versus time $n$, together with the synaptic input $I_{n}$ it receives (green curve). The synaptic current $I_{n}$ is positive, stemming from an excitatory action. Three excitations are visible in the range displayed. Prior to these the ISN spikes regularly with a period of about 242 iterations. The first excitation has the effect of anticipating the successive spike by 35 iterations. In the right panel, which shows the same motion in the single–neuron phase space, regular oscillations appear as a thick limit cycle. The first spike is evident in the thread which emerges from the cycle at the right of the picture, and moves initially by increasing both $x_{n}$ and $y_{n}$, according to the fact that the excitation $I_{n}$ is positive: see Eq. (3). Returning to the left panel, we notice two more excitatory inputs to the neuron that have the paradoxical effect of suppressing further spikes! This can be explained by detecting, in the right panel, two more sharp angles in the trajectory, the first starting form the bottom left of the picture. In both cases $x_{n}$ and $y_{n}$ increase, but this has the effect of pushing the trajectory towards the fixed point of the map, which is unstable but from whose neighborhood the motion requires a long time to escape.

We have so discovered that excitatory actions can sometimes produce an inhibition of spikes, due to the particular nature of the phase space of the system. For a similar reason, inhibitory actions can sometimes induce spikes. This paradox leads to a deeper understanding of the dynamics of this system.

In fact, consider an INSN, like those pictured in the left panels of Figs. 1, 2. Let $(y,x)$ be the initial condition of the motion and let $I_{n}=0$ (isolated neuron). Even if finally attracted by the fixed point, some trajectories will spike before that happens. In the left frame of Fig. 4 we color in red such initial conditions, and in green those that head directly to the fixed point, whose position is marked by a blue asterisk. A limit cycle separates the two regions. In the right frame, the phase–space portrait near the stable fixed point is displayed. It is then clear that an inhibitory action (reducing both $x_{n}$ and $y_{n}$) can move a point from the green to the red region and that an excitatory kick can do the opposite, but it is equally clear that this occurs only in a restricted region of phase space. In fact, increasing $x$ typically moves a phase–space point from the green to the red region, i.e. from the basin of the fixed point to that of the spike. In particular, when the potential $x$ of a NISN is larger than a certain threshold $x_{th}$, the paradoxical effects just described are not possible. In fact, all initial conditions with $x_{n}$ larger than $x_{th}$ lie in the red region and lead to firing, even in the absence of further excitation: the effect of the current $I_{n}$ that puts a NINS in the red region can therefore only be reversed by an inhibitory input. This fact will be of paramount importance in the following. From the analysis of Fig. 3 we deduce that, in the range of parameters of our experiments, we can safely take $x_{th}=-0.7$.

2.3 Network topology and synaptic coupling

According to the way in which they act upon others, neurons can either be classified as excitatory or inhibitory. Following [8] we consider a network of size $N$ with a ratio of inhibitory/excitatory neurons $r$, so that the first $n_{ex}=\lfloor{Nr}\rfloor$ neurons are excitatory and the remainder $N-n_{ex}=n_{in}$ are inhibitory. Furthermore, as described in the previous section, neurons can either be intrinsically spiking or not, so that each group (excitatory and inhibitory) can also be parted in two: out of the $n_{ex}$ ($n_{in}$) excitatory (inhibitory) neurons $m_{ex}$ ($m_{in}$) are chosen to be intrinsically spiking.

Neurons are connected so that each of them is acted upon by $p_{ex}$ excitatory neurons and $p_{in}$ inhibitory neurons, chosen at random in the network—with no regard on whether ISN or NISN, with the only constraint that no neuron can act upon itself. The six integer numbers $n_{ex}$, $m_{ex}$, $n_{in}$, $m_{in}$, $p_{ex}$, $p_{in}$, plus the connection diagram fully specify the network. Analytically, this is formalized by the adjacency matrix $\omega_{ij}$, defined by

[TABLE]

The factor $\psi>0$ serves to differentiate the intensity of inhibitory with respect to excitatory action. We choose $\psi=3$.

Coupling in the network is driven by action potentials generated by firing neurons. Following again [8] the spike of a presynaptic neuron $j$ determines the input $I_{n}^{i}$ received by the neuron $i$ according to the linear difference equation:

[TABLE]

$W$ is a crucial intensity parameter, varying which we can increase the effect of synaptic excitation: recall that $I^{i}_{n}$ enters the dynamical equation of evolution of neuron $i$, Eqs. (2),(3). The form of Eq. (6) can be explained in physiological terms. First observe that the homogeneous equation (with null r.h.s.) implies an exponential decay: $I_{n}\sim\eta^{n}$ (we choose $\eta=0.75$). The spike variables $\xi_{n}^{j}$ at r.h.s. are the source of the currents $I^{i}_{n}$. In Eq. (6) we notice that when the synapse reversal potential $\chi_{ij}$ is larger than the potential $x_{n}^{i}$, its activation (that happens when $\xi_{n}^{j}=1$) increases the synaptic input $I_{n}^{i}$ and also, because of Eq. (3), the membrane potential $x^{i}_{n+1}$. This is equivalent to a depolarization of the neuron and therefore the synapse acts as excitatory. In the opposite case, the activation polarizes the neuron and the action is inhibitory. The numerical values of the above parameters are chosen as in [7, 8], where their physiological relevance is discussed: the reversal potential is null for excitatory connections (i.e $\chi_{ij}=0$ when neuron $j$ is excitatory) and equal to $\chi_{ij}=-1.1$ when neuron $j$ is inhibitory, for all $i$. This ends the description of the laws of motion of the network.

3 Network dynamics

We can now study the dynamical behavior of the full network. We start by discussing a basic problem: the existence of an invariant measure. Then, we turn to the definition of firing cascades.

3.1 Invariant measure

Equations (2),(3),(4),(5) and (6) define a multi-dimensional dynamical system in the phase space $Z=\mathbb{R}^{2N}$ with variables $(x^{i},y^{i}),i=1,\ldots,N$. Let $z=(x,y)$ be such $2N$–component phase–space vector and let $\varphi$ be the dynamical evolution operator. A fundamental ingredient of dynamical theory are the so–called invariant measures, that is, positive Borel measures that are preserved by the dynamical evolution. Lacking a theoretical control of these measures in many–dimensional systems like the one under consideration, we substitute phase–space averages by Birkhoff sums over trajectories. In measure–preserving systems, the two coincide almost everywhere (with respect to the initial point of a trajectory, under the invariant measure).

More precisely, if we consider an intensity function $H(z)$ (for simplicity of notation, from now on we write $H$ for ${\cal H}$, see Eq. (1)), we construct its distribution function $F_{H}(h)$ as follows. We start from a random point $z^{*}\in Z$. We iterate the dynamics $\varphi$ for a number of iterations $n_{0}$ sufficiently large to let the system reach a sort of dynamical equilibrium. At that moment we let $z_{0}=\varphi^{n_{0}}(z^{*})$ be the initial point of a trajectory and we compute the sample values $\{H(\varphi^{n}(z_{0})),n=0,\ldots,L-1\}$, by which the distribution function $F_{H}$ is

[TABLE]

Clearly, in numerical experiments the limit is replaced by the value obtained for a large but finite sample length $L$. We shall also use the complementary distribution function $E_{H}(h)=1-F_{H}(h)$: this is the measure of the tail of the distribution of values of $H$:

[TABLE]

The study of coupled map lattices has revealed puzzling phenomena like stable motion with exponentially long chaotic transients [25, 28], so that the numerical procedure just explained might be theoretically flawed. Nonetheless, the result of [8] on small networks ($N=128$) revealing a transition to positive Lyapunov exponents lead us to believe that convergence in Eq. (7) can be readily obtained.

3.2 Firing cascades as trees of causal excitations

A crucial difference between conventional coupled map lattices and the network under study is that in this latter interaction between different maps happens only at specific instants of time, while for most of the time individual systems evolve independently of each other. When it comes to physiological networks, the first and foremost consequence of this fact is that information can be stored, elaborated and transmitted [20]. In a dynamical perspective this calls attention to collective behaviors like synchronization [21, 29]. Here we focus on a more elusive set of phenomena: neuronal avalanches, which we prefer to call firing cascades.

Most of the studies on the collective behavior of neural networks define avalanches as the uninterrupted activity of the network during a certain time interval [9, 13]. Statistical properties of these collective events have been extensively studied in a large series of works, both experimental, see e.g. [9, 10, 11, 12, 13], and theoretical, with the aid of model systems [22, 23, 27, 30, 19, 31, 32, 33, 18, 34], but statistical analysis alone does not clarify the dynamical detail of the activity network and, in particular, of the impulse propagation. Moreover, it suffers from a degree of arbitrariness in the definition of avalanche. In fact, coarse–graining of time is usually employed, for spiking events appear at particular instants of time and must be grouped to define the quiescence periods separating different avalanches. Figure 6 shows a conventional raster plot of spikes, depicted as red points at coordinates $(n,i)$ where $n$ is the time at which neuron $i$ spikes. One observes regions with different density of spikes. The usual procedure [9] is to partition the time axis into different bins, so that a sequence of adjacent, non–empty bins constitutes a cascade. The arbitrariness in the choice of the bin length is apparent, and its influence on the resulting statistical properties has been well described in [9]. Moreover, different regions of the network even if uncorrelated but spiking at close instants may be classified in the same cascade. In short, the true dynamics of the network is overlooked in this approach. This is justified in the analysis of experimental data, when no access to the network of connections is available.

We pursue here a different approach, based upon the detailed knowledge of the dynamics presented in the previous sections. We have remarked that different types of neuron are present in the network: in particular, excitatory ISN are the origin of propagating excitations. When one of these neuron spikes, say neuron $j$, it acts upon its post-synaptic neurons via Eq. (6). As observed in Section 2, if such postsynaptic neuron is also an ISN, this action will alter its dynamics, anticipating or delaying its next spike. Otherwise, if the post-synaptic neuron $i$ is of the kind NISN, this excitation might eventually drive it to spike, this occurrence following from the dynamics generated by Eq. (2). The induced spiking, when it takes place, is not instantaneous: according to the coordinates $(x_{n}^{i},y_{n}^{i})$ at the time of excitation, the motion requires a certain amount of time to exit the region close to the fixed point and spike. It is nonetheless possible to predict if this event will take place, using the notion of basin of attraction of the spiking event discussed in Sect. 2. Recall that this is defined as the set of initial conditions of the isolated non–spiking neuron which lead to a spike, before returning to approach the fixed point. When $x_{n}$ is pushed by a synaptic excitations into this basin (red in Fig. 4) the NISN will fire, regardless of further excitatory actions: such spike can only be avoided by inhibitory actions.

This reasoning permits to draw a causal relationship between the spiking of a NISN and a single originating excitation from another neuron. Even if it is clear that in a complex network many pre–synaptic neurons can contribute to set the dynamics of $i$ on course to firing, we stipulate to ascribe the cause of the spike of the non–intrinsically spiking neuron $i$ to a single excitatory neuron, $i^{\prime}$, the one which is the last to act on $i$, at time $n^{\prime}$, before it enters the basin of attraction of the spiking event, to eventually fire at time $n$.

Excitations following this time are not considered in this count, because they are not necessary to induce a spike, which would happen nonetheless. Equally not considered are impulses preceding such precipitating event—the straw that broke the camel’s back—because if taken alone they might not be enough to cause firing. As a simpler rule, that can also be applied in experiments, we approximate the basin of attraction by the fact that $x^{i}$ exceeds a threshold value $x_{th}$. This procedure can be formulated as follows: we say that the spike of the excitatory neuron $i^{\prime}$ at time $n^{\prime}$ causes the spike of the INSN $i$ at time $n$, and we write $(n^{\prime},i^{\prime})\rightarrow(n,i)$ if

[TABLE]

The event that two neurons spike simultaneously and excite the same postsynaptic neuron is fully negligible and can be disposed by an ad–hoc prescription. Following rule (9) we can therefore create trees of excitations, which start at the firing of a single excitatory ISN and propagate to INSN’s through the network. Any such tree is defined to be a firing cascade. Return for illustration to Fig. 6, which also shows a realization of this construction. Green lines $(n^{\prime},i^{\prime})\rightarrow(n,i)$ join causally related spikes (blue symbols), on the background of red dots that mark all spiking events. The picture shows how the initial firing of an excitatory, intrinsically spiking neuron ($i^{\prime}=7$ at $n^{\prime}=811$) lead non–intrinsically spiking neurons to fire ($i=772$ at $n=820$, $i=235$ and $i=869$ at $n=823$, $i=2032$ at $n=827$) that, on their turn, stimulate additional NISN’s (for a total of 31 spikes in the figure). Clearly, not all physical connections are dynamically activated following a spike, for otherwise the cascade will have no end. In the example, neuron $i^{\prime}=7$ is also presynaptic to $i=908$, $i=1115$ and $i=2677$ that remain silent.

4 Statistical analysis of firing cascades

The dynamical approach described in the previous section permits to define precisely and measure accurately various quantities usually defined via time binning. These are:

• •

$S$, the size of a cascade, that is, the total number of spikes in a cascade—equally, the number of nodes in a causal tree;

• •

$T$, the time span of a cascade: the length of the time interval between the first spike and the last;

• •

$G$, the generation number, which is the number of levels of the tree—also, the length of the longest chain of successive excitations.

We are so equipped with the tools to analyze the statistics of cascades in the neuronal network: the function $H$ under considerations will be either one of $S,T,G$. As mentioned in Sect. 3.1, we replace phase–space averages by Birkhoff sums. We consider time-segments of length $L=4\cdot 10^{6}$ iterations, following an initial transient (discarded) consisting of 6000 iterations. Occasionally we raise this number to account for the fact that long transient can be observed in extended systems [25]. Using a parallel architecture, we accumulate the results of $32$ realizations of the initial conditions. We performed an extensive investigation, of which we report here a typical instance.

Example 1

$N=3000$* neurons, with $n_{ex}=2400$, $m_{ex}=2$, $n_{in}=600$, $m_{in}=1$, $p_{ex}=4$, $p_{in}=2$, Physical parameters are: $\beta=0.133$, $\alpha=3.6$, $\sigma_{\mbox{\tiny NINS}}=0.09$, $\sigma_{\mbox{\tiny ISN}}=0.103$, $\mu=10^{-3}$, $\psi=3$, $\chi_{\mbox{\small exc}}=0$, $\chi_{\mbox{\small inh}}=-1.1$.*

Since this model neither includes pruning nor adaptation, to achieve criticality one needs to tune the parameters of the network. Following [8], we vary the coupling intensity $W$ in Eq. (6). We first compute the average number of events per-unit time $\tau^{-1}$, observing a rather sharp transition as we increase $W$: when $W=0.084$ one has $\tau^{-1}=3.6\cdot 10^{-2}$, which increases to $\tau^{-1}=5.68\cdot 10^{-1}$ for $W=0.087$ and to $\tau^{-1}=5.65$ for $W=0.09$. At the largest coupling, roughly five neurons (out of the total 3000) spike on average at each single (discrete) time. This transition is clearly reflected in the statistics: in Fig. 7 we plot the fraction $\rho_{S}(S)$ of firing cascades with $S$ spikes and its (complementary) distribution function $E_{S}(S)$ defined in Eq. (8). Two values of the coupling constant $W$ yield the data in Fig. 7: $W=0.087$ and $W=0.09$. The former is visibly an under-threshold case, the latter displays a full-fledged cascade which extends over orders of magnitude. In the first case the distribution $E_{S}(S)$ is well fitted by a stretched exponential: $E_{S}(S)\simeq\exp\{-AS^{B}\}$ with $A\simeq-.1694$, $B\simeq.497$ (red curve). In the second, the power-laws $\rho_{S}(S)\simeq C^{\prime}S^{-\beta_{S}}$ and $E_{S}(S)\simeq CS^{-\beta_{S}+1}$ are well satisfied, with exponent $\beta_{S}=1.5$ ($C=0.75$, green line).

In the case of large coupling, $W=0.09$, we plot in Fig. 8 the statistical distribution $\rho_{G}(G)$ of generation levels $G$ and of time span $T$, $\rho_{T}(T)$. For both of these quantities we observe a power–law decay with exponents $\beta_{G}$ and $\beta_{T}$ close to two. In the same figure we also plot the inverse, for graphical convenience, of the quantity $\bar{S}(G)$, the average number of spikes in a cascade with $G$ levels, which approximately follows a power–law with the same exponent: $\bar{S}(G)\sim G^{\gamma}$, $\gamma\simeq 2$.

The observed exponents $\beta_{S}=3/2$ and $\beta_{T}=2$ are to a good approximation equal to those appearing in physiological experiments [9, 13] and in various models of neuronal networks [35, 30]. To the contrary, these values differ from those obtained in [8] for the system that we are investigating using the conventional definition of cascade. It is also to be quoted the measurement of these exponents in high resolution experimental data [12], still within the time–binning approach.

Critical scaling theory predicts [36, 12] that $\beta_{T}-1=(\beta_{S}-1){\gamma^{\prime}}$, where $\gamma^{\prime}$ is defined in terms of the scaling of the average of $S$ at fixed time length $T$, versus $T$. The exponents derived in [8] verify this scaling relation. Since in our model the time span of a cascade $T$ and its generation level $G$ are proportional, we can identify the exponents $\gamma$ and $\gamma^{\prime}$, providing the relation

[TABLE]

From the data of Fig. 8 we see that approximately $\gamma=2$, so that the relation (10) for the triple $\beta_{S}=3/2$, $\beta_{G}=\gamma=2$ is rather well satisfied by the numerical data.

5 Statistics of extreme events in a finite time–interval

The excitation of an intrinsically spiking neuron propagates sometimes to a large number of neurons and generates a cascade. In the previous section we have studied the relative frequency of these cascades versus their intensity, in a theoretically infinite time–interval. It is now interesting to consider the probability that none of these cascades, within a certain finite observational interval, exceeds a certain intensity threshold—equivalently, the probability that the largest cascade in a time–interval is smaller than such threshold. This investigation belongs to the realm of extreme value theory and goes under the name of statistics of block maxima.

Let in fact $n$ be the time at which an excitatory ISN fires and generates a cascade (possibly, and frequently, limited to such single spike). To this value we associate an intensity, $I(n)$, which can be chosen to be the size of the cascade $S$, its level $G$, or its time span $T$. If more than a cascade is originated at time $n$, $I(n)$ is the maximum of the values corresponding to such cascades. To the contrary, if no cascade starts at time $n$ (i.e. no excitatory ISN spikes at time $n$) $I(n)$ is null. Since the system is deterministic, $I(n)$ is actually also a function of the phase–space point $z$ at which the trajectory is located when one starts recording these intensities:

[TABLE]

Following standard usage we then define the block maxima function ${H}_{L}(z)$, for various values of $L$, as the largest intensity of a cascade that originates in the time–window $n=0,\ldots,L-1$:

[TABLE]

The cascade achieving maximum intensity is therefore appropriately called the extreme event in the time window under consideration; the wider is the time–window, the larger is presumably its intensity. An Extreme Value Law for ${H}_{L}(z)$ is an asymptotic form, for large $L$, of the statistical distribution of this variable. We now show that it applies to the present case.

Let us consider the distribution of the values of ${H}_{L}(z)$ with respect to an invariant measure $\mu$. As discussed before, we assume that the distribution function $F_{{H}_{L}}(h)=\mu(\{z\in Z\mbox{ s.t. }H_{L}(z)\leq h\})$ can be obtained numerically by Birkhoff time averages over trajectories of the system.

Suppose that different cascades arise independently of each other, that is, correlation of the values of $H(\varphi^{j}(z))$ for different $j$ decays quickly when their time difference increases. Under this assumption the probability that ${H}_{L}(z)$ is less than a certain value $h$, which is clearly equal to the probability that the intensity $H(\varphi^{n}(z))$ of all cascades in the time interval $n\in[0,L-1]$ is less than $h$, becomes a product of the probabilities of these individual events. If we let $\lambda$ be the density of cascades per unit time (not to be confused with the density of spikes per unit time), we expect a number $\lambda L$ of cascades to originate in a time–window of length $L$, so that

[TABLE]

Next, taking logarithms and letting $h$ grow, so that $F_{H}(h)$ approaches one, we obtain

[TABLE]

Clearly, this is a heuristic estimate, but it yields the same results predicted by rigorous theory, which requires appropriate conditions [1] to dispose of the correlations of the motion that might render the independence assumption incorrect. We can verify a posteriori the fact that these correlations fade in the limit of large $L$ by a direct evaluation of Eq. (14), when the asymptotic relation symbol is replaced by a much more stringent equality.

To do this, let us consider again the case of Example 1 when $W=0.09$. Dividing the total number of cascades by the total time span of the computation (averaged over $32$ different runs) yields the value $\lambda=8\cdot 10^{-3}$, which is much smaller than $\tau^{-1}=5.65$ and reveals that spikes are organized in cascades of average size $(\lambda\tau)^{-1}$, that is, roughly $22$ spikes per cascade, in this case. In Fig. 9 we plot the quantity $-\log(F_{{H}_{L}}(h))/L$ versus $h$, for various values of $L$, for the three functions $H=G,S,T$. In all three cases data points collapse on a curve, coinciding with $\lambda E_{H}(h)$, as predicted by Eq. (14). Observe that the value of $\lambda$ employed in the figure does not come from a fit of the data, but it is the one just obtained by counting the number of cascades. The extreme value law given by Eq. (14) appears therefore to apply, at least in this regime.

It is to be noted that in certain cases an additional positive factor $\theta$, smaller than one, appears in front of $\lambda E_{H}(h)$. This factor is called the extremal index: it implies a slower decay of $F_{{H}_{L}}(h)$, originating from the tendency of extreme events to appear in clusters. Rigorous theory has studied to a large detail this phenomenon in simple systems [37], which has also been detected in many dimensional systems [38]. The data of Fig. 9 indicate a unit value of this parameter.

In the critical regime the three quantities $H=G,S,T$ are characterized by a power–law decay: $E_{H}(h)\simeq A_{H}h^{-\beta_{H}}$, where $A_{H}$ are suitable constants and the exponents $\beta_{H}$ have been discussed in the previous section. Inserting this information in Eq. (14) yields the relation

[TABLE]

In Fig. 10 we plot the numerically determined function $F_{{H}_{L}}(h)$ when $H=S$ for values of $L$ and $h$ ranging over orders of magnitude and their difference with the analytical relation (15). This latter is remarkably small in all the range of data considered.

We can perform the same analysis in the subcritical case, Example 1 and $W=0.087$. In Fig. 11, which is the analogue of Fig. 9 for this case, we plot again $-\log(F_{{H}_{L}}(h))/L$ and $\lambda E_{H}(h)$, for $H=G,S,T$. Collapse of the data points is again observed, but for larger values of $h$ than in the previous case. The asymptotic law for the distribution of the variable $H_{L}(h)$ becomes here

[TABLE]

The experimental function $F_{{H}_{L}}(h)$ for $H=S$ is plotted versus $L$ and $h$ in the left frame of Fig. 12. Difference with the theoretical expression (16) is plotted in the right frame. As expected from Fig. 11, convergence improves as $L$ and $h$ increase.

Within the framework of Extreme Value Theory, the two cases just discussed differ because they imply different normalizing sequences $h_{L}(t)$, which we now define, while providing the same limiting distribution $F(t)=e^{-t}$. In fact, observe that Eqs. (13,14) imply that

[TABLE]

Let $t$ be a positive, real variable. The sequence of levels $h_{L}(t)$, $L=1,2,\ldots$, parameterized by $t$, is defined so to satisfy the limit behavior

[TABLE]

As a consequence of the above,

[TABLE]

The two previous equations are the core results of EVT, when added to the fact that three different limiting distributions $F(t)$ can be achieved after further rescaling, as we do below in the critical case. They show that $h_{L}(t)$ is a quantile of the asymptotic distribution of the statistical variable $H_{L}(z)$, whose analytical form involves the ratio $t/L$ and the inverse of the complementary distribution of the function $H$:

[TABLE]

In the subcritical case these quantiles grow as the power of a logarithm:

[TABLE]

At criticality a power–law behavior is present:

[TABLE]

and any $q$ moment of the statistical distribution of $H_{L}(z)$ scales as $(\frac{\lambda L}{t})^{\frac{q}{\beta_{H}}}$. In this latter case, the previous relations can be put in one of the three standard forms [2]. In fact, simple algebra permits to rewrite

[TABLE]

as $L\to\infty$, which is the expression of a Fréchet distribution.

6 Conclusions

We have described in this paper an extension of the dynamical theory of extreme events to time–evolving phenomena in complex, multi–dimensional systems. We have applied it with success to a model of neuronal network, that was studied in full detail. This analysis attains wide generality because this example is characterized by features that are frequently encountered: it is composed of many almost identical subsystems, linked in a sparse, directed network; each individual system integrates the input from the network and in turn returns a signal to it; signal propagation can develop a full chain of excitations. This has lead us to introduce a novel approach to define neuronal avalanches, whose statistical characteristics were studied within the framework of extreme value theory. The same definition of avalanche and a similar analysis can be obviously applied to other systems in the class just described.

The presented results are in line with previous theoretical models and experimental findings on neuronal networks. Nonetheless, since the focus of this work is on dynamical theory, we did not investigate various phenomena of physiological relevance, such as the balance between excitatory and inhibitory neurons (that are nonetheless present in our model) and the rôle of pruning and adaptation to sustain criticality [34]. Also, we did not perform rescaling of the shape of cascades at criticality, to reveal a single curve, as in crackling noise, although it is highly likely that this will happen, as shown in [8]. A further topic that could be addressed is the computation of the branching ratio in the chain of excitations, which is clearly possible with the aid of our definition of cascade. And finally, we did not compare among them different network topologies.

While the set-up that we propose to extend the dynamical extreme value theory to complex phenomena, evolving in time, is fully rigorous, our results are formal, even if confirmed by numerical experiments. This poses a new challenge for rigorous mathematical analysis, while at the same time it offers a paradigm for experimental studies in various physical systems.

Acknowledgements Numerical computations for this paper were performed on the Zefiro cluster at the INFN computer center in Pisa. Enrico Mazzoni is warmly thanked for assistance.

7 Appendix: stability of the fixed point of Rulkov map

It is easy to derive the fixed point of Rulkov map and its stability. Consider again an isolated neuron, and assume that we are in the first alternative of Eq. (3), which yields the map

[TABLE]

From this, the unique fixed point $(\bar{x},\bar{y})$ easily follows:

[TABLE]

This result is consistent with Eq. (3) provided that $\bar{x}<0$, that is, $\sigma<1$.

The differential of the transformation, i.e. the linear map from the tangent space at the fixed point to itself is

[TABLE]

Putting for simplicity $A(\alpha,\sigma)=\frac{\alpha}{(2-\sigma)^{2}}$ the linear map $J$ has real eigenvalues provided that $A(1-(\alpha,\sigma))^{2}\geq 4\mu$. In the opposite case, complex conjugate pairs appear. In this case, the square modulus of the eigenvalues equals the determinant of $J$ and is

[TABLE]

This situation is pictured in Fig. 13: $A(\alpha,\sigma)$ and $\mu$ are always positive and are taken as coordinates in the $(A,\mu)$ plane. When the corresponding point lies below the parabola $(1-A(\alpha,\sigma))^{2}=4\mu$, $J$ has two real eigenvalues. If moreover $A<1$ (region I) their absolute value is less than one and the (hyperbolic) fixed point is attractive, while it is repulsive when $A>1$ (region IV). In the region above the parabola and below the line $A(\alpha,\sigma)+\mu=1$ (region II) $J$ has a pair of complex conjugate eigenvalues of modulus less than one, so that the fixed point is still attractive, and trajectories spiral down to it. When finally $A(\alpha,\sigma)$ is above both the parabola and the straight line (region III) two eigenvalues of modulus larger than one appear: the fixed point is unstable, and trajectories spiral out of it.