Paradoxical self-sustained dynamics emerge from orchestrated excitatory and inhibitory homeostatic plasticity rules

Significance Cortical networks have the remarkable ability to self-assemble into dynamic regimes in which excitatory positive feedback is balanced by recurrent inhibition. This inhibition-stabilized regime is increasingly viewed as the default dynamic regime of the cortex, but how it emerges in an unsupervised manner remains unknown. We prove that classic forms of homeostatic plasticity are unable to drive recurrent networks to an inhibition-stabilized regime due to the well-known paradoxical effect. We next derive a novel family of cross-homeostatic rules that lead to the unsupervised emergence of inhibition-stabilized networks. These rules shed new light on how the brain may reach its default dynamic state and provide a valuable tool to self-assemble artificial neural networks into ideal computational regimes.

(C) Example simulation of the HAAA rule over the course of simulated development. The evolution of the firing rate of the excitatory and inhibitory population within a trial in response to a brief external input is shown in every plot. Eset=5 and Iset=14 represent the target homeostatic setpoints. Weights were initialized to WEE=2.1, WEI=3, WIE=4, and WII=2 as in Fig. 2. Note that while an evoked self-sustained activity emerges by Trial 200 the firing rates do not converge to their setpoints, and by Trial 500 the stable activity is no longer observed. (D) Average rate across trials (upper plot) for the excitatory and inhibitory populations for the data shown in C. Weight dynamics (bottom plot) produced by the homeostatic rules across trials for the data shown in C. (E) Average final rate for 100 independent HAAA simulations with different weight initializations.
Those initializations included cases in which the network starts in the sublinear regime (where the initial E firing rate was zero or very low). The weights were initialized uniformly between the following ranges: [1,3], [0.5,1.5], [4,8], [0.2,0.8]. Data represents mean ± SEM. (F) Analytical stability regions of the neural and HAAA plasticity rule subsystems as a function of the free weights WEE and WIE. Here the stability plot is obtained by considering equal learning rates for all four plasticity rules (as used for panels C-E).
(G) Similar to F but with but with E << I. Right of blue line shows the area where the network is in a paradoxical regime (defined by the condition WEE * gE -1 > 0). Contrary to standard homeostatic rules (Fig. 2), the HAAA rule is only stable in the paradoxical region of parameter space (i.e., WEE*gE -1 > 0; note white area to the left of the blue line). This may explain why the rule fails at driving the network to a stable self-sustained activity when starting with developmental-like conditions.    (A) Same example as in Fig. 6B, where a Hebbian change in the connectivity matrix is introduced at Trial 200. A 'memory' is imprinted in the first 10 excitatory neurons by increasing their recurrent weights by a constant factor. This change drives the network away from its setpoints and reengages the two-term cross-homeostatic rules. The rules successfully bring the network rates back to the setpoints, while preserving much of the differential connectivity between the altered weights. Setpoints are shown in dashed lines, Eset=5 and Iset=14. (B) Weight matrix of the network at four different time points labeled in A. a) Initial random weight matrix. b) The weight matrix after two-term cross homeostatic plasticity has driven the network rates to setpoints. c) A Hebbian-like change in connectivity is imposed into the recurrent weights of the first 10 excitatory neurons. d) Weight matrix after two-term cross-homeostatic plasticity re-stabilizes the firing rates. Note that the 'memory' imposed into the weight matrix is preserved. Figure S7. Broader weight initializations for the multi-unit model also lead to convergence of the cross-homeostatic rules.

Numerical simulations of the firing rate model
For all simulations, the weights were updated after the completion of every trial. The trials lasted 2 seconds. Note that the value of and on every rule are implemented as average firing rates. The average of and is computed after every trial and then is low pass filtered by a process with a time constant = 2. The numerical integration time step was 0.1 ms. A minimum weight of 0.1 was set for all weights.
A saturation to the excitatory and inhibitory firing rate (100 and 250 Hz, respectively) was added to prevent the nonbiological scenario in which activity could diverge towards infinity under unstable conditions. Note that at the fixed point the saturation is not necessary for the cross-homeostatic rule because it is inherently stable as proved in the Supplementary Material (Section 1.3).
In Fig. 2D and 4D-G we initialize the weights uniformly between the following ranges: [4,7], [0  Figure S3). We emphasize that many of these initializations resulted in starting conditions with exploding network rates, which were held in check by the saturation imposed on the transfer function. Despite this initial instability, the rules successfully brought the rates to the setpoint values.

