text:
# 1 Introduction

## 1 Introduction

That evolved brains are highly sensitive organs is an everyday observation. Viewed as a dynamical system, a brain may be said to be unusually susceptible to perturbations and initial conditions. This leads one to ask whether brains may be operating near some form of instability, or criticality, a hypothesis related to the notions of computation at the edge of chaos (Langton 1990) and self-organized criticality (Bak et al. 1987). In this paper we propose that while most regulation mechanisms at work in the brain act according to a classical homeostasis schema, i.e., have a stabilizing effect, an opposite effect could result from the regulation of synaptic weights by a specific form of Hebbian covariance plasticity. Such a regulation may bring the system near criticality. We suggest that regulated criticality may be the mechanism whereby sensitivity is maintained throughout life in the face of ongoing changes in brain connectivity.

Hebbian synaptic plasticity (Hebb 1949) plays an important role in the development of the nervous system, and is also believed to underlie many instances of learning in the adult. A covariance rule of Hebbian plasticity roughly states that the change in the efficacy of a given synapse varies in proportion to the covariance between the presynaptic and postsynaptic activities. As noted by many authors (e.g. Sejnowski 1977a, 1977b; Bienenstock et al. 1982; Linsker 1986; Sejnowski et al. 1988), a covariance-type rule is preferable to a rule that uses the mere product of pre- and post-synaptic activities because the covariance rule predicts not only weight increases but also activity-related weight decreases, and as a consequence allows convergence to non-trivial connectivity states. Some forms of covariance plasticity have been shown to be optimal for information storage (Willshaw and Dayan 1990; Dayan and Willshaw 1991; Dayan and Sejnowski 1993). Also, evidence for Hebbian plasticity of the covariance type has been reported in several preparations (Frégnac et al. 1988, 1992; Stanton and Sejnowski 1989; Artola et al. 1990; Dudek and Bear 1992).

We shall investigate, in a simple network including excitatory and inhibitory neurons, the effect of covariance plasticity acting as a mechanism of regulation, rather than supervised learning. Synaptic modification results in changes—quantitative or qualitative—in the activity that reverberates in the network; these changes in turn cause further modification of the weights, thereby creating a feedback loop between activity and connectivity. Studying this loop as such, i.e., independently from any input and output, we demonstrate that, under fairly general conditions, it causes the network to converge to a critical surface in parameter space, the locus of an abrupt transition between different activity modes. In Metzger and Lehmann (1990, 1994) a similar Hebbian rule has been studied in the context of supervised learning of temporal sequences.

Schematically, the convergence to a critical state can be explained as follows. One mode of behavior of a network including excitatory and inhibitory neurons is oscillation; such behavior takes place if the synaptic weights linking excitatory neurons to each other—we will refer to these as E-to-E weights—are high enough but not too high. Oscillation entails high covariance values, hence, according to the covariance rule, results in further increase of the E-to-E weights. If however these weights are higher than a certain critical value—which depends on other parameters of the system—oscillatory behavior is impossible, hence covariance is low or zero, hence, in accord with the covariance rule used, the E-to-E weights decrease. As a result, the E-to-E weights stabilize around the critical surface that separates the region of oscillation from the region(s) of steady firing.

Our study is conducted in the simplest type of network that will support oscillatory activity: all synaptic weights of a given type are given identical values, and so are all firing thresholds of a given type. This results in a system with just six parameters—four synaptic weights and two thresholds—and a limited range of behaviors. Essentially, all neurons fire uniformly, either at a constant rate (the number of possible rates of firing is one or two, depending on parameters) or periodically in time. In the thermodynamic, i.e., large-size, limit, the dynamics of the network is adequately described by a system of differential equations obtained through a classical mean-field approximation.

We first perform a simple bifurcation analysis of this differential system (Guckenheimer and Holmes 1983). We then show that the effect of covariance regulation is to stabilize the parameter state at a surface of transition, where the dynamics exhibits an instability. Such a critical parameter state for a dynamical system may be characterized as degenerate. A generic, i.e., non-exceptional, state is one where one expects to find the system in the absence of further assumptions. Mathematically, a generic state of a dynamical system is in the interior of a parameter region corresponding to a given behavior, and the system in such a parameter state is said to be structurally stable; the set of non-generic parameter states has measure zero. We shall show that a state of higher degeneracy, characterized as a point of intersection of several critical surfaces, can be achieved by the simultaneous regulation of several parameters. In the vicinity of that highly degenerate state, the system displays a range of behaviors, including chaos.

The plan of the paper is as follows. In the next section we study the dynamical properties of our simple network—in the differential-equation formulation—with fixed parameters (synaptic weights and firing thresholds). We characterize the bifurcations which take place at the boundaries between domains corresponding to different modes of behavior. This study is conducted for a reduced system, where the thresholds are eliminated in such a way as to render the dynamics symmetric about the origin. Section 3 describes the regulation equations. Section 4 describes the behavior of the regulated reduced system, and Section 5 that of the regulated full system.

## 2 The fixed-parameter model

This section describes the dynamics of the model with fixed parameters. We first briefly describe a network consisting of a large number ($`2N`$) of binary-valued neurons operating under a stochastic dynamics. However, rather than using this network for our study of plasticity, we make a number of simplifications and approximations, leading to a deterministic two-variable differential system with just six parameters. The two variables are the excitatory and inhibitory population averages of cell activity in the $`2N`$-dimensional model; the six parameters include the four average weights of the synapses within and between these two populations, as well as the average firing thresholds for the two populations. We then study the asymptotic behavior of this differential system for various parameter values. Different types of asymptotic behavior, in different regions of the parameter space, correspond to different phases of the stochastic system, and we pay particular attention to the bifurcations of the solutions, where the bifurcation parameters are the synaptic weights—see Schuster and Wagner (1990) and Borisyuk and Kirillov (1992) for a related bifurcation analysis. Bifurcations correspond to phase transitions in the statistical-physics formulation (the original $`2N`$-dimensional model).

We consider a fully-connected network of $`N`$ excitatory and $`N`$ inhibitory linear-sigmoidal $`\{0,1\}`$-valued neurons,111It is not essential that the numbers of excitatory and inhibitory neurons be the same. operating under a stochastic dynamics. We denote the activity of the $`i`$-th excitatory, resp. inhibitory, neuron by $`x_i^\text{E}(t)`$, resp. $`x_i^\text{I}(t)`$, with $`x_i^\text{E}(t),x_i^\text{I}(t)\{0,1\}`$, $`i=1\mathrm{}N`$, and we denote the synaptic weights by $`w_{ij}^{\text{EE}},w_{ij}^{\text{EI}},w_{ij}^{\text{IE}},w_{ij}^{\text{II}}`$, $`i,j=1\mathrm{}N`$, where $`i`$ is postsynaptic and $`j`$ presynaptic, and the superscripts indicate, for each of the two neurons, whether it is excitatory or inhibitory. Thus, for all $`i`$ and $`j`$, $`w_{ij}^{\text{EE}}`$ and $`w_{ij}^{\text{IE}}`$ are positive or zero, whereas $`w_{ij}^{\text{EI}}`$ and $`w_{ij}^{\text{II}}`$ are negative or zero.

The local field on excitatory neuron $`i`$, i.e., the difference between its membrane potential and its firing threshold $`h_i^\text{E}`$, is $`g_i^\text{E}(t)=_jw_{ij}^{\text{EE}}x_j^\text{E}(t)+_jw_{ij}^{\text{EI}}x_j^\text{I}(t)h_i^\text{E}`$. Similarly, the local field on inhibitory neuron $`j`$ is $`g_i^\text{I}(t)=_jw_{ij}^{\text{IE}}x_j^\text{E}(t)+_jw_{ij}^{\text{II}}x_j^\text{I}(t)h_i^\text{I}`$, where $`h_i^\text{I}`$ is the threshold of inhibitory neuron $`i`$. The network dynamics is defined by: (i) selecting at random, with uniform probability, one of the $`2N`$ neurons; (ii) computing its local field $`g(t)`$, of the form $`g_i^\text{E}(t)`$ or $`g_i^\text{I}(t)`$; and (iii) defining the state of the network at time $`t+\delta t`$ to be equal to the state at time $`t`$ except, possibly, for the selected neuron, whose state becomes—or stays—1 with probability $`\frac{1}{2}(1+\mathrm{tanh}(\beta g(t)))`$. Parameter $`\beta `$ is a fixed non-negative number, an inverse temperature. The temperature $`T=1/\beta `$ measures the amount of noise in the system: the higher the temperature, the noisier the dynamics. The update interval is $`\delta t=1/(2N)`$, so that each neuron is updated on average once every time unit. This asynchronous dynamics, of the Glauber type (Glauber 1963), is widely used in statistical-mechanics models; it lends itself to a convenient mean-field approximation (see below).

A system such as the one just described will exhibit a highly diverse range of behaviors, depending on the values of the synaptic weights and firing thresholds. But we now make the much simplifying assumption that synaptic weights and firing thresholds are uniform across each class. Specifically, for all $`i,j=1,\mathrm{},N`$, we assume that $`h_i^\text{E}=h^\text{E}`$, $`h_i^\text{I}=h^\text{I}`$, $`w_{ij}^{\text{EE}}=w^{\text{EE}}/N`$, $`w_{ij}^{\text{EI}}=w^{\text{EI}}/N`$, $`w_{ij}^{\text{IE}}=w^{\text{IE}}/N`$, and $`w_{ij}^{\text{II}}=w^{\text{II}}/N`$, where $`h^\text{E}`$, $`h^\text{I}`$, $`w^{\text{EE}}`$, $`w^{\text{EI}}`$, $`w^{\text{IE}}`$ and $`w^{\text{II}}`$ are fixed parameters, and $`w^{\text{EE}}`$, $`w^{\text{EI}}`$, $`w^{\text{IE}}`$ and $`w^{\text{II}}`$ are non-negative. The dynamics is thus parameterized by six constants, four synaptic weights and two thresholds; $`\beta `$ is a mere multiplicative factor common to all six parameters, yet it is convenient to use it as a seventh parameter. Unless otherwise mentioned, $`\beta `$ will be 1.

Due to this uniformity assumption, all neurons in any of the two populations experience the same field at any given time. This system exhibits a limited number of fairly simple behaviors, of which Figure 1 is an example. This figure shows the time variation of $`x_i^\text{E}(t)`$ and $`x_i^\text{I}(t)`$, the average activation levels across the excitatory and inhibitory populations. In this example, parameters are: $`N=70`$, $`w^{\text{EE}}=12`$, $`w^{\text{IE}}=8`$, $`w^{\text{EI}}=10`$, $`w^{\text{II}}=2`$, $`h^\text{E}=1`$, $`h^\text{I}=3`$. One unit on the time axis corresponds to $`2N`$ updates, so that each neuron is updated, on average, once every time unit. For these parameter values, the system oscillates. Note that the oscillation is not perfectly regular, a finite-size effect. Note also that the inhibitory activity lags somewhat behind the excitatory activity: the excitatory neurons first trigger the inhibitory ones, which in turn extinguish, for a while, the excitatory population.

(Insert Figure 1 around here)

The presence of oscillations and the amplitude and shape of the waveform depend on the various parameters. However, rather than pursuing this study of the stochastic system, we shall consider the approximation that obtains in the thermodynamic limit, that is, when $`N\mathrm{}`$. The update interval $`\delta t=1/(2N)`$ then goes to 0 and so does each individual synaptic weight. Straightforward approximations (Rubin 1988; Schuster and Wagner 1990) then lead to a continuous-time differential system for the population averages of the excitatory and inhibitory activation levels, which we denote, respectively, by $`s`$ and $`\sigma `$:

$$\{\begin{array}{c}\dot{s}(t)=.5s(t)+.5\mathrm{tanh}[\beta (w^{\text{EE}}s(t)w^{\text{EI}}\sigma (t)h^\text{E})]\hfill \\ \dot{\sigma }(t)=.5\sigma (t)+.5\mathrm{tanh}[\beta (w^{\text{IE}}s(t)w^{\text{II}}\sigma (t)h^\text{I})].\hfill \end{array}$$
(1)
Note that the variables $`s(t)`$ and $`\sigma (t)`$ remain at all $`t`$ within the interval . When $`\beta =0`$ system 1 has a unique attractor, $`(s,\sigma )=(.5,.5)`$. Indeed, in the high-temperature limit, all neurons act independently of each other and fire with probability .5 at each time.

We shall now make a last simplification, whose purpose it is to render $`(.5,.5)`$ a fixed point—though not necessarily stable—at all temperatures and for all values of the synaptic weights. This is easily achieved by letting the thresholds $`h^\text{E}`$ and $`h^\text{I}`$ be determined by the synaptic weights as follows:

$$\begin{array}{c}h^\text{E}=.5(w^{\text{EE}}w^{\text{EI}})\hfill \\ h^\text{I}=.5(w^{\text{IE}}w^{\text{II}}).\hfill \end{array}$$
(2)
It is then convenient to adopt the change of variables: $`ss.5`$, $`\sigma \sigma .5`$, and system 1 becomes:

$$\{\begin{array}{c}\dot{s}(t)=s(t)+.5\mathrm{tanh}[\beta (w^{\text{EE}}s(t)w^{\text{EI}}\sigma (t))]\hfill \\ \dot{\sigma }(t)=\sigma (t)+.5\mathrm{tanh}[\beta (w^{\text{IE}}s(t)w^{\text{II}}\sigma (t))].\hfill \end{array}$$
(3)
In 3, the variables $`s`$ and $`\sigma `$ are in the interval $`[.5,+.5]`$, and the only parameters left are the four synaptic weights and the inverse temperature. For all parameter values, the origin is a fixed point of system 3. A different position for the fixed point could be obtained with an appropriate modification of equations 2, yet in the current version the fixed point is also a center of symmetry. For the moment, this hard-wired symmetry should be regarded as an ad-hoc device, whose purpose is to make the analysis more convenient. We shall refer to system 1 as the full system, and to system 3 as the reduced system. We shall see in Section 5 that under appropriate regulation the two systems behave similarly.

We now discuss some important properties of the reduced system, system 3 (see also Rubin 1988). Consider first Figure 2a (phase diagram), which shows four trajectories of the state $`(s(t),\sigma (t))`$; the starting points of these trajectories are indicated by triangles. The parameters (synaptic weights) used in this example are identical to those used in Figure 1, i.e., $`w^{\text{EE}}=12`$, $`w^{\text{IE}}=8`$, $`w^{\text{EI}}=10`$, $`w^{\text{II}}=2`$. As expected, the asymptotic behavior is periodic; there is a limit cycle which attracts all points of the square $`[.5,.5]^2`$, except the unstable equilibrium $`(0,0)`$. Motion is counterclockwise, for, as mentioned above, $`\sigma (t)`$ lags behind $`s(t)`$.

In addition to these four orbits, Figure 2a shows two curves, the $`s`$\- and $`\sigma `$-nullclines for system 3. These are the loci of the points $`(s,\sigma )`$ such that $`ds/dt`$, resp. $`d\sigma /dt`$, vanish. The equations for the $`s`$\- and $`\sigma `$-nullclines are easily seen to be, respectively:

$$\sigma =\frac{1}{w^{\text{EI}}}(w^{\text{EE}}sT\mathrm{tanh}^1(2s)),$$
(4)
$$s=\frac{1}{w^{\text{IE}}}(w^{\text{II}}\sigma +T\mathrm{tanh}^1(2\sigma )).$$
(5)
The $`\sigma `$-nullcline is an increasing sigmoid-shaped curve, whereas the $`s`$-nullcline generally has the shape of an ‘S’ lying on its side. Of particular interest are the intersection points of the two nullclines; these are the fixed points of the dynamics. In the case illustrated in Figure 2a, the only intersection is $`(0,0)`$, an unstable equilibrium. Trajectories intersect the $`s`$-, resp. $`\sigma `$-, nullcline in a direction parallel to the $`\sigma `$-, resp. $`s`$-, axis.

The study of the nullclines is of interest because it is often possible to predict how a parameter change will affect the dynamics of the system by reasoning about how the nullcline diagram will change; the bifurcation we shall be mostly interested in is associated with a conspicuous change in this diagram. Note that the $`s`$-nullcline is affected by parameters $`w^{\text{EE}}`$ and $`w^{\text{EI}}`$, whereas the $`\sigma `$-nullcline is affected by parameters $`w^{\text{II}}`$ and $`w^{\text{IE}}`$.

(Insert Figure 2 around here)

Let us consider first the changes brought about by letting parameter $`w^{\text{EE}}`$ grow, starting from the point $`w^{\text{EE}}=12`$ for which the system oscillates; other parameters are unchanged. When $`w^{\text{EE}}`$ grows, the slope of the central, quasi-linear, part of the $`s`$-nullcline increases (see equation 4); that part of the curve rotates about the symmetry center (0,0). As a result, the peak of the $`s`$-nullcline to the right approaches the upper part of the sigmoid-shaped $`\sigma `$-nullcline, and the minimum of the $`s`$-nullcline to the left approaches the lower part of the $`\sigma `$-nullcline. Eventually, at a certain critical value $`\widehat{w}_{\text{sn}}^{\text{EE}}`$ (subscript ‘sn’ stands for ‘saddlenode’—see below), the two curves become tangent to each other. This happens in two points at once, near the upper right-hand corner and near the lower left-hand corner, due to the symmetry of the system. This situation is depicted in Figure 2b: $`w^{\text{EE}}`$ is exactly equal to the critical value $`\widehat{w}_{\text{sn}}^{\text{EE}}`$ (with parameters as above, $`\widehat{w}_{\text{sn}}^{\text{EE}}14.22`$), and the nullclines are just tangent to each other.

When $`w^{\text{EE}}`$ grows a little further, each point of contact splits into two intersection points, of which one is an attractor. Figure 2c shows this situation, with $`w^{\text{EE}}=15`$, somewhat above the critical value $`\widehat{w}_{\text{sn}}^{\text{EE}}`$. Four trajectories are shown, in addition to the two nullclines. The system has five fixed points, three unstable ones and two stable ones (attractors). Only the stable fixed points are of interest to us; they are very near the upper right-hand and lower left-hand corners of the square, corresponding to high, respectively low, excitatory and inhibitory activities.

The bifurcation occurring at $`\widehat{w}_{\text{sn}}^{\text{EE}}`$ is of the saddlenode type. It results in a drastic change of behavior of the system: the periodic attractor disappears and is ‘siphoned’ into the two new point attractors. These two points attract the entire square, except a set of measure 0 which includes the three unstable fixed points. Thus, altough this bifurcation is caused by a mere local change, namely the intersection of the nullclines, it results in a reorganization of the dynamics that is both abrupt and global.222As mentioned, two distinct saddlenode bifurcations take place simultaneously. Such a double bifurcation is not generic; it occurs here due to the symmetry that we introduced when reducing system 1 into system 3.

Having described the breakdown of oscillations when parameter $`w^{\text{EE}}`$ is increased, we now consider the opposite change, that is, we let $`w^{\text{EE}}`$ decrease. This results in a decrease of the slope of the central, increasing, portion of the $`s`$-nullcline (equation 4). Eventually, the curve becomes monotonically decreasing; this does not alter the number of intersections of the nullclines, point $`(0,0)`$ remaining the sole equilibrium. However, the amplitude of the limit cycle decreases along with $`w^{\text{EE}}`$. The cycle eventually collapses to a point; the equilibrium $`(0,0)`$ has then become stable. This can be seen in a linear stability analysis of system 3 around point $`(0,0)`$. It is easily shown that, in case there are two complex conjugate eigenvalues,333The condition for this is $`4w^{\text{EI}}w^{\text{IE}}>(w^{\text{EE}}+w^{\text{II}})^2`$. the real part of these eigenvalues is negative if and only if $`w^{\text{EE}}<w^{\text{II}}+4T`$. Thus, $`w^{\text{II}}+4T`$ is a critical value for parameter $`w^{\text{EE}}`$. We define $`\widehat{w}_{\text{hopf}}^{\text{EE}}\stackrel{\mathrm{def}}{=}w^{\text{II}}+4T`$ (with the current parameter setting, $`\widehat{w}_{\text{hopf}}^{\text{EE}}=6`$). The change of behavior occurring at $`\widehat{w}_{\text{hopf}}^{\text{EE}}`$ is a normal444That is, supercritical. However, for very large values of $`w^{\text{IE}}`$, the bifurcation is subcritical—see footnote 5. Hopf bifurcation.

So far, we studied the behavior of system 3 for different values of parameter $`w^{\text{EE}}`$, all other parameters being fixed. In other words, we described the system’s behavior on a particular 1-dimensional subspace of the 4-dimensional parameter space. We now extend this study to a 2-dimensional subspace, the $`(w^{\text{EE}},w^{\text{IE}})`$ plane. Figure 3a is the bifurcation diagram of system 3 in that plane, with other parameters as before ($`w^{\text{EI}}=10`$, $`w^{\text{II}}=2`$). This diagram shows three distinct regions, corresponding to three different attractor configurations; unstable fixed points and unstable limit cycles are ignored in this diagram. In the middle region—which we call region $`𝒫`$, for Periodic—the system oscillates. The boundary of this region to the right is the saddlenode bifurcation curve, which we denote $`𝒮`$; as discussed above, the rightmost region has two point attractors, and we call it region $`𝒯`$. The leftmost region, which we call $`𝒪`$, has only one point attractor, the center of symmetry $`(0,0)`$; it is separated from region $`𝒫`$ by the Hopf bifurcation curve, a vertical line of equation $`w^{\text{EE}}=\widehat{w}_{\text{hopf}}^{\text{EE}}`$. The curve in the lower left of the diagram, separating region $`𝒪`$ from region $`𝒯`$, is the locus of a pitchfork bifurcation. This bifurcation diagram, obtained for one particular set of values of the parameters $`w^{\text{EI}},w^{\text{II}}`$ and $`\beta `$, is representative of the general case.555It is however simplified in two ways. First, the transition from region $`𝒫`$ to region $`𝒯`$ is of the saddlenode type only for large enough values of $`w^{\text{IE}}`$; this range of values corresponds roughly to the straight portion of curve $`𝒮`$ (Figure 3a). To see why this is so, consider again Figure 2b, the nullcline diagram at the bifurcation, with $`w^{\text{IE}}=8`$. Note that the points of contact between the nullclines appear near the corners of the square, far from the origin; this is due to the fact that $`w^{\text{IE}}`$ is relatively large, hence the slope of the $`\sigma `$-nullcline at the origin is larger than the slope of the $`s`$-nullcline. The bifurcation is then of the saddlenode type, as described. If however $`w^{\text{IE}}`$ is small, hence so is the slope of the $`\sigma `$-nullcline at the origin, the transition from $`𝒫`$ to $`𝒯`$ as $`w^{\text{EE}}`$ is increased takes place differently. A pair of intersection points between the nullclines first split off from the origin; these are unstable equilibria. As $`w^{\text{EE}}`$ increases, these two equilibria move away from the origin, while remaining inside the large stable limit cycle. At a certain critical value for $`w^{\text{EE}}`$ they become stable—a (double) subcritical Hopf bifurcation—and almost immediately thereafter the large limit cycle disappears. Thus, the transition from region $`𝒫`$ to region $`𝒯`$ really takes place in two steps, giving rise to a three-attractor behavior: the system has one large limit-cycle attractor as well as two point attractors, the latter being inside the cycle. The region of the $`(w^{\text{EE}},w^{\text{IE}})`$ plane where this behavior takes place is a strip extending along the lower, curved, part of the $`𝒫/𝒯`$ boundary; it is too narrow to be seen in Figure 3a. (With parameters $`w^{\text{EI}}`$ and $`w^{\text{II}}`$ as above and $`w^{\text{IE}}=2.75`$, the three-attractor behavior occurs for $`w^{\text{EE}}`$ between $`8.993`$ and $`9.030`$. For some other values of $`w^{\text{EI}}`$ and $`w^{\text{II}}`$ this behavior does not occur at all, and the transition from $`𝒫`$ to $`𝒯`$ is always of the saddlenode type.) For the purpose of this paper (see footnotes 8 and 11) it is important to note that the point attractors appear either exactly or almost at the same time as the periodic attractor disappears. The second approximation in the bifurcation diagram, mentioned only for the sake of completeness, concerns the $`𝒪`$-to-$`𝒫`$ transition. This is generally a smooth, supercritical, Hopf bifurcation. However, as mentioned in footnote 4, this Hopf bifurcation becomes subcritical for very large values of $`w^{\text{IE}}`$. There is thus a narrow region to the left of the bifurcation line $`w^{\text{EE}}=\widehat{w}_{\text{hopf}}^{\text{EE}}`$ where the limit-cycle attractor coexists with the point attractor (0,0); for instance, at $`w^{\text{IE}}=100`$, the width of this region is $`0.63`$.

(Insert Figure 3 around here)

In sum, the $`(w^{\text{EE}},w^{\text{IE}})`$ bifurcation diagram for system 3 is characterized by a central periodic-attractor region, a large vertical patch extending to $`+\mathrm{}`$ in the $`w^{\text{IE}}`$ direction (phase $`𝒫`$), flanked by point-attractor regions on each side (phases $`𝒪`$ and $`𝒯`$). The transition from $`𝒫`$ to $`𝒯`$ is abrupt ($`𝒮`$ line), while the transition from $`𝒪`$ to $`𝒫`$ is smooth. As mentioned in the Introduction, system 1—the full system—is not amenable to such a thorough analysis; however, we shall see in Section 5 that the two systems behave in much the same way under the plasticity rules that we shall now introduce.

## 3 The regulation equations

Whereas in the previous section the synaptic weights $`w^{\text{EE}}`$ and $`w^{\text{IE}}`$ were fixed parameters, they will now be made to evolve. Their evolution will obey a Hebbian covariance rule, hence be a function of second-order temporal averages of the dynamic variables $`s`$ and $`\sigma `$. Synaptic plasticity creates a regulation loop: changing the parameters affects the dynamics of the system, which in turn alters the second-order moments of $`s`$ and $`\sigma `$. Formally, the regulation is implemented by introducing additional differential equations, coupled to system 3 (or to system 1—see Section 5). The rate of change of $`w^{\text{EE}}`$ and $`w^{\text{IE}}`$ will typically be several orders of magnitude slower than that of $`s`$ and $`\sigma `$.

Let us first define, for any function of time $`r(t)`$, a moving time average:

$$\overline{r}(t)=\rho _{\mathrm{}}^tr(u)e^{\rho (ut)}𝑑u.$$

Parameter $`\rho `$ is a positive constant, physically an inverse time; the larger $`\rho `$, the narrower the averaging kernel. Equivalently, $`\overline{r}(t)`$ may be defined by a differential equation, more convenient for simulation purposes:

$$\frac{d\overline{r}(t)}{dt}=\rho (r(t)\overline{r}(t)).$$

Consider now, with reference to the original stochastic model (Section 2), the instantaneous covariance between two excitatory neurons $`i`$ and $`j`$, defined as: $`c_{ij}^{\text{EE}}(t)\stackrel{\mathrm{def}}{=}(x_i^\text{E}(t)\overline{x}_i^\text{E}(t))(x_j^\text{E}(t)\overline{x}_j^\text{E}(t))`$. If we take the population average $`c_{ij}^{\text{EE}}(t)`$ of this instantaneous covariance, we obtain, in the thermodynamic limit $`N\mathrm{}`$, the instantaneous variance of $`s(t)`$:

$$c^{\text{EE}}(t)\stackrel{\mathrm{def}}{=}(s(t)\overline{s}(t))^2.$$
(6)
It is this quantity $`c^{\text{EE}}`$ that we use to regulate the excitatory-to-excitatory synaptic weight $`w^{\text{EE}}`$. The regulation equation is linear in $`c^{\text{EE}}`$:

$$\frac{dw^{\text{EE}}(t)}{dt}=\epsilon ^{\text{EE}}(c^{\text{EE}}(t)\theta ^{\text{EE}}).$$
(7)
Parameters $`\epsilon ^{\text{EE}}`$ and $`\theta ^{\text{EE}}`$ are positive. Note that the quantity $`c^{\text{EE}}(t)`$ is always non-negative; the term $`\theta ^{\text{EE}}`$ is therefore necessary to allow for decreases of $`w^{\text{EE}}`$.

We shall also consider a regulation for $`w^{\text{IE}}`$, the synaptic weight from excitatory to inhibitory neurons, although this regulation will play a less important role than that of $`w^{\text{EE}}`$. The modification rule for $`w^{\text{IE}}`$ has the same form as equation 7, yet it uses the excitatory-to-inhibitory instantaneous covariance, defined as:

$$c^{\text{IE}}(t)\stackrel{\mathrm{def}}{=}(s(t)\overline{s}(t))(\sigma (t)\overline{\sigma }(t)).$$
(8)
The regulation equation for $`w^{\text{IE}}`$ then reads:

$$\frac{dw^{\text{IE}}(t)}{dt}=\epsilon ^{\text{IE}}(c^{\text{IE}}(t)\theta ^{\text{IE}}).$$
(9)
In equation 9, $`\theta ^{\text{IE}}`$ is a positive constant, as $`\theta ^{\text{EE}}`$ in equation 7. However, the modification rate constant $`\epsilon ^{\text{IE}}`$ is negative. The main reason for this will be given in the next section; for now, note that this choice is consistent with the spirit of Hebb’s principle, for, when considered postsynaptically to the target neuron, the effect of synapse reinforcement if that target neuron is inhibitory is the opposite of the effect obtained if the target neuron is excitatory.

## 4 Behavior of the regulated reduced system

This section describes the behavior of the regulated reduced system. We demonstrate that each of the two regulation loops introduced in Section 3, when acting separately, brings the system to the critical surface $`𝒮`$, the locus of an abrupt phase transition (saddlenode bifurcation). We then examine the behavior of the system with the two regulation loops active simultaneously; we show that under some conditions the state converges to a point on $`𝒮`$ with a remarkable nullcline configuration.

Before we consider the regulation proper, let us examine how the covariances change across the $`(w^{\text{EE}},w^{\text{IE}})`$ plane. Figure 3b shows the values of $`\overline{c}^{\text{EE}}`$, the time average of the instantaneous variance of $`s(t)`$,666This corresponds, in the original system, to the population- and time-average of the covariance, $`\overline{c}_{ij}^{\text{EE}}`$; the latter becomes $`\overline{c}^{\text{EE}}`$ in the thermodynamic limit $`N\mathrm{}`$. In the regulation equation, we use the instantaneous covariance $`c^{\text{EE}}(t)`$ rather than its time average $`\overline{c}^{\text{EE}}`$ (see Discussion). The time-averaged variance $`\overline{c}^{\text{EE}}`$ is used here for illustration purposes only. In order to obtain an essentially constant value for $`\overline{c}^{\text{EE}}`$ rather than an oscillating function of time, different values of $`\rho `$ are used for the two averaging operations: the kernel used to average $`c^{\text{EE}}`$ into $`\overline{c}^{\text{EE}}`$ is ten times broader than the kernel used to compute $`\overline{s}`$ from $`s`$. along several horizontal lines in the $`(w^{\text{EE}},w^{\text{IE}})`$ plane. As expected, $`\overline{c}^{\text{EE}}`$ is positive only in region $`𝒫`$, where the dynamics is periodic;777In general, positive average covariance across a neuronal population indicates collective fluctuations; in our simplified two-dimensional system, the only possible nontrivial asymptotic behavior is periodic oscillation. although not shown, the same is true of $`\overline{c}^{\text{IE}}`$, the time average of the E-to-I covariance. Note that as $`w^{\text{EE}}`$ crosses the $`𝒪`$-to-$`𝒫`$ boundary (Hopf bifurcation) from left to right, $`\overline{c}^{\text{EE}}`$ increases smoothly from 0 to positive values: as discussed above, the amplitude of the limit cycle at this bifurcation is infinitesimal. In contrast, the change in $`\overline{c}^{\text{EE}}`$ and in $`\overline{c}^{\text{IE}}`$ at $`𝒮`$ (saddlenode bifurcation) is a sharp one, as the system undergoes there a transition from a large limit-cycle regime to a fixed-point attractor.

We now start our study of covariance plasticity by regulating parameter $`w^{\text{EE}}`$ in system 3 while all other parameters, including $`w^{\text{IE}}`$, remain fixed. The system under study then consists of coupled equations 3, 6, 7. Equation 7 prescribes an increase of $`w^{\text{EE}}`$ when $`c^{\text{EE}}>\theta ^{\text{EE}}`$, and a decrease when $`c^{\text{EE}}<\theta ^{\text{EE}}`$. Referring to Figure 3b, we see that to the left of $`𝒮`$, where $`c^{\text{EE}}`$ is high, the first of the two conditions applies; in this region $`w^{\text{EE}}`$ increases. To the right of $`𝒮`$ the covariance vanishes, and $`w^{\text{EE}}`$ decreases. Therefore, $`w^{\text{EE}}(t)`$ is attracted to the transition line $`𝒮`$.888The control parameter $`\theta ^{\text{EE}}`$ should be smaller than the value of $`\overline{c}^{\text{EE}}`$ immediately to the left of $`𝒮`$. The portion of the boundary line where the bifurcation is a subcritical Hopf rather than a saddlenode (footnote 5) yields similar behavior, since the disruption of the large-amplitude limit cycle occurs very near the emergence of point attractors (see also footnote 11).

(Insert Figure 4 around here)

The behavior of this $`w^{\text{EE}}`$ regulation loop is illustrated in Figure 4a for the following setting of parameters: $`w^{\text{EI}}=10`$, $`w^{\text{II}}=6`$, $`\rho =.1`$, $`\theta ^{\text{EE}}=.01`$, $`\epsilon ^{\text{EE}}=.01`$. This figure focuses on a small region of the $`(w^{\text{EE}},w^{\text{IE}})`$ plane, and shows the projection of the trajectory of $`(s,\sigma ,w^{\text{EE}},w^{\text{IE}})`$. Several trajectories are shown, all horizontal since $`w^{\text{IE}}`$ is a constant, These trajectories terminate on the critical line $`𝒮`$, and the behavior of the $`s`$ and $`\sigma `$ components on them is as follows. On the trajectories coming from the left, in the $`𝒫`$ region, $`(s,\sigma )`$ moves along a cyclic orbit, whose amplitude grows as $`w^{\text{EE}}`$ increases and approaches the bifurcation line. On the trajectories coming from the right, in the $`𝒯`$ region, $`(s,\sigma )`$ stays in one of the two point attractors while $`w^{\text{EE}}`$ decreases until it reaches the bifurcation curve. When $`𝒮`$ is reached, either from the left or from the right, motion does not really stop. Rather, $`w^{\text{EE}}`$ sets in a periodic oscillation of small amplitude synchronized with a large-amplitude periodic motion of $`(s,\sigma )`$; the frequency of this oscillation is several orders of magnitude slower than in $`𝒫`$, hence covariance is small—it matches, on average, the control parameter $`\theta ^{\text{EE}}`$. When in this regime, the system spends a long time in one of the two almost-attracting corners of the $`[.5,+.5]^2`$ box before leaving it and moving rapidly to the other corner. This results in an almost-square wave, a behavior that is intermediate between the fast periodic motion observed in $`𝒫`$ and the bistable situation prevailing in $`𝒯`$. The period of this oscillation and the amplitude of the oscillation of $`w^{\text{EE}}`$ depend on parameters $`\rho `$, $`\theta ^{\text{EE}}`$, and $`\epsilon ^{\text{EE}}`$.999Not shown on Figure 4a is the leftmost part of region $`𝒫`$, near the Hopf bifurcation, where the limit cycle is of small amplitude hence the condition $`\overline{c}^{\text{EE}}>\theta ^{\text{EE}}`$ is not realized. When initialized there, the system does not converge to $`𝒮`$. However, in both the $`w^{\text{EE}}`$ and the $`w^{\text{IE}}`$ directions, the domain of attraction of $`𝒮`$ extends to $`+\mathrm{}`$.

We next consider the $`w^{\text{IE}}`$-regulated system, where $`w^{\text{EE}}`$ and all other parameters remain fixed. This system consists of coupled equations 3, 8, 9. As noted, the E-to-I covariance $`c^{\text{IE}}`$ vanishes outside region $`𝒫`$, just like $`c^{\text{EE}}`$; within $`𝒫`$ it varies, in a first approximation, like $`c^{\text{EE}}`$. Since we chose $`\epsilon ^{\text{IE}}`$ to be negative, $`w^{\text{IE}}`$ decreases in $`𝒫`$ and increases in $`𝒯`$, whereas the opposite was true of $`w^{\text{EE}}`$ when it was regulated. Figure 4b shows this $`w^{\text{IE}}`$ dynamics in the same region of the $`(w^{\text{EE}},w^{\text{IE}})`$ plane as before. Parameters are $`w^{\text{EI}}=10`$, $`w^{\text{II}}=6`$, $`\rho =.1`$, $`\theta ^{\text{IE}}=.01`$ and $`\epsilon ^{\text{IE}}=.01`$. The trajectories are now parallel to the $`w^{\text{IE}}`$ axis, and $`(w^{\text{EE}},w^{\text{IE}})`$ is again attracted to the critical line $`𝒮`$ separating region $`𝒫`$ from region $`𝒯`$. This is true only to the left of the vertical asymptote of that curve; trajectories to the right of that line go to $`+\mathrm{}`$.

In sum, regulation of either one of the two parameters $`w^{\text{EE}}`$, $`w^{\text{IE}}`$ has the effect of bringing the system to the critical surface $`𝒮`$ separating the region of oscillation from the region of bistable steady firing; the nullcline diagram is then as in Figure 2b. Note that when the system is on $`𝒮`$, a small perturbation in the weights will elicit either oscillation, constant firing at near-maximum rate, or constant firing at near-minimum rate.

We now turn to the behavior of the system when the two regulation loops act simultaneously; we thus study the system of coupled equations 3, 6, 7, 8, 9. Figure 4c shows the $`(w^{\text{EE}},w^{\text{IE}})`$ dynamics for the same parameters as before, i.e., $`w^{\text{EI}}=10`$, $`w^{\text{II}}=6`$, $`\rho =.1`$, $`\theta ^{\text{EE}}=.01`$, $`\epsilon ^{\text{EE}}=.01`$, $`\theta ^{\text{IE}}=.01`$ and $`\epsilon ^{\text{IE}}=.01`$. It appears from this diagram that the evolution proceeds in two clearly distinct stages. In the first stage, which could be predicted from the study of the regulation loops acting separately, $`(w^{\text{EE}},w^{\text{IE}})`$ moves toward line $`𝒮`$.101010The direction of this linear motion is roughly parallel to the line $`w^{\text{EE}}=w^{\text{IE}}`$. This is because $`\epsilon ^{\text{EE}}=\epsilon ^{\text{IE}}`$, $`\theta ^{\text{EE}}=\theta ^{\text{IE}}`$, and the two covariances $`c^{\text{EE}}`$ and $`c^{\text{IE}}`$ are nearly the same. Another choice of parameters would result in a different slope, but otherwise similar behavior. When this line is reached, motion slows down considerably—typically by several orders of magnitude—and proceeds along the critical line, eventually converging to a point on $`𝒮`$ denoted $`G`$ in Figure 4c. As before, attractor $`G`$ is in reality a slow limit cycle, of small amplitude in $`w^{\text{EE}}`$ and $`w^{\text{IE}}`$, and large amplitude in $`s`$ and $`\sigma `$. All four variables, $`s(t)`$, $`\sigma (t)`$, $`w^{\text{EE}}(t)`$, and $`w^{\text{IE}}(t)`$, are now synchronized; the distinction between slow and fast variables has thus vanished. The basin of attraction of $`G`$ in the $`(w^{\text{EE}},w^{\text{IE}})`$ plane roughly consists of the union of the two domains of attraction of $`𝒮`$ for the separate $`w^{\text{EE}}`$ and $`w^{\text{IE}}`$ regulation dynamics; only the region to the left of and around the Hopf line is not attracted to the saddlenode line $`𝒮`$ and eventually to $`G`$.

(Insert Figure 5 around here)

The location of $`𝒮`$ in the $`(w^{\text{EE}},w^{\text{IE}})`$ plane depends on the values of the fixed parameters $`w^{\text{EI}}`$ and $`w^{\text{II}}`$. The location of the attractor $`G`$ on $`𝒮`$ further depends on the control parameters $`\theta ^{\text{EE}}`$ and $`\theta ^{\text{IE}}`$. When the latter are given identical values, as in the case illustrated in Figure 4c, the attractor $`G`$ has the remarkable property that the $`s`$\- and $`\sigma `$-nullclines stand in near overlap over a large portion of the interval \[-.5,+.5\] (Figure 5); the flow of the system in this configuration nearly vanishes on a large one-dimensional manifold in the two-dimensional phase space. Further, $`s(t)`$ and $`\sigma (t)`$ remain nearly identical at all times.111111Giving different values to parameters $`\theta ^{\text{EE}}`$ and $`\theta ^{\text{IE}}`$ mostly affects the behavior of the system after it has reached $`𝒮`$; if $`\theta ^{\text{EE}}`$ is larger, resp. smaller, than $`\theta ^{\text{IE}}`$, the state moves downward, resp. upward, on $`𝒮`$. When $`(w^{\text{EE}},w^{\text{IE}})`$ is on $`𝒮`$ but above point $`G`$, the nullclines are tangent to each other but do not overlap; such a situation is illustrated in Figure 2b. When $`(w^{\text{EE}},w^{\text{IE}})`$ is on $`𝒮`$ but below point $`G`$, the nullclines do overlap, but over a smaller domain. With $`\theta ^{\text{EE}}=.0118`$ and $`\theta ^{\text{IE}}=.0100`$, the state stabilizes in the narrow three-attractor region described in footnote 5. The state $`(s,\sigma )`$ then visits each of the three ‘attractors’ in turn: its motion consists of a succession of large-amplitude oscillations (periodic attractor) and of spiraling orbits around two symmetric points in the interior of the large cycle (point attractors). The amplitude of the motion of $`(w^{\text{EE}},w^{\text{IE}})`$ remains small. This is a mildly chaotic behavior; a more pronounced chaotic behavior will be described in the next section for the full system.

## 5 Behavior of the regulated full system

Recall that system 3, which we used so far, was derived from system 1 by eliminating the firing thresholds $`h^\text{E}`$ and $`h^\text{I}`$ (equations 2) in such a way as to make $`(.5,.5)`$$`(0,0)`$ in system 3—a center of symmetry of the dynamics. While easier to analyze, the reduced system is less realistic. There is no clear biological justification for this hard-wired symmetry; moreover, when the system is in phase $`𝒯`$, i.e., to the right of the critical surface $`𝒮`$, it can stay for arbitrarily long periods of time in one of the two fixed point attractors, e.g. in the high-activity one; this is unrealistic.

In this section we consider a biologically more plausible way of introducing symmetry in the dynamics. Rather than eliminating the thresholds according to equations 2, we regulate them, thereby implementing a form of ‘soft’ symmetry. Regulating the firing thresholds in a neural network is a simple way to maintain the mean activity around an intermediate, useful, value. This may be viewed as a simplification of the regulation mechanisms at work in real brains, which, in all likelihood, involve systems of inhibitory neurons acting on various time scales.

The simultaneous regulation of four parameters results in a complex dynamics, which makes a thorough analysis impractical. We shall proceed as follows. We first consider, in system 1, the regulation of $`w^{\text{EE}}`$ and $`h^\text{E}`$ for a given setting of all other parameters. We show that the system converges to the intersection of two critical curves, each of which corresponds to the establishement of one point of contact between the nullclines. We next consider the system with all four parameters $`h^\text{E}`$, $`h^\text{I}`$, $`w^{\text{EE}}`$ and $`w^{\text{IE}}`$ regulated, and study the projection of the dynamics on the $`(w^{\text{EE}},w^{\text{IE}})`$ plane. There are again two stages; the first essentially reproduces the behavior observed with the sole $`(w^{\text{EE}},h^\text{E})`$ regulation, while the second is analogous to that observed when regulating $`w^{\text{EE}}`$ and $`w^{\text{IE}}`$ in the reduced system; this applies for a broad range of the remaining fixed parameters $`w^{\text{EI}}`$ and $`w^{\text{II}}`$.

Figure 6a is the bifurcation diagram of system 1 in the $`(w^{\text{EE}},h^\text{E})`$ plane, for the following values of the fixed parameters: $`w^{\text{EI}}=10`$, $`w^{\text{IE}}=10`$, $`w^{\text{II}}=1`$, $`h^\text{I}=5`$. As before, we ignore unstable equilibria and unstable limit cycles. As before there are three regions, denoted respectively by $`𝒪`$, $`𝒯`$ and $`𝒫`$, corresponding to three types of asymptotic behavior: single fixed-point attractor; two fixed-point attractors (high and low activity); one periodic attractor. We now however subdivide region $`𝒪`$—somewhat arbitrarily—according to the location of the fixed-point attractor in the phase space: the three subregions denoted $`𝒪_h`$, $`𝒪_m`$, and $`𝒪_l`$, correspond, respectively, to high, middle, and low activity for this attractor. The transition between region $`𝒫`$ and region $`𝒪_m`$ takes place through the familiar, smooth, Hopf bifurcation. The transition between $`𝒫`$ and $`𝒪_h`$, as well as its continuation between $`𝒪_l`$ and $`𝒯`$, takes place through a saddlenode bifurcation. We denote by $`𝒮_h`$ the locus of this transition; it marks the appearance of a point of contact between the nullclines near the high-activity corner, and is thus similar to the $`𝒮`$ transition in the reduced system. However, due to the symmetry of that system, another point of contact appeared simultaneously near the low-activity corner, giving rise to a double bifurcation. In system 1 this is no more the case, and the intersection of the nullclines near the low-activity corner gives rise to a distinct saddlenode bifurcation line, the transition between $`𝒫`$ and $`𝒪_l`$, which we denote $`𝒮_l`$.

(Insert Figure 6 around here)

When regulating $`w^{\text{EE}}`$ according to equation 7 and leaving all other parameters fixed, the behavior of system 1 is as follows. When starting in region $`𝒫`$ to the left of the critical line $`𝒮_h`$, the system oscillates, covariance is high, hence $`w^{\text{EE}}`$ increases until it reaches the critical line $`𝒮_h`$. A point of contact is then established near the high-activity corner of the square. The system settles in a slow periodic attractor, of small amlitude in $`w^{\text{EE}}`$ and large amplitude in $`(s,\sigma )`$, whereby nearly all the time is spent in the high-activity state.

We now regulate the threshold $`h^\text{E}`$ as well, in such a way as to stabilize $`\overline{s}`$, the time average of $`s`$, around a given target value $`\theta ^\text{E}`$:

$$\frac{dh^\text{E}(t)}{dt}=\epsilon ^\text{E}(\overline{s}(t)\theta ^\text{E}).$$
(10)
The rate constant $`\epsilon ^\text{E}`$ is positive and small, and the control parameter $`\theta ^\text{E}`$ is chosen well in the interior of the interval $`(0,1)`$, e.g. between .2 and .8 (remember that in system 1 the activity variables $`s`$ and $`\sigma `$ lie in the interval $`(0,1)`$). To see how equation 10 achieves the desired regulation, note for instance that, if $`\overline{s}(t)>\theta ^\text{E}`$, $`h^\text{E}`$ will increase, which in turn will result in a decrease of $`\overline{s}(t)`$.

When both $`w^{\text{EE}}`$ and $`h^\text{E}`$ are regulated, the system converges to the intersection of the two critical lines $`𝒮_h`$ and $`𝒮_l`$. In effect, we saw that the full system, when at a generic point of $`𝒮_h`$, stays nearly all the time in the high-activity state; this results in a high value of $`\overline{s}`$. To achieve the condition $`\overline{s}\theta ^\text{E}`$, the equilibrium for equation 10, the system can only be on $`𝒮_l`$ as well.

The joint $`(w^{\text{EE}},h^\text{E})`$ dynamics is illustrated in Figure 6b, for parameters $`w^{\text{EI}}`$, $`w^{\text{IE}}`$, $`w^{\text{II}}`$ and $`h^\text{I}`$ as above, and $`\rho =.2`$, $`\epsilon ^\text{E}=.001`$, $`\theta ^\text{E}=.5`$, $`\epsilon ^{\text{EE}}=.01`$, $`\theta ^{\text{EE}}=.01`$. The intersection of $`𝒮_h`$ and $`𝒮_l`$, denoted $`F`$ in Figure 6b, is reached from all directions in the $`(w^{\text{EE}},h^\text{E})`$ plane. When coming from low $`w^{\text{EE}}`$ values, the system oscillates and converges to $`F`$ through region $`𝒫`$. When coming from high $`w^{\text{EE}}`$ values, the system reaches $`F`$ through region $`𝒯`$, where it bounces back and forth between the high- and low-activity point attractors (an oscillation much slower than in $`𝒫`$).

The nullcline diagram for point $`F`$ of Figure 6b is illustrated in Figure 6c. There are now two points of contact between the nullclines, a situation more degenerate than the one that obtains from regulating $`w^{\text{EE}}`$ only, but ‘equivalent’ to the situation obtained in the reduced system by regulating a single parameter, $`w^{\text{EE}}`$ or $`w^{\text{IE}}`$ (compare Figure 6b to Figure 2b). What characterizes the dynamics at point $`F`$ is that the system is on the verge of oscillation and on the boundary of each of the two steady-firing phases.

We finally consider the system with the four parameters $`h^\text{E}`$, $`h^\text{I}`$, $`w^{\text{EE}}`$ and $`w^{\text{IE}}`$ regulated. We thus include, in addition to equations 1, 6, 7, 8, 9 and 10, a regulation equation for the inhibitory threshold $`h^\text{I}`$:

$$\frac{dh^\text{I}(t)}{dt}=\epsilon ^\text{I}(\overline{\sigma }(t)\theta ^\text{I}).$$
(11)
As in equation 10, the rate constant $`\epsilon ^\text{I}`$ is positive and small, and $`\theta ^\text{I}`$ is chosen in the interval $`(.2,.8)`$, with $`\theta ^\text{I}\theta ^\text{E}`$. The variables now include the activity state $`(s,\sigma )`$ as well as the four regulated parameters $`h^\text{E}`$, $`h^\text{I}`$, $`w^{\text{EE}}`$ and $`w^{\text{IE}}`$.

Figure 7 illustrates the behavior of this system projected on the $`(w^{\text{EE}},w^{\text{IE}})`$ plane, for the following parameter values: $`w^{\text{EI}}=10`$, $`w^{\text{II}}=6`$, $`\rho =.05`$, $`\theta ^{\text{EE}}=.01`$, $`\epsilon ^{\text{EE}}=.01`$, $`\theta ^{\text{IE}}=.01`$, $`\epsilon ^{\text{IE}}=.005`$, $`\theta ^\text{E}=.5`$, $`\epsilon ^\text{E}=.005`$, $`\theta ^\text{I}=.5`$, $`\epsilon ^\text{I}=.002`$. In the sequel, this parameter setting will be referred to as standard. In a first stage, the system converges to a doubly critical point $`F`$ as described above; each such point $`F`$ belongs to the common boundary of the regions of oscillation, high steady firing, and low steady firing. Although we cannot thoroughly characterize the surface of $`F`$ points in the four-dimensional $`(w^{\text{EE}},w^{\text{IE}},h^\text{E},h^\text{I})`$ space as we did in the $`(w^{\text{EE}},w^{\text{IE}})`$ plane for the reduced system, there is, as remarked above, a functional equivalence with the $`𝒮`$ surface. Note that the projection of the $`F`$ surface on the $`(w^{\text{EE}},w^{\text{IE}})`$ plane has a shape quite similar to that of $`𝒮`$ in the reduced system. As before, when the system reaches a point $`F`$, all variables settle in a slow, synchronous, almost-periodic motion. The oscillation of $`s`$ and $`\sigma `$ is a nearly rectangular wave, the system spending nearly all its time in the two corners of the square, where the relative amount of time spent in each corner is determined according to the value of parameter $`\theta ^\text{E}(\theta ^\text{I})`$. As before too, the first stage, which consists of the convergence to a doubly critical point $`F`$, is robust against parameter changes; most parameters can be individually varied over several orders of magnitude without qualitatively affecting this part of the behavior.

The second stage, consisting of a much slower motion on the $`F`$ surface, depends on the values of the various parameters. For most parameter settings, including the standard set (see above), the behavior on this critical surface is a slow, simple, periodic motion, of large amplitude in $`(s,\sigma )`$ and very small amplitude in $`(w^{\text{EE}},w^{\text{IE}})`$. The system eventually settles in a periodic attractor of this simple type, denoted again $`G`$ in Figure 7. Figure 8a shows the $`(s,\sigma )`$ projection of this attractor for the standard parameter set; its $`(w^{\text{EE}},w^{\text{IE}})`$ projection is a small cycle around point $`G`$, whose nullcline diagram is similar to the one shown in Figure 5 (largely overlapping nullclines).

There exists however a small region of parameter space, mostly around $`\epsilon ^\text{I}\epsilon ^\text{E}`$, for which a variety of more complex behaviors are observed during the second stage. The following two cases are examples of such complex behavior. For parameters as above (standard) except that $`\epsilon ^\text{E}=.0051`$, $`\epsilon ^\text{I}=.0046`$, and $`\theta ^{\text{EE}}=.011`$, the system settles in a complex quasi-periodic motion (Figure 8b). For parameters as standard except that $`\epsilon ^\text{I}=\epsilon ^\text{E}=.005`$, the system displays strongly chaotic behavior (Figure 8, c–e).121212This behavior takes place only for some initial values in the $`(w^{\text{EE}},w^{\text{IE}})`$ plane; other initial values converge to a point attractor. Both of these behaviors are actually attractors, reached after considerable time, yet similar behaviors also take place while the system is still moving slowly on the critical surface.

To summarize, both in the reduced and in the full system, convergence to a doubly critical surface between the regions of fast oscillations and of high and low steady firing takes place reliably for a broad range of parameters. Once this doubly critical surface is reached, motion becomes slow, depends on parameters, and, when examined in detail, reveals a variety of behaviors, ranging from simple periodic firing to chaos.

## 6 Discussion

This paper proposes that a regulation mechanism underlies criticality in brain dynamics. In such a scheme, regulation stabilizes the dynamics near an instability. The force driving the system towards criticality is a covariance-governed modification of synaptic efficacies in a recurrent network. Although it has been argued that criticality in some physical systems may be self-organized (Bak et al. 1987), this phenomenon may not be very widespread. The nervous system is moreover regulated homeostatically to withstand perturbations of various sorts. It is then of interest to explain how criticality in brain dynamics may nevertheless arise, from a specific, well-documented, mechanism of synaptic plasticity.

The chief motivation for viewing the dynamics of the nervous sytem as critical is the observation that brains are very sensitive organs. Not only do our brains draw distinctions between stimuli that differ only in minute details, they also allow us to establish subtle yet clear-cut boundaries between cognitive categories at different levels. The various manifestations of hyperacuity in sensory systems may be no more than elaborate forms of signal amplification; however, in higher cognitive functions such as language, abstract thinking, or, for instance, artistic composition, the ability to create a new category by drawing a fine line—an ability which manifests itself early in life and stays with us for a long time—argues in favor of regulated criticality. Such a mechanism appears to be necessary in order to explain how sensitivity is maintained in the face of the profound changes that affect the connectivity of the brain throughout development and learning. If no such mechanism were present, one would expect that the ongoing modification of the networks which carry mental representations would soon bring these networks to generic states; as mentioned in Section 1, a dynamical system in a generic state does not show high susceptibility to external influences.

The emergence of new cognitive categories may in effect be likened to a process of morphogenesis in embryology, or differentiation in cell biology. A biological structure that is about to undergo differentiation is at that particular instant of time unstable, and, as the well-studied mechanism of induction shows, highly susceptible to external signals. From a dynamical-system perspective, the emergence of qualitatively new behavior, e.g., the splitting of one attractor into two, is a bifurcation. The complexification of an individual’s cognitive apparatus in the course of his or her life may be viewed as an open-ended sequence of such bifurcations. Such an interpretation has been defended by René Thom (1975), and related ideas have been expressed by several authors (e.g. van der Maas and Molenaar 1992). Thom (1975) also suggested that structurally stable non-generic singularities may arise from a process he termed the stabilization of thresholds; this process itself would result, in various biological contexts, from the reinforcment of homeostatic mechanisms.131313We thank Jean Petitot for pointing out to us that regulated criticality as proposed here is closely related to Thom’s ideas.

The covariance plasticity rule we use is linear and straightforward. Equation 7 may be viewed as a mean-field version of the covariance rule used in the associative-memory literature (see e.g. Willshaw and Dayan 1990). However, we make a rather different use of this rule. In an associative-memory model, pre- and post-synaptic activities are generally assumed to be independent, yielding a zero expected value for the covariance. Weights are modified according to the instantaneous covariance, and, as noted in Dayan and Sejnowski (1993), storage is marked by the departures of the empirical average of this quantity from its expected value, which is zero. In our model, the expected covariance is positive in the oscillatory phase. The regulation mechanism acts on a slow time scale, and, although we use the instantaneous covariance in the modification rule (Equation 7), we might as well have used the time-averaged covariance; fast variations of the instantaneous covariance are actually smoothed out in the integration of the differential equation. Of course by the very principle of regulation proposed, the system does not dwell in the oscillatory phase; in the regulated state, the average covariance is low.

The other major difference between the situation studied here and the associative-memory paradigm is the assumption of uniform weights. As noted in Section 3, the covariance in our uniform-weight network is simply the variance of the population-averaged activity about its mean, and is always non-negative. This makes it necessary to subtract from it a positive constant $`\theta ^{\text{EE}}`$ in order to allow for decreases of the weights. Thus, whereas in associative-memory models a synaptic weight decreases as a result of negative instantaneous covariance between the pre- and post-synaptic neurons, the condition for weight decrease in our model is that the mean covariance be small or zero, which happens when the system is at rest in a point attractor, of either low or high activity.

The uniform-weight network used in the present study lends itself to a detailed mathematical/numerical analysis. We performed a bifurcation analysis of the continuous-time differential system that describes the behavior of this network in the thermodynamic limit. This analysis (Section 5) reveals, among other features, the existence of a critical surface $`𝒮_h`$ in parameter space, where the system undergoes an abrupt transition from oscillatory behavior to high-rate steady firing. We showed (Section 5) that Hebbian modification of the E-to-E synaptic weights drives the system toward this surface $`𝒮_h`$; this is the main mechanism of regulated criticality proposed.

However, when the system is at a general position on $`𝒮_h`$, it spends most of its time in the high-activity state, which is undesirable. A more realistic situation results from introducing some form of symmetry between the high- and low-activity phases, making use of the firing thresholds. We investigated two ways—formally different but functionally equivalent—to do this. The mathematically simpler way is to enforce an accurate symmetry on the dynamics, by imposing an appropriate relationship between the firing thresholds and the synaptic weights (equations 2). This results in a reduced system, with only four parameters; in this system, there occurs a double bifurcation when the system traverses the critical surface—now denoted $`𝒮`$—separating the oscillatory phase from the bistable, high/low, steady-firing phase (Section 2). Regulation of the sole E-to-E weight brings the system to this doubly critical surface $`𝒮`$ (Section 4).

A biologically more satisfactory solution is to regulate one or both of the firing thresholds so as to control the mean firing rates (Section 5). Thus, when we regulate the threshold for the excitatory neurons in addition to the E-to-E weight, the system converges to the intersection of $`𝒮_h`$ with another critical surface, $`𝒮_l`$, which separates the oscillatory phase from the low-activity fixed-point region. Intersection points between $`𝒮_h`$ and $`𝒮_l`$ are again doubly critical, and they attract the system for a wide range of parameter values.

When the system is on this doubly critical surface, it takes only a small weight perturbation to induce either of the three behaviors: intrinsic oscillation (region $`𝒫`$ of section 5), high activity, quiescence. It is easily seen that, when in this state, the network can also be efficiently driven by a small-amplitude time-varying signal, i.e., an external field; it is thus highly sensitive to input. The stabilization at the boundary of a region of oscillatory behavior appears to be consistent with at least one conclusion that can be safely drawn from the recent literature on cortical oscillations, namely the fact that the precise conditions under which these oscillations occur are difficult to pin down.

We further investigated the effect of regulating the E-to-I weight in addition to the E-to-E weight, according to a similar covariance rule. We showed that regulating these two weights as well as the two firing thresholds results, under appropriate parametric conditions, in convergence to an even more degenerate state. When the system is in that state, its flow vanishes on an entire one-dimensional curve in the two-dimensional phase space, instead of on isolated points. This convergence is slow and parameter-dependent, yet it is interesting to note that when the system is in or near this highly degenerate state it exhibits a range of diverse behaviors, including chaos (Section 5). The chaotic behavior shown in Figure 8, c–e consists of an irregular sequence of spontaneous transitions between the three fundamental phases of the system: oscillatory, high-activity, low-activity.

While the uniform-weight network studied in this paper lends itself to a convenient mathematical analysis, it would be interesting to know whether critical behavior may arise from local covariance plasticity, where synaptic changes are made to depend on pre- and post-synaptic activities relative to individual synapses. This question should be focused by considerations about the elaborate forms of input sensitivity that could play a role in higher brain functions.

Acknowledgments Important contributions in the early stages of this work were made by Howard Gutowitz. We thank Rob de Boer for making available to us the use of the software GRIND, which proved an extremely valuable tool in this study. Various suggestions were made by a number of colleagues, from whom we wish to thank in particular Gérard Weisbuch, Jean Petitot, Gérard Toulouse, Jean-Pierre Nadal, Claude Meunier, David Hansel, Haïm Sompolinsky, and Stuart Geman.

References

Artola, A., Bröcher, S., and Singer, W. 1990. Different voltage-dependent thresholds for inducing long-term depression and long-term potentiation in slices of rat visual cortex. Nature, 347, 69–72.

Bak, P., Tang, C., and Wiesenfeld, K. 1987. Self-Organized Criticality: An Explanation of $`1/f`$ Noise. Phys. Rev. Lett., 59, 381–384.

Bienenstock, E., Cooper, L.N., and Munro, P. 1982. Theory for the development of neuron selectivity: Orientation specificity and binocular interaction in visual cortex. J. Neurosci., 2, 32–48.

Borisyuk, R.M., and Kirillov, A.B. 1992. Bifurcation analysis of a neural network model. Biol. Cybernetics 66, 319–325.

Dayan, P., and Willsahw, D. 1991. Optimising synaptic learning rules in linear associative memories. Biol. Cybernetics 65, 253–265.

Dayan, P., and Sejnowski, T.J. 1993. The Variance of Covariance Rules for Associative Matrix Memories and Reinforcement Learning. Neural Computation, 5, 205–209.

Dudek, S.M., and Bear, M.F. 1992. Homosynaptic long-term depression in area CA1 of hippocampus and effects of N-methyl-D-aspartate receptor blockade. Proc. Natl. Acad. Sci. USA, 89, 4363–4367.

Glauber, R.J. 1963. Time-dependent Statistics of the Ising Model. J. Math. Phys. 4, 294–307.

Hebb D.O. 1949. The Organization of Behavior. Wiley, New York.

Frégnac, Y., Shulz, D., Thorpe, S., and Bienenstock, E. 1988. A cellular analogue of visual cortical plasticity. Nature, 33, 367–370.

Frégnac, Y., Shulz, D., Thorpe, S., and Bienenstock, E. 1992. Cellular analogs of visual cortical epigenesis. I–Plasticity of orientation selectivity. Journal of Neuroscience, 12, 1280–1300.

Guckenheimer, J., and Holmes, P. 1983 Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Springer-Verlag, New York.

Langton, C.R. 1990. Computation at the edge of chaos: Phase transitions and emergent computation. Physica D, 42, 12–37.

Linsker, R. 1986. From basic network principles to neural architecture: Emergence of spatial opponent cells. Proc. Natl. Acad. Sci. USA, 83, 7508–7512.

Metzger, Y., and Lehmann, D. 1990. Learning Temporal Sequences by Local Synaptic Changes. Network: Computation in Neural Systems, 1(3), 169–188.

Metzger, Y., and Lehmann, D. 1994. Learning Temporal Sequences by Excitatory Synaptic Changes Only. Network: Computation in Neural Systems, 5, 89–99.

Rubin, N. 1988. Equilibrium and oscillations in stochastic neural networks. Master’s thesis, The Hebrew University, Jerusalem, Israel (in Hebrew).

Schuster, H.G., and Wagner, P. 1990 A model for neuronal oscillations in the visual cortex. 1. Mean-field theory and derivation of the phase equations. Biological Cybernetics, 64, 77–82.

Sejnowski, T.J. 1977a. Storing covariance with nonlinearly interacting neurons. J. Math. Biol., 4, 303–321.

Sejnowski, T.J. 1977b. Statistical constraints on synaptic plasticity. J. Theor. Biol., 69, 385–389.

Stanton, P.K. and Sejnowski, T.J. 1989. Associative long-term depression in the hippocampus induced by hebbian covariance. Nature, 339, 215, May 1989.

Thom, R. 1975. Structural Stability and Morphogenesis: an Outline of a General Theory of Models. W. A. Benjamin, Reading, Mass.

van der Maas, H.L.J., and Molenaar, P.C.M. 1992. Stagewise Cognitive Development: An Application of Catastrophe Theory. Psychological Review, 99(3), 395–417.

Willshaw, D., and Dayan, P. 1990. Optimal Plasticity from Matrix Memories: What Goes Up Must Come Down. Neural Computation, 2, 85–93.

> Figure 1: Mean activities of excitatory and inhibitory populations in a moderate-size uniform-weight system exhibiting oscillatory behavior ($`N=70`$; Glauber dynamics).

> Figure 2: Behavior of mean-field system for different values of the E-to-E synaptic weight $`w^{\text{EE}}`$. Diagrams show trajectories and nullclines. (a) $`w^{\text{EE}}=12`$ (parameters are as in Figure 1); all trajectories converge to a limit cycle. (b) $`w^{\text{EE}}=\widehat{w}_{\text{sn}}^{\text{EE}}14.22`$; the system is at the saddlenode bifurcation: nullclines are tangent to each other (no trajectories shown). (c) $`w^{\text{EE}}=15`$; nullclines intersect, the periodic attractor has vanished, two point attractors have appeared.

> Figure 3: Bifurcation diagram in $`(w^{\text{EE}},w^{\text{IE}})`$ plane. (a) Diagram shows three regions, characterizing different attractor configurations. Region $`𝒪`$: single point attractor, of intermediate activity level; region $`𝒫`$: periodic attractor, as depicted in Figure 2a; region $`𝒯`$: two point attractors, of high and low activity, as depicted in Figure 2c. Transitions between regions occur through bifurcations, of Hopf type between $`𝒪`$ and $`𝒫`$, of saddlenode type between $`𝒫`$ and $`𝒯`$ (curve $`𝒮`$), and of pitchfork type between $`𝒪`$ and $`𝒯`$. (b) Average covariance along four different lines of constant $`w^{\text{IE}}`$ in the $`(w^{\text{EE}},w^{\text{IE}})`$ plane. Note the sharp variation of the covariance on the critical line $`𝒮`$ separating $`𝒫`$ from $`𝒯`$.

> Figure 4: Regulation of system 3 by covariance plasticity. (a) $`w^{\text{EE}}`$ is regulated, $`w^{\text{IE}}`$ is constant: state converges to critical surface $`𝒮`$. (b) $`w^{\text{IE}}`$ is regulated, $`w^{\text{EE}}`$ is constant: state converges to critical surface $`𝒮`$. (c) both $`w^{\text{EE}}`$ and $`w^{\text{IE}}`$ are regulated: state converges to a particular point, denoted $`G`$, on critical surface $`𝒮`$.

> Figure 5: Nullcline diagram at point $`G`$ (see figure 4c). Nullclines overlap almost perfectly over much of the interval $`[.5,.5]`$.

> Figure 6: Regulation of two parameters in system 1. (a) Bifurcation diagram in $`(w^{\text{EE}},h^\text{E})`$ plane. (b) Regulation of $`w^{\text{EE}}`$ and $`h^\text{E}`$ causes convergence to point $`F`$, the intersection of critical lines $`𝒮_h`$ and $`𝒮_l`$. (c) Nullcline diagram at $`F`$.

> Figure 7: Behavior of full system under simultaneous regulation of four parameters. Diagram shows projection on $`(w^{\text{EE}},w^{\text{IE}})`$ plane, illustrating the similarity of behavior with reduced system (compare with Figure 4c, but note difference of scales). Limits of the attraction basin to the left are roughly indicated by the starting points of the trajectories shown; attraction basin is unbounded in all other directions.

|  |  |
| --- | --- |
|  |  |
|  |  |

> Figure 8: Various behaviors of regulated full system after it has reached critical surface (Figure 7). Diagrams show $`s\left(t\right)`$ for three slightly different parameter settings (see text); in all cases, the projection of the motion on the $`(w^{\text{EE}},w^{\text{IE}})`$ plane remains of small amplitude. (a) Simple periodic attractor, point $`G`$ of Figure 7; similar periodic attractors are reached for most parameter settings. (b) Complex quasi-periodic attractor. (c–e) Chaotic attractor; for a given parameter setting, three diagrams corresponding to different instants of time and different lengths of time; characteristic are the irregular transitions between the high-activity, low-activity, and oscillatory phases.


format: text

text:
# Long-range magnetic interaction due to the Casimir effect

## Abstract

The zero-point quantum fluctuations of the electromagnetic field in vacuum are known to give rise to a long-range attractive force between metal plates (Casimir effect). For ferromagnetic layers separated by vacuum, it is shown that the interplay of the Casimir effect and of the magneto-optical Kerr effect gives rise to a long-range magnetic interaction. The Casimir magnetic force is found to decay as $`D^1`$ in the limit of short distances, and as $`D^5`$ in the limit of long distances. Explicit expressions for realistic systems are given in the large and small distance limits. An experimental test of the Casimir magnetic interaction is proposed.
published in: Phys. Rev. Lett. 88, 240401 (2002)

Since the discovery of magnets by the ancient Greeks, long-range magnetic interactions have been an object of fascination. It is usually considered that there exists essentially two kinds of magnetic interactions between magnetic moments or magnetized bodies: (i) the dipole-dipole magnetostatic interaction, and (ii) the electron-mediated exchange interaction. The latter have recently received a renewed attention, with the discovery of a spectacular oscillatory behavior of the interlayer exchange coupling between ferromagnetic layers separated by a non-magnetic metal spacer , due to a spin-dependent quantum size effect .

For the case of two uniformly magnetized ferromagnetic plates (of infinite lateral extension) held parallel to each other in vacuum, the two above mentioned magnetic coupling mechanisms yield a magnetic interaction which decreases exponentially with interplate distance $`D`$: (i) the stray field due to a uniformly magnetized plate decreases exponentially (with a characteristic decay length of the oder of the interatomic distance) with the distance from the plate, and so does also the interplate dipolar interaction; (ii) the interplate exchange interaction also decays exponentially with $`D`$, since it is mediated by electrons tunneling through vacuum between the two plates. The aim of the present Letter is to point out the existence of a novel, so far overlooked, mechanism of magnetic interaction between magnetized (i.e., ferro- or ferrimagnetic) bodies. This interaction arises from the Casimir effect, and gives rise to a long-range (i.e., with power-law decay) magnetic interaction; at sufficiently large distance $`D`$ it is therefore the dominant mechanism of magnetic interaction between the two ferromagnetic plates.

As Casimir pointed out in a seminal paper , the zero-point quantum fluctuations of the electromagnetic (EM) field in vacuum leads to observable effects: changing the boundary condition of the EM field (e.g., by moving a body with respect to another) yields a finite change of the (infinite) zero-point energy of the system, and therefore results in an observable force. In particular, Casimir predicted the existence of a long-range attractive force between mirrors in vacuum. The Casimir effect is currently attracting considerable interest , in particular with respect to micromechanical devices, and has deep implications in many fields of physics.

The boundary condition of the EM field can, however, also be modified without any mechanical displacement, but rather, by changing the order parameter of a collective ordering phenomenon such as ferromagnetism. When the two mirrors are ferromagnetic, the magneto-optical Kerr effect influences the boundary condition of the EM field, so that the Casimir effect manifests as a energy difference (per unit area) $`\mathrm{\Delta }_{AF}_{FM}`$ between the configurations in which the two mirrors have their magnetizations antiparallel (AF) or parallel (FM) to each other, i.e., as a magnetic interaction. Alternatively, the above effect manifests as a dependence of the (mechanical) Casimir force (per unit area) among the mirrors upon the relative orientation of their magnetizations: $`\mathrm{\Delta }_{AF}_{FM}=\mathrm{d}\mathrm{\Delta }/\mathrm{d}D`$. The Casimir magnetic force $`\mathrm{\Delta }`$ has been calculated for ferromagnetic mirrors described by a Drude model: in the limit of very large distance, the magnetic force decays as $`D^5`$; in the limit of small distances, it decays as $`D^1`$; in the intermediate regime, it decays as $`D^4`$. For equivalent ferromagnetic materials on both sides, the Casimir magnetic interaction is always antiferromagnetic. For realistic systems, the explicit expression of the Casimir magnetic force is found to be

$$\mathrm{\Delta }\frac{3\zeta (3)}{16\pi ^3}\frac{\mathrm{}c^2}{D^5}\frac{\theta _A\theta _B}{\sqrt{\sigma _A\sigma _B}},$$
(1)
in the limit of large distance, where $`\sigma _{A(B)}`$ and $`\theta _{A(B)}`$ are, respectively, the dc conductivity and anomalous Hall angle of ferromagnetic plate $`A`$ (resp. $`B`$), and

$$\mathrm{\Delta }\frac{1}{4\pi ^2}\frac{\mathrm{}}{c^2D}_0^+\mathrm{}\frac{\omega ^2\epsilon _{xy}^A(\mathrm{i}\omega )\epsilon _{xy}^B(\mathrm{i}\omega )}{\left[1+\epsilon _{xx}^A(\mathrm{i}\omega )\right]\left[1+\epsilon _{xx}^B(\mathrm{i}\omega )\right]}d\omega ,$$
(2)
in the limit of short distances, where $`\epsilon _{xx}^{A(B)}(\mathrm{i}\omega )`$ and $`\epsilon _{xy}^{A(B)}(\mathrm{i}\omega )`$ are, respectively, the diagonal and off-diagonal elements of the dielectric tensor of ferromagnetic plate $`A`$ (resp. $`B`$), evaluated at imaginary frequency $`\mathrm{i}\omega `$.

The Casimir interaction energy (per unit area) between two mirrors can be conveniently expressed in terms of their reflection coefficients as ,

$$=_0^+\mathrm{}d\omega \frac{\mathrm{}}{2}\mathrm{coth}\left(\frac{\mathrm{}\omega }{k_BT}\right)\frac{1}{(2\pi )^2}\mathrm{d}^2𝐤_{}\frac{1}{\pi }\text{Im Tr}\mathrm{ln}\left(1𝖱_A𝖱_B\mathrm{e}^{2ik_{}D}\right),$$
(3)
where $`𝐤_{}`$ and $`k_{}`$ are the components of the wavevector, respectively, parallel and perpendicular to the mirrors. The above expression is completely analogous to the one derived independently for the case of interactions mediated by fermions (electrons) . The $`2\times 2`$ matrix of reflection coefficients on mirror $`A`$ (resp. $`B`$) is given by

$$𝖱_{A(B)}\left(\begin{array}{cc}r_{ss}^{A(B)}r_{sp}^{A(B)}& \\ r_{ps}^{A(B)}r_{pp}^{A(B)}& \end{array}\right),$$
(4)
where the index $`s`$ (resp. $`p`$) corresponds to a polarization with the electric field perpendicular (resp. parallel) to the incidence plane. The off-diagonal matrix elements $`r_{sp}`$ and $`r_{ps}`$ are responsible for the magneto-optical effects. With the usual convention that the $`s`$ axis remains unchanged upon reflection, one has $`r_{sp}=r_{ps}`$, and, for perpendicular incidence, $`r_{ss}=r_{pp}`$. Performing the change of variables $`(\omega ,𝐤_{})(\omega ,k_{})`$ and using complex plane integration methods, one can rewrite Eq. (3) at $`T=0`$ as

$$=\frac{\mathrm{}}{4\pi ^2}_0^+\mathrm{}k_{}dk_{}_0^{k_{}c}d\omega \text{Re Tr }\mathrm{ln}\left[1𝖱_A(\mathrm{i}\omega ,\mathrm{i}k_{})𝖱_B(\mathrm{i}\omega ,\mathrm{i}k_{})\mathrm{e}^{2k_{}D}\right],$$
(5)
where the reflection coefficients are evaluated at imaginary values of the frequency and normal wavevector.

For a mirror magnetized along its normal (pointing outwards), the reflection coefficients are given by

$`r_{ss}(\mathrm{i}\omega ,\mathrm{i}k_{})`$ $`=`$ $`{\displaystyle \frac{k_{}c\sqrt{\omega ^2\left(\epsilon _{xx}(\mathrm{i}\omega )1\right)+\left(k_{}c\right)^2}}{k_{}c+\sqrt{\omega ^2\left(\epsilon _{xx}(\mathrm{i}\omega )1\right)+\left(k_{}c\right)^2}}},r_{pp}(\mathrm{i}\omega ,\mathrm{i}k_{})={\displaystyle \frac{\epsilon _{xx}(\mathrm{i}\omega )k_{}c\sqrt{\omega ^2\left(\epsilon _{xx}(\mathrm{i}\omega )1\right)+\left(k_{}c\right)^2}}{\epsilon _{xx}(\mathrm{i}\omega )k_{}c+\sqrt{\omega ^2\left(\epsilon _{xx}(\mathrm{i}\omega )1\right)+\left(k_{}c\right)^2}}}`$ (7)
$`r_{sp}(\mathrm{i}\omega ,\mathrm{i}k_{})`$ $`=`$ $`r_{ps}(\mathrm{i}\omega ,\mathrm{i}k_{})={\displaystyle \frac{k_{}c\omega \epsilon _{xy}(\mathrm{i}\omega )}{\left[k_{}c+\sqrt{\omega ^2\left(\epsilon _{xx}(\mathrm{i}\omega )1\right)+\left(k_{}c\right)^2}\right]\left[\epsilon _{xx}(\mathrm{i}\omega )k_{}c+\sqrt{\omega ^2\left(\epsilon _{xx}(\mathrm{i}\omega )1\right)+\left(k_{}c\right)^2}\right]}}.`$ (8)
For the sake of simplicity, the arguments $`(\mathrm{i}\omega ,\mathrm{i}k_{})`$ will be omitted below. If the magnetization point inwards, then the sign of $`r_{sp}`$ and $`r_{ps}`$ is reversed.

As the magneto-optical reflection coefficients $`r_{sp}`$ are usually much smaller than 1 and than the usual reflection coefficients $`r_{ss}`$ and $`r_{pp}`$, one can expand the Casimir magnetic energy to lowest order in the magneto-optical coefficients, yielding

$$\mathrm{\Delta }=\frac{\mathrm{}}{\pi ^2}_0^+\mathrm{}k_{}dk_{}_0^{k_{}c}d\omega \text{ Re}\left[\frac{r_{sp}^Ar_{sp}^B\mathrm{e}^{2k_{}D}}{\left(1r_{ss}^Ar_{ss}^B\mathrm{e}^{2k_{}D}\right)\left(1r_{pp}^Ar_{pp}^B\mathrm{e}^{2k_{}D}\right)}\right].$$
(9)
Eq. (9) together with Eqs. (7,b) allow to calculated the Casimir magnetic energy and force.

In order to illustrate the above result, let us calculate the Casimir magnetic interaction for the case of (equivalent) ferromagnetic mirrors with a dielectric tensor approximated by a Drude model:

$$\epsilon _{xx}(\mathrm{i}\omega )=1+\frac{\omega _{p}^{}{}_{}{}^{2}\tau }{\omega (1+\omega \tau )},\epsilon _{xy}(\mathrm{i}\omega )=\frac{\omega _{p}^{}{}_{}{}^{2}\omega _c\tau ^2}{\omega (1+\omega \tau )^2}.$$
(10)
The plasma frequency $`\omega _p`$ is given by $`\omega _{p}^{}{}_{}{}^{2}4\pi ne^2/m^{}`$; $`\omega _ceB_{\mathrm{eff}}/m^{}c`$ is the cyclotron frequency, where $`B_{\mathrm{eff}}`$ is the effective magnetic field experienced by conduction electrons as a result of the combined effect of the exchange and spin-orbit interactions; $`\tau `$ is the relaxation time. It is assumed that $`\omega _c\tau 1\omega _p\tau `$, which constitutes the usual situation.

One can distinguish three different regimes: (i) $`Dc\tau `$, (ii) $`c/\omega _pDc\tau `$, (iii) $`Dc/\omega _p`$. In regime (i) (i.e., at long distances), the integral in Eq. (9) is dominated by the range $`\omega k_{}cc/D1/\tau `$, for which one has

$$\epsilon _{xx}(\mathrm{i}\omega )\frac{\omega _{p}^{}{}_{}{}^{2}\tau }{\omega }1,\epsilon _{xy}(\mathrm{i}\omega )\frac{\omega _{p}^{}{}_{}{}^{2}\omega _c\tau ^2}{\omega },$$
(11)
so that

$$r_{ss}r_{pp}1,r_{sp}\frac{\omega _c}{\omega _p}\sqrt{\omega \tau }.$$
(12)
One then finds that

$$\mathrm{\Delta }\frac{3\zeta (3)}{16\pi ^2}\frac{\mathrm{}c^2}{D^4}\frac{\omega _{c}^{}{}_{}{}^{2}\tau }{\omega _{p}^{}{}_{}{}^{2}},\mathrm{\Delta }\frac{3\zeta (3)}{4\pi ^2}\frac{\mathrm{}c^2}{D^5}\frac{\omega _{c}^{}{}_{}{}^{2}\tau }{\omega _{p}^{}{}_{}{}^{2}},$$
(13)
where $`\zeta (x)_{n=1}^{\mathrm{}}n^x`$ is the Riemann zeta function, with $`\zeta (3)1.202\mathrm{}`$

In the intermediate distance regime (ii) ($`c/\omega _pDc\tau `$), the integral in Eq. (9) is dominated by the range $`1/\tau \omega k_{}cc/D\omega _p`$, for which one has

$$\epsilon _{xx}(\mathrm{i}\omega )\frac{\omega _{p}^{}{}_{}{}^{2}}{\omega ^2}1,\epsilon _{xy}(\mathrm{i}\omega )\frac{\omega _{p}^{}{}_{}{}^{2}\omega _c}{\omega ^3},$$
(14)
so that

$$r_{ss}r_{pp}1,r_{sp}\frac{\omega _c}{\omega _p}.$$
(15)
One then finds that

$$\mathrm{\Delta }\frac{1}{24}\frac{\mathrm{}c}{D^3}\frac{\omega _{c}^{}{}_{}{}^{2}}{\omega _{p}^{}{}_{}{}^{2}},\mathrm{\Delta }\frac{1}{8}\frac{\mathrm{}c}{D^4}\frac{\omega _{c}^{}{}_{}{}^{2}}{\omega _{p}^{}{}_{}{}^{2}}.$$
(16)
In the short distance regime (iii) ($`Dc/\omega _p`$), one needs to consider separately the range with $`\omega k_{}c\omega _p`$, for which the reflection coefficients are given by Eqs. (15), and the range with $`\omega _p\omega k_{}c`$. For $`\omega \omega _p`$, $`\epsilon _{xy}`$ is given by Eq. (14) and

$$\epsilon _{xx}(\mathrm{i}\omega )1\frac{\omega _{p}^{}{}_{}{}^{2}}{\omega ^2}1,$$
(17)
so that, for $`k_{}c\omega _p`$,

$$|r_{ss}|1,r_{pp}1,r_{sp}\frac{\omega _{p}^{}{}_{}{}^{2}\omega _c}{2k_{}c(2\omega ^2+\omega _{p}^{}{}_{}{}^{2})}.$$
(18)
One then finds that

$`\mathrm{\Delta }`$ $``$ $`{\displaystyle \frac{1}{16\pi \sqrt{2}}}{\displaystyle \frac{\mathrm{}}{c^2}}\mathrm{ln}\left({\displaystyle \frac{c}{\omega ^{}D}}\right)\omega _{c}^{}{}_{}{}^{2}\omega _p,`$ (20)
$`\mathrm{\Delta }`$ $``$ $`{\displaystyle \frac{1}{16\pi \sqrt{2}}}{\displaystyle \frac{\mathrm{}}{c^2D}}\omega _{c}^{}{}_{}{}^{2}\omega _p,`$ (21)
where $`\omega ^{}`$ is a cut-off frequency of the order of the plasma frequency $`\omega _p`$.

For realistic systems, it is in general necessary to perform a detailed calculation. However, in the limit of large and small distances, explicit expressions can be obtained. At large distances (i.e., for $`Dc/\tau `$), the Casimir magnetic force is essentially determined by the dielectric tensor $`\epsilon _{ij}(\mathrm{i}\omega )=\delta _{ij}+4\pi \sigma _{ij}(\mathrm{i}\omega )/\omega `$ at low imaginary frequency. In this regime, one can safely approximate the conductivity tensor $`\sigma _{ij}(\mathrm{i}\omega )`$, in the above expression, by its dc value: $`\sigma _{xx}(\mathrm{i}\omega )\sigma `$, $`\sigma _{xy}(\mathrm{i}\omega )\sigma \theta `$, where $`\sigma `$ is the dc conductivity, and $`\theta `$ the anomalous Hall angle of the ferromagnetic mirror. Proceeding as for the Drude model, one then obtains,

$$\mathrm{\Delta }\frac{3\zeta (3)}{64\pi ^3}\frac{\mathrm{}c^2}{D^4}\frac{\theta _A\theta _B}{\sqrt{\sigma _A\sigma _B}},$$
(22)
from which Eq. (1) follows immediately.

For short distances, the Casimir magnetic interaction is dominated by imaginary wavevectors $`\mathrm{i}k_{}`$ with $`\omega ^{}k_{}cc/D`$, where the cut-off frequency $`\omega ^{}`$ is of the order of the plasma frequency, or the typical frequency of interband transitions. In this regime, one has

$$|r_{ss}|1,r_{pp}1,r_{sp}\frac{\omega \epsilon _{xy}(\mathrm{i}\omega )}{2k_{}c\left[1+\epsilon _{xx}(\mathrm{i}\omega )\right]}.$$
(23)
One eventually obtains,

$$\mathrm{\Delta }\frac{1}{4\pi ^2}\frac{\mathrm{}}{c^2}\mathrm{ln}(\frac{c}{\omega ^{}D})_0^+\mathrm{}\frac{\omega ^2\epsilon _{xy}^A(\mathrm{i}\omega )\epsilon _{xy}^B(\mathrm{i}\omega )}{\left[1+\epsilon _{xx}^A(\mathrm{i}\omega )\right]\left[1+\epsilon _{xx}^B(\mathrm{i}\omega )\right]}d\omega ,$$
(24)
from which Eq. (2) follows immediately.

Let us now discuss whether the novel Casimir magnetic interaction can be observed experimentally. Obviously, the regime of potential experimental interest is the short distance limit. To obtain a rough estimate of the magnitude of the effect, it is sufficient to approximate the (magneto-)optical absorption spectrum by a single absorption line at frequency $`\omega _0`$ containing all the spectral weight, i.e., we write:

$`\text{Im }\epsilon _{xx}(\omega )`$ $``$ $`\omega _0\epsilon _{xx}^{\mathrm{eff}}\delta (\omega \omega _0),`$ (26)
$`\text{Re }\epsilon _{xy}(\omega )`$ $``$ $`\omega _0\epsilon _{xy}^{\mathrm{eff}}\delta (\omega \omega _0).`$ (27)
This is expected to be a good approximation in the limit of small distances (i.e., high frequencies) where the details of the (magneto-)optical spectra should not matter too much. The dielectric tensor at imaginary frequency is then obtained from the causality relations

$`\epsilon _{xx}(\mathrm{i}\omega )`$ $`=`$ $`1+{\displaystyle \frac{2}{\pi }}{\displaystyle _0^+\mathrm{}}d\omega ^{}{\displaystyle \frac{\omega ^{}\text{ Im }\epsilon _{xx}(\omega ^{})}{\omega _{}^{}{}_{}{}^{2}+\omega ^2}},`$ (29)
$`\epsilon _{xy}(\mathrm{i}\omega )`$ $`=`$ $`{\displaystyle \frac{2}{\omega \pi }}{\displaystyle _0^+\mathrm{}}d\omega ^{}{\displaystyle \frac{\omega _{}^{}{}_{}{}^{2}\text{ Re }\epsilon _{xy}(\omega ^{})}{\omega _{}^{}{}_{}{}^{2}+\omega ^2}},`$ (30)
and one eventually obtains

$`\mathrm{\Delta }`$ $``$ $`{\displaystyle \frac{1}{16\pi ^3}}{\displaystyle \frac{\mathrm{}}{c^2}}\mathrm{ln}\left({\displaystyle \frac{c}{\omega _0D}}\right){\displaystyle \frac{\omega _{0}^{}{}_{}{}^{3}\epsilon _{xy}^{\mathrm{eff}}{}_{}{}^{2}}{\left(1+\epsilon _{xx}^{\mathrm{eff}}/\pi \right)^{3/2}}},`$ (32)
$`\mathrm{\Delta }`$ $``$ $`{\displaystyle \frac{1}{16\pi ^3}}{\displaystyle \frac{\mathrm{}}{c^2D}}{\displaystyle \frac{\omega _{0}^{}{}_{}{}^{3}\epsilon _{xy}^{\mathrm{eff}}{}_{}{}^{2}}{\left(1+\epsilon _{xx}^{\mathrm{eff}}/\pi \right)^{3/2}}}.`$ (33)
By simple inspection of the (magneto-)optical absorption spectra of transition metal ferromagnets , one finds that the model parameters assume the typical values $`\omega _06\times 10^{15}\mathrm{s}^1`$, $`\epsilon _{xx}^{\mathrm{eff}}10`$, $`\epsilon _{xy}^{\mathrm{eff}}1.5\times 10^2`$. Experimentally, it is usually not convenient to maintain two plates accurately parallel to each other, so that a configuration with a planar mirror and a lens-shaped mirror is usually adopted. Even in this configuration, the parasitic magnetostatic interaction can be made as small as needed by taking a uniformly magnetized plate of sufficiently low thickness and sufficiently large lateral extension. The net resulting Casimir magnetic force $`\mathrm{\Delta }F`$ (not to be confused with the Casimir magnetic force per unit area $`\mathrm{\Delta }`$) is then obtained by means of the “proximity force theorem” : $`\mathrm{\Delta }F=2\pi R\mathrm{\Delta }(D),`$ where $`R`$ is the curvature radius of the lens-shaped mirror and $`D`$ the shortest distance. For $`R=100\mu `$m and $`Dc/\omega _050`$ nm, one finds that $`|\mathrm{\Delta }F|10`$ fN, with only a weak (logarithmic) dependence upon $`D`$.

The Casimir magnetic force $`\mathrm{\Delta }F`$ is several orders of magnitude weaker than the non-magnetic Casimir force $`F`$ ($`|F|0.5`$ nN, for $`D=50`$ nm); however, it can be detected independently of $`F`$ by using a resonant differential method, as done in “magnetic resonant force microscopy” : in this approach, one measures the mechanical force between two magnetic samples, one of them being fixed and magnetically hard, the other one being magnetically soft and attached to a high-$`Q`$ mechanical resonator (of resonance frequency $`\omega _r`$) consisting of a micro-cantilever. By modulating the magnetization of the soft sample by means of an ac magnetic field (the hard sample remaining unchanged) at $`\omega =\omega _r`$, one can detect the magnetic force among the two samples (i.e., the difference $`\mathrm{\Delta }FF_{AF}F_{FM}`$) with a considerably higher sensitivity than allowed by a dc measurement of $`F_{AF}`$ or $`F_{FM}`$ separately. In addition, when measuring the non-magnetic Casimir effect, great care has to be taken to eliminate parasitic electrostatic interactions, which is done automatically in the approach proposed here. A detailed discussion of the sensitivity limitations of ”magnetic resonant force microscopy” is given in Ref. : it is limited on one hand by the sensitivity in measuring the deflection of cantilever, which yields $`|\mathrm{\Delta }F_{\text{min}}|k\delta x/Q`$ (where $`k`$ is the cantilever spring constant, $`Q`$ the quality factor, and $`\delta x`$ the deflection sensitivity), and by the thermomechanical noise on the other hand, which yields $`|\mathrm{\Delta }F_{\text{min}}|\sqrt{4k_BT\mathrm{\Delta }\nu k/(\omega _rQ)}`$, for a bandwidth $`\mathrm{\Delta }\nu `$. For the cantilever used in Ref.  ($`k=1`$ mN/m, $`Q=3000`$, $`\omega _r=1.4`$ kHz), a force sensitivity of 0.5 fN at room temperature was obtained, to be compared with the sensitivity of at best 1 pN reported in Ref.  for the dc measurement of the non-magnetic Casimir effect. Prospects for further improvement in the force sensitivity in ”magnetic resonant force microscopy” up to $`3\times 10^{17}`$ N$`/\sqrt{\text{Hz}}`$ at room temperature and $`4\times 10^{18}`$ N$`/\sqrt{\text{Hz}}`$ at 4.2 K are discussed in Ref. . Indeed, force sensitivity in the attonewton ($`10^{18}`$ N) range has recently been demonstrated .

The cantilever used in Ref.  consisted of a 50 $`\mu `$m-long, 5 $`\mu `$m-wide and 90 nm-thick Si beam terminated by a square paddle of 30 $`\mu `$m side-length. Such a cantilever would be appropriate for the proposed experiment. By depositing a droplet of polymer on the paddle, it would be possible to produce a lens-shaped substrate of suitable curvature radius ($``$100 $`\mu `$m), which could then be covered by evaporation with a thin ($``$10 nm) layer of a soft ferromagnet such as permalloy. Great care has to be taken to minimize the parasitic magnetostatic interaction: the hard magnetic plate should be as wide, thin, and magnetically uniform as possible. One should therefore chose a material with a high coercivity and 100% remanence. A thin ($``$ 10 nm) layer of CoPt alloy or a Co/Pt multilayer would be suitable. For a plate radius of 1 cm, the parasitic magnetostatic force can be estimated to be below 1 attonewton, which is sufficient for the present purpose. Finally, the parasitic force due to the interaction of the soft ferromagnet with the ac field would yield a signal at 2$`\omega `$ and would therefore be filtered out by lock-in detection.

To conclude, the above discussion suggests that the experimental test of the novel Casimir magnetic interaction would indeed be possible by using currently available techniques.