Analytical stability analyses of the firing rate model
We analyzed the entire dynamical system (composed of the neural subsystem and the learning rule subsystem) for every synaptic learning rule considered in this work, and analyzed its stability. In every case, the general prescription is: a) Take the combined neural and learning rule subsystems and nondimensionalize all variables, so that the two different time scales are evident (fast neural, slow synaptic plasticity). For the description of the learning rule subsystem we switch from discrete-time dynamics to continuous-time dynamics: ∆W → τ0 dW/dt b) Make a quasi-steady state (QSS) approximation of the neural subsystem.
This means we will consider the neural subsystem is fast enough so that it converges "instantaneously" (when compared to the synaptic plasticity subsystem) to its corresponding fixed point. For this we will require that the stability conditions of the neural subsystem are satisfied (see below). c) Find the steady-state solution of the synaptic plasticity subsystem, i.e. the self-sustained activity fixed point; compute the Jacobian of the synaptic plasticity subsystem at the fixed point; compute the eigenvalues of the Jacobian. Two out of the four eigenvalues are expected to be zero because the solution is not an isolated fixed point of the system but a continuous 2D plane in 4D weight space. d) Address (linear) stability. If both nonzero eigenvalues have negative real parts, then the fixed point is stable under the learning rule; if at least one of the nonzero eigenvalues has positive real part, then the fixed point is unstable. (A note on abuse of notation: we might say indistinctly "the fixed point is stable/unstable" and "the learning rule is stable/unstable".) For a detailed explanation see Section 2 in the Supplementary Material.

Implementation of the Multi-unit firing rate model
A rate-based recurrent network model containing Ne = 80 excitatory and Ni = 20 inhibitory neurons was implemented with all-to-all connectivity (without selfconnections). The activation of the neurons followed equations (1), (2) and (3) in the main Methods. The same parameters as for the population model were used, where WXY represents now a matrix of synaptic weights from population X to population Y. A minimum weight of 0.1/Nx for WEI and WIE and 0.1/(Nx-1) for WEE and WII was set for all weights.
The synaptic plasticity rules were implemented as follows.
Cross-homeostatic family of rules: where i and j represent the post-and presynaptic neurons, respectively, and k denotes the presynaptic inhibitory neurons targeting the excitatory neurons (or the presynaptic excitatory neurons targeting an inhibitory neuron). NE and NI denote the total number of excitatory and inhibitory neurons, respectively. The weights are therefore updated following the average presynaptic error of the crossed E/I population classes. Note as stated above that this formulation can be implemented in a local manner (see Discussion). A learning rate of = 0.00002 was used for all simulations.
Two-term cross-homeostatic family of rules: here the first term represents the standard homeostatic rule, and the second term cross-homeostatic plasticity (as implemented above). A learning rate of = 0.00001 was used for all simulations.
In Fig. 5G-H and 6G-H we initialize the mean weights of the population uniformly in between the following ranges: [1,6], [0.5,2], [5,7], [0. 5,2]. The weights within each class were then normally distributed around that mean (and normalized by the number of neurons) with a standard deviation of 10% of the mean. Note that this initialization led to multiple initial conditions with exploding network rates (which were held in check by the saturation cutoff of the neurons). Those initial rates are not displayed in Fig. 5-6H for visualization purposes, but the rules successfully brought all those cases to the corresponding setpoints (final rates are displayed). Simulations were run for 1000 trials to assess stability of the convergence. Note that broader initializations for all weights [0,10], [0,10], [0,10], [0,10] also result in convergence of the rules, although many more of those combinations display initial exploding network rates (Supplementary Fig. S7). In the example shown in Fig. 5A-F and 6A-F individual weights were normally distributed with equal mean across weights classes and standard deviation (0.1±0.04). For Supplementary Fig. S5G-H a lognormal distribution of weights was used. The mean weights were distributed uniformly as in Figure 5G-H, and then weights within each class were distributed log-normally, with arithmetic standard deviation of 0.05. We note that while increasing the standard deviation of the log-normal led to more complex time-varying neural dynamics (1), the rules still converged to the average set-points.

Implementation of the Spiking model
The spiking model was designed based on previous work (2). Units in the model were leaky integrate-and-fire neurons with a spike adaptation current. The membrane potential of each unit was represented as The noise term √ ( ) represents an Ornstein-Uhlenbeck process with zero mean, standard deviation , and a time constant equal to the membrane time constant = / . When ( ) ≥ ℎ ℎ , the unit emitted a spike, its voltage was reset to , and its adaptation current was incremented by / . After spiking, the unit entered an absolute refractory period . Default values for unit parameters can be found in Table S1.
Synapses were current-based, and the total synaptic current ( ) was summed across each unit's incoming synapses with distinct synaptic weights determined by the matrices , , , and . Total synaptic current to a postsynaptic excitatory or inhibitory unit was given by each of the following two equations, respectively The kinetics of the synaptic currents were determined by the function ( , , ) for each presynaptic unit y and postsynaptic unit x. When a presynaptic spike occurred in unit y at time * , ( , , ) was incremented by an amount described by a delayed difference of exponentials equation (3) where indicates the postsynaptic membrane time constant. Thus, synaptic kinetics were determined by the delay , the rise time , and the decay time .
The synaptic delay was uniformly distributed between 0 and 1 ms across all excitatory (inhibitory) synapses. Synaptic parameters can be found in Table S2.
Networks consisted of 1600 E units and 400 I units with probability of connection = 0.25 and no autapses (self-connections). Connectivity was uniformly random, and weights for each synaptic class were initialized from normal distributions with a coefficient of variation equal to 0.2. For the example training session shown in Figure 7, the initial weights were ̅̅̅̅̅̅ = 80 , ̅̅̅̅̅ = 100 , ̅̅̅̅̅ = 350 , ̅̅̅̅̅ = 225 . Network simulations were evaluated using forward Euler integration using a time step of 0.1 ms.
During each trial (1.5 s) of training a brief external current large enough to cause a spike ( ⇒ + 0.98 ) was injected into 100 E units. This constituted a "kick" (2, 4) that provided the possibility for recurrent excitation to ignite a self-sustaining Up-state. After each trial, the contiguous time period of nonzero FR was detected, and each unit's FR during that time period was calculated. FRs for each unit in each trial contributed to a moving average vector with a time constant of 2 trials, which we refer to as ⃗ and ⃗ ℎ . Accordingly, we refer to the firing rate setpoints as and ℎ , which were set to 5 and 14 Hz respectively, just as in the firing rate model.
The combined plasticity equations used in the spiking model were each formulated as a sum of homeostatic (first) and local cross-homeostatic (second) terms where is a learning rate constant set to 0.0025 2 . For the local crosshomeostatic term, each element of ⃗ represents the average FR of the units in the opposite population that synapse onto that unit. This is calculated by multiplying the unit FRs by the connectivity matrix and dividing by the vector that results from summing its columns, which we refer to as ⃗ . Note that the ⊘ symbol refers to element-wise division For a vector of presynaptic FRs ( ⃗ , ⃗ ℎ ) we imposed a minimum value such that each element was at least 1 Hz (otherwise networks can get stuck in the downstate). Additionally, all synaptic weights were bounded to stay within minimum and maximum weight values of 10 pA and 750 pA respectively.
For the paradoxical effect analysis in Figure 7I, the adaptation current was disabled for all units in order to allow for a long and stable Up-state. In each trial, a kick was given at 100 ms. From 3 to 4 s, a small positive current was injected into all inhibitory units. 40 trials were conducted at each value of the injected current, and a PSTH for each value was constructed using the inhibitory population spiking activity across all trials.
For the analysis in Figure 7L, nine networks were trained from distinct mean weight values for the four synaptic classes as shown in Table S3. We measured the error of unit firing rates with respect to their setpoints at the beginning and end of training, as quantified by the mean-squared error (MSE) across all individual excitatory and inhibitory units at each trial where ⃗ represented the FR of a unit in that population (excitatory or inhibitory) and represented the corresponding set-point.

Summary of results
In this section we describe the general results of the analytical stability analyses of the joint neural and synaptic plasticity subsystems. We express results in terms of the "free weights" W EE and W IE . Subscript "up" is used to identify values at the nontrivial fixed point where E and I are larger than zero (as opposed to "down" where E = I = 0 which is the other possible solution). In Section 2 we provide a detailed description of the approach.

Homeostatic plasticity
In continuous-time dynamics, the equations for the Homeostatic plasticity rule are The condition for the fixed point to be stable (i.e., the two nonzero eigenvalues to have negative real parts, see Section 2) under this rule is: It is difficult to determine whether the stability condition of Eq. 2 is satisfied for a general set of parameter values (see numerical analysis below). However, this condition can be re-expressed in a more useful form in terms of W EE and W II : where R = E set /I set As an illustration of the results above, in Figure  S8(top) we plot the stability condition Eq. 2 with parameter values as in Table 1 and learning rates α XY = 0.02. It is clear that the plasticity rule is stable in a region with little overlap with the stability region of the neural subsystem. The stability region can be enlarged by making the dynamics of the weights onto the inhibitory neuron slower, as in Figure S8(bottom) where α EE = α EI = 0.02 and α IE = α II = 0.0002. See Section 2.3 for a detailed analysis.

Homeo-antiHomeo variations
The stability condition in the previous section was obtained by assuming all learning rates are positive. Interestingly, if some of them are negative then the fixed point may still be stable. A negative learning rate can be interpreted as the corresponding equation being anti-homeostatic, i.e. if the neural activity (E or I) departs from its setpoint then the rule will drive it even farther away. While this kind of behavior would be usually deemed undesired, it is worth considering due to its relationship with the paradoxical regime.
In this section we consider the Homeostatic rule, Eq. 28, and let the learning rates α XY be either positive or negative. The particular case where all learning rates are positive corresponds to the original Homeostatic plasticity rule.
Once we free the signs of the learning rates, the fixed point needs two conditions to be stable:

Cross-Homeostatic plasticity
In continuous-time dynamics, the equations for the Cross-Homeostatic plasticity rule are and its stability condition in terms of the free weights W EE and W IE reads: This stability condition can be put in a simpler form by switching to W EI and W IE :

Two-Term plasticity
The equations for the Two-Term plasticity rule in continuous-time dynamics are and its stability condition in terms of the free weights In Figure S10 we plot the stability condition of this rule, Eq. 10, for three different parameter values: α β (the "Cross-Homeostatic" terms dominate over the "Homeostatic" terms, and the rule is stable with the largest stability region); α = β (the two terms are of comparable size); and α β (the "Homeostatic" terms dominate instead, and the stability region of the rule is as small as that of the Homeostatic plasticity).
In order to determine the validity of the stability condition, Eq. 10, in a more general situation, we rewrite it in a more useful form: Note that the following is satisfied for a biologically backed set of parameter values: and thus it is likely that c > 0. In addition, b and b are positive definite, and a, a > 0 in the paradoxical regime (W EE g E − 1 > 0). All this makes the stability condition likely satisfied, and thus the plasticity rule stable. Finally, a small enough β would make the condition more likely to hold.

SynapticScaling plasticity
The equations for the SynapticScaling plasticity rule in continuous-time dynamics are and the condition for the fixed point to be stable under this rule is (where α 2,3,4 are defined as in previous subsections). This stability condition does not hold for biologically backed parameter values unless the dynamics of the weights onto the inhibitory neuron are slow enough (and in a few fine-tuned cases). To show this, we express the stability condition in terms of the free weights W EE and W IE and plot it with parameter values as in Table 1 and equal rates (α XY = 0.02; Figure S11 top left). The stability condition is a homographic function (i.e. a hyperbola) with stability regions in its upper-left and lower-right quadrantsentirely outside the stability region of the neural sub-system. If the dynamics of the weights onto the excitatory neuron are made slower, the homographic function is even steeper (bottom left); if the weights onto the inhibitory neuron are made slower instead, the stability regions switch and overlap with the stability region of the neural subsystem, making the fixed point stable (top right). It is illustrative to consider the particular case where all learning rates are equal. In this case the stability condition, Eq. 13, doesn't depend on the learning rates and takes the simpler form: Note that in a biologically backed set of parameter values the following is true: This makes b > 0 and likely a > a (in addition, a is a sum of positive terms while a is a difference). Then in the paradoxical regime (W EE g E − 1 > 0) it seems likely that the stability condition is not satisfied, because the right-hand side is a sum of positive terms and one of them is likely greater than the left-hand side. The SynapticScaling rule is then likely unstable when the learning rates are equal. A more general case with different learning rates can be analyzed by grouping terms in the following way: If (W EE g E − 1) > 0 (paradoxical regime), then decreasing α 3 and/or α 4 (slow dynamics of the weights onto the inhibitory neuron) helps satisfying the condition and thus making the rule stable.
See Section 2.4 for a detailed analysis.

ForcedBalance plasticity
The equations for the ForcedBalance plasticity rule are and the conditions for the fixed point to be stable under this rule are where In Figure S12 we plot the stability condition of this rule, Eq. 16, for three different parameter values: α 1 = α 3 , α 1 α 3 (inhibitory plasticity slower); and α 1 α 3 (excitatory plasticity slower). In order to decide whether conditions Eq. 16 are satisfied in a more general case, note that b 1 and b 2 on the left-hand side are subtractions whereas b 1 and b 2 on the right-hand side are sums of positive definite terms, which helps satisfying the condition. On the other hand, one of the stability conditions of the neural subsystem might counter the effect: (W II up g I + 1)τ E > (W EE up g E − 1)τ I (see Section 2.2 below) but for biologically backed parameter values it is τ E > τ I thus leaving room for the condition to hold. See Section 2.4 for a detailed analysis.    2. Make a quasi-steady state (QSS) approximation of the neural subsystem [1,2]. This means we will consider the neural subsystem is fast enough so that it converges "instantaneously" (when compared to the synaptic subsystem) to its cor-responding fixed point. For this we will require that the stability conditions of the neural subsystem are satisfied (see below).
3. Find the steady-state solution of the synaptic subsystem, i.e. the fixed point; compute the Jacobian of the synaptic subsystem at the fixed point; compute the eigenvalues of the Jacobian [2,3]. Two out of the four eigenvalues are expected to be zero because the fixed point is not an isolated fixed point of the system but a continuous 2D plane in 4D weight space.
4. Address (linear) stability. If both nonzero eigenvalues have negative real part, then the fixed point is stable under the plasticity rule; if at least one of the nonzero eigenvalues has positive real part, then the fixed point is unstable [2,3]. (A note on abuse of notation: we might say indis-tinctly "the fixed point is stable/unstable" and "the plasticity rule is stable/unstable".) Eigenvalues and stability in the presence of continuous, i.e. non-isolated, attractors have been discussed in the context of neural networks for eye position control [4,5] (keep in mind that the eigenvalues' critical value in these references is 1 instead of zero because they consider eigenvalues of the connectivity matrix alone, whereas we consider eigenvalues of the whole linear part). As the fixed point is a collection of non-isolated fixed points that form a 2D plane, there is no dynamics along the plane, and the linear stability analysis is enough to fully address stability-we do have two zero eigenvalues, but there is no need to compute the center manifold [3] because the other two eigenvalues represent the whole dynamics around the fixed point and have nonzero real part.
In order to apply the tools from Dynamical Systems' theory for flows in a unified way for both the neural and synaptic subsystems, we will switch from a discrete-time description of synaptic weight dynamics (where the change in weight W is represented by ∆W applied every certain time interval) to a continuoustime description (where the weights are continuously evolving albeit with a long time scale τ 0 ): In the following we first define the neural subsystem and compute its stability conditions (next subsection). Then we consider every plasticity rule in detail (following subsections).
Paradoxical regime. In this text we show detailed calculations of the stability conditions for the Homeostatic plasticity in the paradoxical regime only; see Section 2.5 for the non-paradoxical case.

Neural dynamics
For the neural+synaptic system in the QSS approximation to be stable under a specific synaptic plasticity rule, it is necessary that the neural subsystem is stable so it remains in its QSS solution as the weights evolve. In this section we define the neural subsystem and compute its stability conditions.

System's equations and fixed points
We consider a two-subpopulation model with firingrate units E and I with ReLU activation functions (gain g X , threshold Θ X , with X = E, I). The dynamics for synaptic currents above threshold is given by: All variables and parameters are definite positive. In this subsection the synaptic weights W XY are fixed.
Neural fixed point. The fixed point of the neural subsystem in the supreathreshold regime is the solution of dE/dt = dI/dt = 0: We named it with the subscript "up" to distinguish it from the trivial solution "down" where E and I are zero (and the neural subsystem is below threshold). The activity of the excitatory and inhibitory subpopulations at the nontrivial fixed point, E up and I up , depend on all weight values. Only some of the combinations, however, lead to a stable steady state. We compute the stability conditions in the following subsection.

Stability of the nontrivial neural fixed point
The Jacobian matrix, that is the matrix of first derivatives, gives information regarding the stability 8 of fixed points: if the real parts of its eigenvalues are all negative, then the fixed point is stable. The Jacobian of the neural system (Eq. 17) is 20) Its eigenvalues can be expressed as: where Tr and Det are the trace and determinant of the matrix, respectively. For eigenvalues either complex or purely real, their real parts are negative (and thus the fixed point is stable) when Det > 0 and Tr < 0, that is: Note that the positive determinant condition, Eq. 22, is equivalent to C > 0 (Eq. 19).
In the following, we will require that the stability conditions of the neural subsystem, Eqs. 22 and 23, are satisfied.

Weight values consistent with the neural fixed point
The fixed point relationships, Eq. 18, are expressed as the E and I values resulting from a given set of weight values. If we set instead E and I to their target values E set and I set and solve for the weights, we get the weight values that are consistent with a given fixed point activity: Note first that any stable plasticity rule for the evolution of the weights for the neural subsystem (Eq. 17) must converge to weight values in accordance with these relationships (either in the form Eq. 24 or Eq. 18).
Second, note that the system is underdetermined and that is why two of the weights are free (chosen to be W EE and W EI ). Note also that all weight values must be positive; specifically, requiring W EI up > 0 and W II up > 0 leads to We refer to these expressions as the "positive W EI " and the "positive W II " conditions, respectively.

Paradoxical effect
The paradoxical effect arises when an external depolarization of the inhibitory subpopulation (increase of I) produces an actual decrease of I. In this model, an external depolarization of I can be mimicked by a decrease of its threshold Θ I , thus there is a paradoxical effect whenever the coefficient of Θ I in the numerator of I up is positive. The coefficient is g I (W EE g E −1)/C and thus there is paradoxical effect if The paradoxical effect can also be seen in a plot of the fixed point values E up and I up (Eq. 18) as a function of each individual weight. Specifically, from a naive point of view I up should increase when W IE is increased, and decrease when W II is increased; however, it does the opposite in either case (see Figure  S13).  Figure S13: Paradoxical effect in the neural subsystem (W EE = 5; parameter values as in Table 1). E up behaves as expected when each weight is varied. I up , however, shows paradoxical behavior when either W IE or W II are varied. Dashed lines are the vertical asymptote of every case.

Homeostatic plasticity: Detailed calculation
In this section we show in detail the calculation of the stability condition for the Homeostatic plasticity rule.

Definition of the plasticity rule
In continuous-time dynamics, the Homeostatic plasticity rule reads: where α XY (X, Y = E, I) are the learning rates (with appropriate units) setting the time scales of the weight dynamics, and E set and I set are the set points of the excitatory and inhibitory subpopulations, respectively.
The fixed points of the system (i.e. steady states) are determined by setting all derivatives to zero. There is a non-trivial fixed point compatible with the neural subsystem being above threshold: it is the set of weight values such that: The values of the weights corresponding to the nontrivial neural fixed point are given by the (underdetermined) system defined by equating Eqs. 29 and 18. Since it is a two-equation system for a set of four unknown weights, there are two free weights that we choose to be W EE up and W IE up . The values of the other two are given by Eq. 24. This means that the fixed point is actually a continuous set of nonisolated fixed points forming a 2D plane in 4D weight space. In other words, there is an infinite number of weight values compatible with the nontrivial neural fixed point (possibly not all stable, though).

Nondimensionalization
Next we nondimensionalize all variables in order to have a simpler system and make the QSS approximation in a safe way. We define new (nondimensional) variables e, i, τ , w EE , w EI , w IE , and w II , and their corresponding scaling parameters. We substitute the new variables into the full system (neural+synaptic, Eqs. 17 and 28) and choose the values of the scaling parameters such that all nondimensional variables are of order 1 (see attached SageMath code). With this, the full system reads: where we defined the new parameters

Quasi-steady state approximation
Neural dynamics evolves in a much shorter time scale (τ E and τ I ) than synaptic dynamics (τ 0 ): thus we can safely assume e and i very quickly reach quasi-equilibrium values, i.e. practically instantaneous convergence to quasi-steady state (QSS) values as if the weights were fixed, while the synaptic weights evolve according to their slow dynamics. This allows us to reduce the system's dimensionality from six to four.
In the QSS approximation, the values of the nondimensionalized excitatory and inhibitory activities instantaneously track the slow dynamics of the plasticity rule. They are determined by applying Eq. 31 to the first two rows of Eq. 30; solving for e and i leads to The full system in the QSS approximation reads where e qss and i qss are nonlinear functions of the weights as defined by Eq. 32. Note that the nontrivial neural fixed point, defined by making all derivatives equal to zero, can be expressed as e qss = 1 which is the nondimensionalized version of Eq. 29.
The weight values compatible with this condition are defined by equating Eqs. 32 and 34: (w EE up and w IE up are free). This is the nondimensionalized version of Eq. 24.

Stability condition
The program for assessing linear stability of the fixed point is as follows: a) compute the Jacobian (the matrix of first derivatives) of Eq. 33 and evaluate it at the fixed point; b) compute the eigenvalues of the Jacobian (two of them will be zero because the fixed points form a continuous 2D plane in phase space); c) If the real part of the two nonzero eigenvalues is negative then the fixed point is stable; if at least one of the nonzero eigenvalue has positive real part then the fixed point is unstable.
Jacobian matrix. Let the full system in the QSS approximation (Eq. 33) be written as where e qss and i qss are functions of the weights as defined by Eq. 32. By applying the chain rule the elements J ij (i, j = 1 . . . 4) of the Jacobian matrix can be expressed as In order to have the Jacobian specialized in the fixed point, these expressions are to be substituted by Eqs. 32-35.
Eigenvalues of the Jacobian matrix. The Jacobian matrix has two zero eigenvalues and two nonzero eigenvalues. The nonzero eigenvalues have the form: Sign of the eigenvalues. To determine the sign of the real part of Eq. 36, first note that the factor D is positive definite. Second, C must be positive because it is related to one of the stability conditions of the neural subsystem (Eq. 22, after substituting back to dimensionalized quantities). Note next that A 2 − DC is less than A 2 (since C and D are positive), and thus the square root is either real and less than |A| or imaginary, both cases leading to The plasticity rule is then stable (both eigenvalues have negative real part) if A < 0, which in terms of the original parameters and free weights W EE and W IE reads:

Analysis of the stability condition
It is hard to determine whether the stability condition Eq. 38 is satisfied for a general set of parameter values (see numerical analysis below). However, by using the fixed point relationship Eq. 24, this condition can be re-expressed in a more useful form in terms of W EE and W II : Note that learning rates values of the same order lead to α 2,3,4 ∼ 1 and that biologically backed parameter values satisfy: both likely preventing the condition to hold.
On the other hand, small enough values of α 3 and α 4 (by making the dynamics of the weights onto the inhibitory neuron W IE and W II slower) would help satisfy the condition thus making the system stable.

Relationship between the synaptic stability and the paradoxical condition
The boundary of the stability condition for this plasticity rule, Eq. 38, is a linear function in the (W EE , W IE ) space with a slope that tends to infinity as the excitatory learning rates (α EE ,EI ) tend to zero: set α EE + I 2 set α EI )E set while its root is a complicated expression (see Sage-Math notebook) that tends to W EE = 1/g E . The region of stability is to the left of the line. Thus, the boundary of stability in this limit coincides exactly with the boundary of the paradoxical condition (W EE > 1/g E ). This can be construed as an inconsistency/contradiction between the stability of the rule and the existence of the paradoxical effect.

Detailed calculations for the other rules
The stability calculations for the rest of the rules follow very similar paths. They can be found in the corresponding SageMath-Jupyter notebooks: upstates-CrossHomeostatic stability.ipynb upstates-TwoTerm stability.ipynb upstates-SynapticScaling stability.ipynb upstates-ForcedBalance stability.ipynb

Stability of the rules in a nonparadoxical regime
All results above were developed with the neural subsystem set in the paradoxical regime-that is, the region in (W EE , W IE ) leading to a stable fixed point was completely within the paradoxical region (W EE g E > 1). In order to show the importance of the paradoxical behavior for the stability of the plasticity rules, we also computed the stability conditions of every plasticity rule in a more general setting where the excitatory subpopulation in the neural subsystem has an external, constant, excitatory input current I ext . This allows the neural subsystem to display both paradoxical and non-paradoxical stable behavior (in the second case, at the expense of the fixed point not being an inhibition-stabilized fixed point; see upstates-Neural subsystem stability-with Iext.ipynb).

Homeostatic with I ext
The stability condition doesn't depend on I ext and it reads the same as Eq. 2: (SageMath code in upstates-Homeostatic stability-with Iext.ipynb)

CrossHomeostatic with I ext
The stability condition with I ext is: (41) which is very similar to Eq. 7 except that it has (Θ E − I ext ) instead of just Θ E . From this it should be evident that the condition will still hold for any positive value of I ext (right-hand side decreases).
The validity of the condition can also be seen after switching to W IE and W EI , leading to exactly the same condition as Eq. 8: which holds for any value of I ext .

TwoTerm with I ext
The stability condition with I ext is: (43) which is very similar to Eq. 10 except that it has (Θ E − I ext ) instead of just Θ E . From this it should be evident that the larger the value of I ext (righthand side decreases) the larger the stability region.

SynapticScaling with I ext
When I ext is included in the dynamics of E, the stability condition for the rule reads: which is very similar to Eq. 13 except that it has (Θ E − I ext ) instead of just Θ E . From this it should be evident that including a positive I ext will increase the chances that the condition holds (right-hand side increases).

ForcedBalance with I ext
The stability conditions when I ext is included in the neural subsystem are: where Here we show how to compute a plasticity rule for the weights starting from a loss function. Then we make an approximation by considering that the weight values are close to the values corresponding to the fixed point.

General prescription
We consider the full neural+synaptic system in the QSS approximation (see e.g. Section 2.3). In this approximation the neural subsystem is represented by the quasi-steady-state values where the functions E up and I up are defined by Eq. 18 (see [7] for a related discussion on quasi-steady state, synaptic plasticity, and gradient descent). The synaptic subsystem, that is the plasticity rule, will be obtained as a result of considering a specific loss function, and the general prescription to compute the plasticity rule from a loss function L is the following: 1. Consider a loss function depending on E and I (which in turn depend on all weights): Conditions to be satisfied by the loss function are, for instance, to be smooth enough (i.e. continuous and differentiable) and to have a minimum when the activities E and I are at the set points E set and I set (i.e. homeostatic plasticity).
2. The dynamics of the weights is such that it follows a gradient descent on the loss function towards its minimum. In vector notation: with a single learning rate α for simplicity. The unfolded plasticity rules, that is the equations that govern the weights' dynamics, are then to be taken from the quasi-steady-state values of E and I, Eq. 46. We will, however, compute the partial derivatives from the implicit expressions given by setting dE/dt = dI/dt = 0 in Eq. 17 without solving for E and I.

Exact plasticity rules
Loss function. We choose a very general loss function that depends homeostatically on both E and I activities: This loss function is an elliptic paraboloid in (E, I) space with a global minimum at (E set , I set ) so a gradient descend working on E and I should converge to that minimum (see Liapunov function and gradient systems: [3, Section 1.1B][8, Sections 9.3 and 9.4][2, Section 7.2]). Keep in mind, however, that L has a different shape when expressed as a function of the weights, and that E and I are not necessarily monotonic functions of the weights (particularly for a paradoxical system), so the conditions for the set point of L to be stable or a global minimum or even unique are not necessarily satisfied.
Partial derivatives of L. The partial derivatives of L with respect to E and I are simply Partial derivatives of E and I. We compute the partial derivatives ∂X/∂W XY (X, Y = E, I) by first equating the neural subsystem (Eq. 17) to zero: where C = W EI W IE g E g I − (W II g I + 1)(W EE g E − 1) Exact plasticity rules. Putting everything together, the plasticity rules Eq. 48 are: ∆W EE = − α C ((I set − I)EW IE g e g I + (E set − E)E(W II g I + 1)g E ) ∆W EI = + α C ((I set − I)IW IE g e g I + (E set − E)I(W II g I + 1)g E ) ∆W IE = + α C ((E set − E)EW EI g e g I + (I set − I)E(W EE g E − 1)g I ) ∆W II = − α C ((E set − E)IW EI g e g I + (I set − I)I(W EE g E − 1)g I ) (56) Note that these are very complicated, nonlinear expressions because both E and I depend on all weights via Eq. 53. Also the denominator C depends on all weights (see previous paragraph).

Approximation
We want simpler expressions for the plasticity rules. Note that the exact expressions above all have a homeostatic factor (either E − E set or I − I set ) and a presynaptic factor (either E or I), while the rest are complicated expressions coming from the derivatives ∂E/∂W XY and ∂I/∂W XY . We want to keep the homeostatic and presynaptic factors as they are while simplifying the rest of the expressions (explicit dependence on the weights including C) by performing a lowest-order Taylor series expansion of the explicit dependence of Eqs. 55 on the weights. Although this is not a textbook Taylor expansion of the full expressions, it is very informative because the results can be much more easily interpreted (for a similar approach see [7]).
We perform a zeroth-order approximation of the derivatives ∂E/∂W XY and ∂I/∂W XY as functions of the weights (i.e. while holding the presynaptic factors E and I constant) around the fixed point. In this approximation the weights are not small but close to their target values, represented by the relationships Eq. 24. By substituting the result in Eq. 49, we get the following approximated plasticity rules: where Analysis. Note that α E , α I , β E , and β I are all constant. Furthermore, note that • A > 0 as it is equal to the "positive W EI " condition, Eq. 25; • B > 0 as it is part of the "positive W II " condition, Eq. 26; • D > 0 as it is equal to the numerator of I up , Eq. 18 (up to a positive factor), which must be positive because the denominator is.
Interestingly, note that the learning rate β I can be either negative or positive depending on whether the fixed point where the dynamics is converging to is paradoxical (W EE up g E − 1 > 0) or not (W EE up g E − 1 < 0). Note that the terms with α E,I in the approximated plasticity rules, Eq. 57, are exactly equal to the Cross-Homeostatic rules, Eq. 6. Additionaly, the terms with β E,I are exactly equal to the Homeostatic rules, Eq. 28, unless β I < 0 which would make the plasticity rule a Cross-Homeo-antiHomeo hybrid.