Stable Control of Firing Rate Mean and Variance by Dual Homeostatic Mechanisms

Homeostatic processes that provide negative feedback to regulate neuronal firing rates are essential for normal brain function. Indeed, multiple parameters of individual neurons, including the scale of afferent synapse strengths and the densities of specific ion channels, have been observed to change on homeostatic time scales to oppose the effects of chronic changes in synaptic input. This raises the question of whether these processes are controlled by a single slow feedback variable or multiple slow variables. A single homeostatic process providing negative feedback to a neuron’s firing rate naturally maintains a stable homeostatic equilibrium with a characteristic mean firing rate; but the conditions under which multiple slow feedbacks produce a stable homeostatic equilibrium have not yet been explored. Here we study a highly general model of homeostatic firing rate control in which two slow variables provide negative feedback to drive a firing rate toward two different target rates. Using dynamical systems techniques, we show that such a control system can be used to stably maintain a neuron’s characteristic firing rate mean and variance in the face of perturbations, and we derive conditions under which this happens. We also derive expressions that clarify the relationship between the homeostatic firing rate targets and the resulting stable firing rate mean and variance. We provide specific examples of neuronal systems that can be effectively regulated by dual homeostasis. One of these examples is a recurrent excitatory network, which a dual feedback system can robustly tune to serve as an integrator.


Introduction
Homeostasis, the collection of slow feedback processes by which a living organism counteracts the effects of external perturbations to maintain a viable state, is a topic of great interest to biologists [1,2]. The brain in particular requires a precise balance of numerous state variables to remain properly operational, so it is no surprise that multiple homeostatic processes have been identified in neural tissues [3]. Some of these processes appear to act at the level of individual neurons to maintain a desirable rate of spiking. When chronic changes in input statistics dramatically lower a neuron's firing rate, multiple slow changes take place that each act to increase the firing rate again, including the collective scaling of afferent synapses [4,5] and the adjustment of intrinsic neuronal excitability through adding and removing ion channels [5][6][7]. These changes suggest the existence of multiple independent slowly-adapting variables that each integrate firing rate over time and provide negative feedback.
Here we undertake an analytical investigation of the dynamics of homeostasis via two independent slow mechanisms ("dual homeostasis"). Our focus on dual homeostasis is partially motivated by the rough breakdown of firing rate homeostatic mechanisms into two categories, synaptic and intrinsic. In our analytical work, we maintain sufficient generality to describe a broad class of firing rate control mechanisms, but we illustrate our results using examples in which homeostasis is governed by one synaptic mechanism acting multiplicatively on neuronal inputs and one intrinsic mechanism acting additively to increase or decrease neuronal excitability. We limit our scope to dual homeostasis to allow us to derive strong analytical results.
It is not immediately clear that dual homeostasis should even be possible. When two variables independently provide negative feedback to drive the same signal toward different targets, one possible outcome is "wind-up" [2], where each variable perpetually ramps up or down in competition with the other to drive the signal toward its own target.
In a recent publication [8], we perform numerical simulations of dual homeostasis (intrinsic and synaptic) in biophysically detailed neurons. We show empirically that this dual homeostasis is stable across a broad swath of parameter space and that it serves to restore not only a characteristic mean firing rate but also a characteristic firing rate variance after perturbations.
Here, we demonstrate analytically that stable homeostasis occurs in a broad family of dual control systems. Further, we find that dual homeostatic control naturally preserves both the mean and the variance of the firing rate, a task impossible for a homeostatic system with a single slow feedback mechanism. We identify broad conditions under which a dually homeostatic neuron possesses a stable homeostatic fixed point, and we derive estimates of the characteristic firing rate mean and variance maintained by homeostasis in terms of homeostatic parameters. We use rate-based neurons and Poisson-spiking neurons for illustration, but our main result is sufficiently general to apply to any model neuron.
One specific application in which such a control system could play an essential role is in tuning a recurrent excitatory network to serve as an integrator. This task is generally considered one that requires biologically implausible precise calibration of multiple parameters [9,10] and is not well understood (though various solutions to the fine tuning problem have been proposed in [11][12][13][14]). In [8], we show empirically that a heterogeneous network of dually homeostatic neurons can tune itself to serve as an integrator. Here, we demonstrate analytically in a simple model of a recurrent excitatory network that integrating behavior can be stabilized by single-cell dual homeostasis and that this stability is robust to the homeostatic parameters of the neurons in the network.
In Sect. 2, we introduce our generalized model of dual homeostasis with the simple but informative example of synaptic/intrinsic firing rate control, and we discuss the reasons that stable homeostatic control is possible for this system. In Sect. 3, we pose a highly general mathematical model of dual homeostatic control. We derive an estimate of the firing rate mean and variance that characterize the fixed points in a given dual homeostatic control system and conditions under which these fixed points are stable. In Sect. 4, we give further specific examples that are encompassed by our general result. In Sect. 5, we use our results to explore dual homeostasis as a strategy for integrator tuning in recurrent excitatory networks. In Sect. 6, we summarize and discuss our results.

Preliminary Examples
Throughout this manuscript, we consider a homeostatic neuronal firing rate control system with slow homeostatic control variables that serve as parameters for neuronal dynamics. Each of these variables represents a biological mechanism that provides slow negative feedback in response to a more rapidly varying neuronal firing rate r. In this section, we focus on an example in which the two control variables are (1) a factor g describing the collective scaling of the strengths of afferent synapses and (2) the neuron's "excitability" x, which represents a horizontal shift in the mapping from input current to firing rate. An increase in x can be understood as an increase in excitability (or a decrease in firing threshold) and might be implemented in vivo by a change in ion channel density as suggested in [7]. The choice of x and g as homeostatic control variables is motivated by the observation that synaptic scaling and homeostasis of intrinsic excitability operate concurrently in mammalian cortex [5]. We write this dual control system in the form where r is a neuronal firing rate, r x and r g are the "target firing rates" of the two homeostatic mechanisms, f x and f g are increasing functions describing the effect of deviations from the target rates on the two control variables, and τ x and τ g are time constants assumed to be long on the time scale of variation of r. An extra factor of g multiplies the second ODE because g acts as a multiplier and must remain nonnegative. As a result, g increases/decreases exponentially (or ln(g) increases/decreases linearly) if the firing rate r is below/above the target rate r g . This extra g multiplier is not essential for any of the results we derive here.
In general, r may represent the firing rate of any type of model neuron, or a correlate of firing rate such as calcium concentration. Likewise, the target rates r x and r g may represent firing rates or calcium levels at which the corresponding homeostatic mechanisms equilibrate. We assume that r changes quickly relative to τ x and τ g and that r is "distribution-ergodic." This term is defined precisely in the next section; intuitively, it means that over a sufficiently long time, r behaves like a series of independent samples from a stationary distribution. This allows us to approximate the right-hand sides in (1) by averages over the distributions of f x (r) and f g (r). We will use · to represent the mean of a stationary distribution. Since the dynamics of the firing rate depends on control variables x and g, the distributions we consider here also depend on these variables. Averaging (1) over r, we can write In this section, we assume that the neuron is a standard linear firing rate unit with time constant τ r receiving synaptic input I (t). This input is scaled by synaptic strength g, and the neuron's response is additively shifted by the excitability x. Thus, the firing rate is described by the equation In a later section, we will consider a similar system with spiking dynamics.

Constant Input
First, let us assume that I (t) is equal to a constant φ. In this case, r assumes its asymptotic value r = gφ + x on time scale τ r and closely tracks this value as g and x change slowly. Thus, we have f x (r) = f x (gφ + x) and f g (r) = f g (gφ + x). To find the x-nullcline (the set of points (x, g) whereẋ = 0), we setẋ = 0 in (2). Since f x is increasing, it is invertible over its range, so we find that the x-nullcline in the (x, g) phase plane consists of the set gφ + x = r x . Similarly, the g-nullcline is the line gφ + x = r g , plus the set g = 0. Fixed points of this ODE are precisely the set of intersections of the nullclines. We are interested primarily in fixed points with g > 0, so we ignore the set g = 0. Representative vector fields, nullclines, and trajectories for (2) with r = gφ + x are illustrated in Fig. 1. From the nullcline equations, it is clear that if r x = r g , there are no fixed points with g > 0. If r g > r x , g increases and x decreases without bound (Fig. 1A); if r g < r x , then g goes to zero (Fig. 1B). Intuitively, this is because the two mechanisms are playing tug-of-war over the firing rate, each ramping up (or down) in a fruitless effort to bring the firing rate to its target. In control theory, this phenomenon is called "wind-up." In the (degenerate) case r x = r g , the nullclines overlap perfectly, forming a line of fixed points (Fig. 1C). This situation is undesirable because it leaves an extra degree In neurons with constant firing rate, dual homeostasis fails to converge on a set point. A firing rate unit receives constant input I (t) = 1. It is controlled by the homeostatic x (intrinsic homeostasis) and g (synaptic homeostasis) as described by (2) with f x (r) = r and f g (r) = r 2 . Other parameters are listed in Appendix 2. Vector fields of the control system are illustrated with arrows in the (x, g) phase plane. The x-and g-nullclines are plotted with sample trajectories in the phase plane (above), and these sample trajectories are plotted over time (below). (A) If the target firing rate r x of the excitability-modifying homeostatic mechanism is lower than the target firing rate r g of the synaptic scaling mechanism (in this case, r x = 2.5 and r g = 3.5), then g increases and x decreases without bound. This phenomenon is called "controller wind-up." (B) If r x > r g (in this case, r x = 2.5 and r g = 3.5), then g → 0, i.e., all afferent synapses are eliminated. (C) If r x = r g (in this case, r x = r g = 3), then the nullclines lie on top of each other, creating a one-parameter set of fixed points that collectively attract all control system trajectories of freedom: homeostasis has no unique fixed point, so the neuron could reach a set point with any synaptic strength, including arbitrarily strong or weak synapses, depending on initial conditions. Further, this state is destroyed by any perturbation to the target rates, so it could not be easily sustained in a biological system.
These results might lead us to believe that a control system consisting of two homeostatic control mechanisms cannot drive a neuron's firing rate toward a single stable set-point. However, we shall find that this is only because we have posed the problem in the context of a perfectly constant input I (t). When I (t) varies, the resulting picture is very different.

Varying Input
We now consider an input I (t) that is not constant at φ, but instead fluctuates randomly around φ. One simple example is I (t) = φ + σ ξ(t), where ξ(t) is white noise with unit variance. In this case, r is an Ornstein-Uhlenbeck (OU) process and is described by the stochastic differential equation An OU process approaches a stationary Gaussian distribution from any initial condition after time T τ r . In this case, this distribution has mean gφ + x and variance g 2 σ 2 2τ r .
Why does the introduction of variations in I (t) change the situation at all? This is closely connected with the basic insight that the mean value of a function f over some distribution of arguments r, written f (r) , may not be the same as the function f applied to the mean of the arguments, written f ( r ). The mean value of f (r) is affected by the spread of the distribution of r and by the curvature of f . Only linear functions f have the property that f (r) = f ( r ) for all distributions of r.
As a consequence, "satisfying" both homeostatic mechanisms may not require the condition r x = r = r g . The value ofẋ averaged over time may be zero even when the average rate r is not exactly r x , and the average value ofġ may be zero when r is not exactly r g . The conditions required to satisfy each mechanism depend on the entire distribution of r, including the mean r and the variance var(r). As long as at least one of the homeostatic mechanisms controls var(r) and r , the system has two degrees of freedom and therefore may satisfy the two fixed-point equations nondegenerately, that is, at a single isolated fixed point.
Example 1 (Rate model with linear and quadratic feedback) In order to more clearly see the influence of input variation on the control system, we explore the case in which f x (r) := r and f g (r) := r 2 . Substituting into equation (2) to describe the averaged dynamics of x and g, we have For the OU process r, we have r = gφ + x and r 2 = var(r) + r 2 = g 2 σ 2 2τ r + (gφ + x) 2 , so A vector field, nullclines, and trajectories for this system are plotted for σ 2 = 0 (constant input) in Fig. 1. The same system with σ 2 = 0.001 (varying input) is represented in Fig. 2.
The fixed points of this ODE can be studied using basic dynamical systems methods. Settingẋ =ġ = 0, we find that this equation has fixed points at (x * , g * ) = (r x , 0) . We are not interested in the first fixed point because it has g * = 0. Of the next two fixed points, we are interested only in the one with nonnegative g * . This fixed point exists with g * = 0 if and only if the term under the square root is positive, that is, if and only if r g > r x . It is asymptotically stable (i.e., attracts trajectories from all initial conditions in its neighborhood) if the Jacobian of the ODE at the fixed point has two eigenvalues with negative real part. If one or more eigenvalues have positive real part, then it is asymptotically unstable. At this fixed point, the Jacobian is J = , and it is easy In neurons with fluctuating firing rate, dual homeostasis is effective under certain conditions. A firing rate unit receives a variable synaptic input I (t). The parameters are listed in Appendix 2. Vector fields of equation (5) in the (x, g) phase plane are illustrated with arrows. The x-and g-nullclines are plotted with sample trajectories in the phase plane (above), and these sample trajectories are plotted over time (below). (A) If the target firing rate r x of the intrinsic homeostatic mechanism is lower than the target firing rate r g of the synaptic scaling mechanism (in this case, r x = 2.5 and r g = 3.5), then the nullclines cross, and all trajectories are attracted to the fixed point at their intersection. (B) If r x > r g (in this case, r x = 3.5 and r g = 2.5), then the nullclines do not cross, and g goes to zero. (C) If r x > r g and f x is exchanged with f g (r x = 3.5, r g = 2.5, f x (r) = r 2 , f g (r) = r), then the nullclines do cross, but the resulting fixed point is unstable to check that both eigenvalues have negative real part. We conclude that this (averaged) system possesses a stable homeostatic set-point if and only if r g > r x . At such a fixed point, the firing rate has mean r = r x and variance (r − r ) 2 = r 2 g − r 2 x . Conversely, given any firing rate mean μ * and variance ν * ≥ 0, we can choose targets to stabilize the neuron with this firing rate mean and variance by setting r x = μ * and r g = ν * + r 2 x . Note that this equation is not dependent on φ or σ , the parameters of the input I (t). Thus, if φ or σ changes, then the homeostatic control system will return the neuron to the same characteristic mean and variance.
In Fig. 2A, r g > r x . In this case, a stable set point exists, and it is evident that it is reached by the following process: 1. If the mean firing rate is below r x , then g and x increase, both acting to increase the mean firing rate until it is in the neighborhood of r x . If the mean firing rate is above the targets, then g and x both decrease to lower the mean firing rate to near r x . 2. If g is now small, then the firing rate variance var(r) is small, and the second moment r 2 = var(r) + r 2 is close to r 2 x . Once r slightly exceeds r x , the averaged control system haṡ x is still less than r 2 g , sȯ Alternatively, if g is large, then var(r) is large, so the second moment r 2 exceeds r 2 g , whereas r is still below r x ,ġ is negative, andẋ is positive. 3. g slowly seeks out the intermediate point between these extremes, where the variance of r makes up the difference between r x and r g . As it does so, x changes in the opposite direction to keep the mean firing rate near r x .
In Fig. 2B, it is evident that when r g < r x , no such equilibrium exists. In Fig. 2C, we show that if we exchange f x with f g and τ x with τ g such that g is linearly controlled and x is quadratically controlled, then an equilibrium exists for r g < r x , but it is not stable. These observations suggest that certain general conditions must be met for there to exist a stable equilibrium in a dually homeostatic system. We explore these conditions in the next section.
Note that if only x were dynamic but not g, the firing rate variance var(r) would be fixed at g 2 σ 2 2τ r , and therefore the variance at equilibrium would be sensitive to changes in σ , the variance of the input current. If only g were dynamic but not x, the firing rate mean and variance could both change over the course of g-homeostasis, but the two would be inseparably linked: using the expressions above for firing rate mean and variance, we can see that no matter how g changed we would always have r = φ √ 2τ r var(r) σ + x. Thus, the neuron could only maintain a characteristic firing rate mean and variance if they satisfied this constraint.

Analysis
Now we shall consider the general case in which two control variables a and b evolve according to arbitrary control functions f a and f b and control the distribution of a neuron's firing rate r. We make the dependence of this distribution on a and b explicit by writing the distribution of r as P (r; a, b). We address several questions to this model. First, what fixed points exist for a given control system, and what characterizes these fixed points? Second, under what circumstances are these fixed points stable?
In this section, we answer these questions under the simplifying assumption that f a (r) and f b (r) are constant on any domain where P (r; a, b) > 0. In Appendices 1 and 2, we show that our results persist qualitatively for nonconstant f a and f b .
In Theorem 1, we write expressions for the firing rate mean μ * and variance ν * that characterize any fixed point (a * , b * ). From this result we find that the difference between the two target firing rates plays a key role in establishing the characteristic variance at a control system fixed point.
In Theorem 2, we present a general condition that ensures that a fixed point (a * , b * ) of the averaged control system is stable. This condition takes the form of a relationship between convexity of the control functions and the derivatives of the first and second moments of P (r; a, b) with respect to a and b.

Definitions
Consider a pair of homeostatic variables a and b whose instantaneous rates of change are functions of a firing rate variable r: where 0 < 1 is a multiplier separating the fast time scale of firing rate variation and the slow time scale of homeostasis, τ a and τ b are homeostatic time constants (in units of slow time), r a and r b are the target firing rates of the two homeostatic mechanisms, and f a and f b are smooth increasing bounded functions with bounded derivatives. Note that we have introduced the small parameter representing the separation of the time scales of homeostasis and firing rate dynamics rather than assuming that τ a and τ b are large. This form is sufficiently general to encompass a wide range of different feedback strategies.

Remark 1
In order to describe the evolution of a homeostatic variable a that acts multiplicatively and must remain positive (e.g., the synaptic scaling multiplier g used in many of our examples), we can instead set τ a ȧ = a(f a (r a ) − f a (r)). We can then put this system into the general form above by replacing a withã := log(a), whose evolution is described by the ODE τ a ȧ = f a (r a ) − f a (r).
We assume that, for fixed a and b, the firing rate r(t; a, b) (written as a function of time and control variables) is distribution-ergodic with stationary firing rate distribution P (r; a, b), that is, lim T →∞ 1 T T 0 f (r(t; a, b)) dt = R f (r)P (r; a, b) dr with probability 1 for all integrable functions f . For brevity of notation, we let · (a,b) := E(·|a, b) denote the expected value of a function of r over the stationary distribution P (r; a, b) (or, equivalently, the time average of this function over time T → ∞), given a control system state (a, b). Let μ(a, b) and ν(a, b) denote the mean and variance of P (r; a, b), respectively. Averaging (6) over the invariant distribution, we arrive at the "averaged equations": We use the averaged equations to study the behavior of the unaveraged system (6). Since r may be constantly fluctuating, a and b may continue to fluctuate even once a fixed point of the averaged system has been reached, so we cannot expect stable fixed points in the classical sense. Instead, we define a weaker form of stability.
We shall call a point (a * , b * ) "stable in the small-limit" if there exists a continuous increasing function α with α(0) = 0 and an * > 0 such that for all 0 < < * , the (a, b) trajectory initialized at (a * , b * ) remains within a radius-α( ) ball centered at (a * , b * ) for all time with probability 1. Intuitively, a point is stable in the smalllimit if trajectories become trapped in a ball around that point, and the ball is smaller when homeostasis is slower. (7) is a stable fixed point of the original system (6) in the small-limit.

Lemma 1 Any exponentially stable fixed point of the averaged system
Proof This follows from Theorem 10.5 in [15].

Fixed Points
Given a homeostatic control state (a * , b * ), it is straightforward to find the target firing rates that make that state a fixed point in terms of the average values of f a and f b . By settingȧ =ḃ = 0 in (7) we find that r a = f −1 a ( f a (r) (a * ,b * ) ) and . (These expressions are well defined because f a is increasing and hence invertible, and f a (r) (a * ,b * ) must fall within the range of f a ; likewise for b.) Given a pair of target firing rates r a and r b and functions f a and f b , we can ask what states (a * , b * ) become fixed points of the averaged system. We shall answer this question in order to show that (1) when f a and f b are constant, the fixed points are exactly the points at which P (r; a, b) attains a certain characteristic mean and variance, (2) the relative convexities of the control functions determine whether r a or r b must be larger for fixed points to exist, and (3) fixed points with high firing rate variance are achieved by setting r b far from r a .

Theorem 1
Consider a dual control system as described in Sect. 3.1 with target firing rates r a and r b and control functions f a and f b .
We consider a domain of control system states (a, b) on which each distribution P (r; a, b) has constant f a (r) and f b (r) on its support. The fixed points of the averaged control system in this domain are exactly the points We will henceforth call μ * and ν * the "characteristic" mean and variance of any neuron regulated by this control system.

Remark 2
Note that this result is a/b symmetric. If a and b are swapped, then the signs of r b − r a , K b − K a , and k reverse; however, since these terms all occur in pairs, this reversal leaves the expressions for μ * and ν * unchanged.

Remark 3
In Appendix 1, we show that this result persists in some sense for nonconstant f a and f b . Specifically, if variation in f a and f b over the appropriate domain is small, the mean and variance at any fixed point are close to μ * and ν * , and every point at which the mean is μ * and the variance is ν * is close to a fixed point.
Proof We abbreviate μ(a, b) as μ and ν(a, b) as ν. Since f a and f b have constant second derivatives on the domain of interest, we can write Taking the expected values of both sides, we have At a fixed point (a * , b * ) of the averaged control system with firing rate mean μ * and variance ν * , Deriving a similar expression by expanding f b (r) around r b , we have Multiplying (8) by K b and (9) by K a and taking the difference, the ν * terms cancel, leaving Solving for μ * , we have Taking the difference of (8) and (9), we get or, substituting for μ * and solving for ν * , Given the parameters of the control system (including a pair of target firing rates), this theorem shows that achieving a specific firing rate mean and variance is necessary and sufficient for the time-averaged control system to reach a fixed point. If P (r; a, b) (the distribution of the firing rate as a function of a and b) changes, as it might as a result of changes in the statistics of neuronal input, then the new fixed points will be the new points at which this firing rate mean and variance are achieved. Conversely, given a desirable firing rate mean and variance, we could tune the parameters of the control system to make these the characteristic mean and variance of the neuron at control system fixed points.
Whether any fixed point (a * , b * ) actually exists depends on whether the characteristic firing rate mean and variance demanded by Theorem 1 can be achieved by the neuron, that is, fall within the range of μ(a * , b * ) and ν(a * , b * ). If the mapping from (a, b) to (μ, ν) is not degenerate, then there exists a nondegenerate (two-parameter) set of reachable values of μ and ν for which control system fixed points exist. In the degenerate case that neither μ nor ν depend on b, the set of reachable values of μ and ν are a degenerate one-parameter family in a two-dimensional space. This corresponds to the case of a single-mechanism control system. In this case, a control system possesses a fixed point with a given firing rate mean and variance only if they are chosen in a particular relationship to each other. A perturbation to neuronal parameters would displace this one-parameter family in the (μ, ν)-space, likely making the preperturbation firing rate mean and variance unrecoverable.
We now prove a corollary giving a simpler form of Theorem 1, which holds if r b − r a is sufficiently small.
| are sufficiently small, then the characteristic mean and variance given in Theorem 1 are arbitrarily well approximated by Proof For k defined in Theorem 1, we can write . This gives us an approximation of μ * . We can also use it to write All of these terms are bounded in norm by multiples of either K|r b − r a | or so this expression is arbitrarily small. This gives us an approximation for ν * .
The range of r b − r a for which this result holds is determined by K a and K b , measures of the convexities of the control functions. Informally, we say that this corollary holds if r b − r a is "small on the scale of the convexity of the control functions." From the corollary we draw two important conclusions that hold while r b − r a remains small on the scale of the convexity of the control functions: 1. Since a negative firing rate variance can never be achieved by the control system, there can only be a fixed point if r b − r a and K b − K a take the same sign. 2. Increasing |r b − r a | causes a proportionate increase in control system's characteristic firing rate variance.

Fixed Point Stability
Next, we address the question of whether a fixed point of the averaged control system is stable. We again make the simplifying assumption that f a and f b are constant and then drop this assumption in Appendix 2.
Theorem 2 Let (a * , b * ) denote a fixed point of the averaged control system described above. We assume the following: Let μ * = μ(a * , b * ) and ν * = ν(a * , b * ) denote the firing rate mean and variance at this fixed point. Below, all derivatives of μ and ν with respect to a and b are evaluated at (a * , b * ).
The fixed point (a * , b * ) of the averaged system is stable if Remark 4 Note that this result is a/b symmetric: if a and b are swapped, then the signs of both terms reverse and these sign changes cancel, leaving the stability condition unchanged.
Proof A fixed point (a * , b * ) of the averaged system is exponentially stable if the Jacobian has two negative eigenvalues. By Assumption 2, the Jacobian of the dual control system at (a * , b * ) has negative trace. A matrix has two negative eigenvalues if it has a negative trace and positive determinant. Therefore, the fixed point (a * , b * ) of the averaged control system is exponentially stable if Below, we abbreviate μ(a, b) as μ and ν(a, b) as ν. In order to write useful expressions for the terms in det(J ), we Taylor-expand f a (r) about μ out to second order, writing f a (r) = f a (μ) + f a (μ)(r − μ) + r μ (r − s)f a (s) ds. We similarly expand f b (r) and average these expressions at (a, b) to rewrite the averaged control equations: (a,b) . Differentiating these expressions and evaluating them at (a * , b * ), we calculate the terms in det(J ): where all derivatives are evaluated at (a * , b * ). We have assumed that f a (r) ≡ f a (μ * ) over the support of P (a * ,b * ) , so we can write Likewise for the other three terms ∂F a ∂a , ∂F b ∂a , and ∂F b ∂b . Calculating the determinant of J and canceling like terms, we have Thus, Remark 5 In Appendix 2, we drop the assumption that f a and f b are constant over the range of r and derive a sufficient condition for stability of the form ∂μ ∂a for a Δ ≥ 0 that is close to zero if f a and f b do not vary too widely over most of the range of r.
Remark 6 A similar result to Theorem 2 could be proven for a system with a single slow homeostatic feedback. The limitation of such a system would be in reachability.
As the single homeostatic variable a changed, the firing rate mean μ and variance ν could reach only a one-parameter family of values in the (μ, ν)-space. Thus, most mean/variance pairs would be unreachable. A perturbation to neuronal parameters would displace this one-parameter family in the (μ, ν)-space, likely making the original mean/variance pair unreachable. Thus, a single homeostatic feedback could only succeed in recovering its original firing rate mean and variance after perturbation in special cases.
By Lemma 1, fixed points of the averaged system that satisfy the criteria for stability under Theorem 2 are stable in the small-limit for the full, un-averaged control system.

Further Single-Cell Examples
We will focus on examples in which the two homeostatic variables are excitability x and synaptic strength g, respectively, as discussed before.
The generality of the main results above (which require that f a and f b be constant) allows us to investigate a range of different models of firing rate dynamics even if we do not have an explicit expression for the rate distribution P . We only need to know dependence of the firing rate mean μ and variance ν on the control variables to use Theorem 1 to identify control system fixed points and to use Theorem 2 to determine whether those fixed points are stable.
In the more general case addressed in Appendix 2, where the second derivatives are not assumed to be constant, the left side of (10) must be positive and sufficiently large to guarantee the existence of a stable fixed point, where the lower bound Δ for "sufficiently large" is close to zero if f a and f b are nearly constant over most of the distribution P (r; a, b). We will further discuss the simpler case, but all of our analysis can be applied to the more general case by replacing "positive" with "sufficiently large." Returning to Example 1 We can use Theorem 2 to generalize the results for the dually controlled OU process from Sect. 2. We consider the rate model described by the differential equation where I (t) is any second-order stationary process, that is, a process with stationary mean φ = I (t) and stationary autocovariance function R(w) = I (t)I (t + w) − φ 2 , both independent of t. We assume that r is ergodic. Let μ * and ν * denote the characteristic firing rate mean and variance determined from the parameters of the control system parameters using Theorem 1. This firing rate process is the output of a stationary linear filter applied to a secondorder stationary process, so according to standard results, we have μ(x, g) = gφ + x and ν(x, g) Thus, as long as ν * > 0, there exists a fixed point of the control system at At this fixed point, we have ∂μ ∂x = 1, ∂μ ∂g = φ, ∂ν ∂x = 0, and ∂ν ∂g = g * C τ r . This gives us  Figure 3 shows simulation results for this system under conditions sufficient for stability. Note that if the statistics of the input I (t) change, the fixed point changes so that the system maintains its characteristic firing rate mean and variance at equilibrium. In Fig. 4A and Fig. 4B, we alter this system by increasing and by giving I (t) correlations on long time scales, respectively. In both these cases, trajectories fluctuate widely about the fixed point of the averaged system but remain within a bounded neighborhood, consistent with the idea of stability in the small-limit.
The slopes f g and f x can be understood as measures of the strength of the homeostatic response to deviations from the target firing rate, and the second derivatives f g and f x can be understood as measures of the asymmetry of this response for (C) Mean firing rate r (above) and firing rate variance var(r) (below) are calculated as a function of homeostatic state (x, g) and represented by color in x/g parameter space. The fixed point is marked in white. (D) Mean firing rate r (above) and firing rate variance var(r) (below) are plotted over time for both initial conditions. (E)-(H) When the mean and variance of the input current are increased, the system seeks out a new homeostatic fixed point. Note in G and H that, in spite of the new input statistics, a fixed point is reached with the same firing rate mean and variance upward and downward rate deflections. If x and g are rescaled to set f g ≈ f x , then Theorem 2 predicts that dual homeostasis stabilizes a fixed point with a given characteristic mean and variance if f g (μ * ) > f x (μ * ), thats is, if the (signed) difference between the effects of positive rate deflections and negative rate deflections is greater for the synaptic mechanism than for the intrinsic mechanism.
Example 2 (Intrinsic noise) In Example 1, we assumed that all of the rate fluctuation was due to fluctuating synaptic input. If we introduce an intrinsic source of noise (e.g., channel noise), then the picture becomes slightly more complicated. We set where ξ(t) is unit-variance white noise independent of I (t), and η sets the magnitude of the noise. The same calculations as before show that the conditions for stability under Theorem 2 are met at any fixed point for  Fig. 1A, but fluctuate randomly within that neighborhood. By Lemma 1 these trajectories converge in the small-limit, so this neighborhood represents the ball of radius α( ) that traps all trajectories and shrinks to zero as → 0. (C)-(D) The system described in Fig. 1A is modified by introducing long temporal correlations into the time course of the input current: I (t) is an Ornstein-Uhlenbeck process described by the SDE τ I dI = −I dt + dξ , where ξ is white noise with unit variance and τ I = 10 s. Again, trajectories enter and remain within a large neighborhood of the fixed point observed in Fig. 1A but fluctuate randomly within that neighborhood firing rate variance includes the noise variance: ν(x, g) = g 2 C 2τ r + η 2 2τ r . Under Theorem 1, a fixed point only exists if control system parameters are chosen to establish a characteristic variance of ν * > η 2 2τ r . This neuron cannot be stabilized with variance less than η 2 2τ r because a variance that low cannot be achieved by the inherently noisy neuron.
In Fig. 5, we show the behavior of this system when ν * > η 2 2τ r (the mean and variance necessary for a fixed point are in the ranges of μ and ν) and when ν * < η 2 2τ r (the necessary variance is not in the range of ν).
Example 3 (Poisson-spiking neuron with calcium-like firing rate sensor) In some biological neurons, firing rate controls homeostasis via intracellular calcium concentration [4]. Intracellular calcium increases at each spike and decays slowly between spikes, and it activates signaling pathways that cause homeostatic changes. Our dual homeostasis framework is general enough to describe such a system. We let ρ represent the concentration of some correlate of firing rate, such as intracellular calcium, and use it in place of firing rate r. We model neuronal spiking as a Poisson process with rate λ(t) = gI (t) + x, where I (t) is a stationary synaptic input. We let ρ in- x winds up without bound, and g winds down toward zero. (D) Note that, due to intrinsic noise, the firing rate variance everywhere in parameter space is larger than the characteristic variance reached at equilibrium in B. Thus, the characteristic variance of this system at equilibrium is unreachable, and dual homeostasis does not converge crease instantaneously by δ at each spike and decay exponentially with time constant τ d between spikes. We assume that ρ is ergodic.
We show in Appendix 4 that, after sufficient time, y = ρ assumes a stationary distribution with mean μ(x, g) = δτ d (gφ + x) and variance ν(x, g) = δ 2 τ d (Cg 2 + gφ+x 2 ), where φ is the stationary mean of I (t), and C is a positive constant determined by the stationary autocovariance of I (t Note that the conditions for stability in this model are the same as the conditions in the firing rate models. In [8], we show the same result empirically for biophysically detailed model neurons. What all these models have in common is that changes in g significantly affect the firing rate variance in the same direction, whereas x controls mainly the firing rate mean and has little or no effect on the variance. These results suggest that is a general, model-independent condition for stability of synaptic/intrinsic dual homeostasis. In Appendix 2, where control function second derivatives are not assumed to be constant, this condition is replaced by the condition of sufficiently large As in Example 2, not all mean/variance pairs can be achieved by the control system: no matter how small g is, we still have ν(x, g) ≥ δ 2 τ gφ+x 2 = δμ(x,g) 2 due to the inherently noisy nature of Poisson spiking, which acts as a restriction on the range of ν. We also must have r > 0, so the range of μ is constrained to μ(x, g) > 0. If r x and r g are chosen such that the characteristic firing rate mean μ * and variance ν * defined in Theorem 1 obey these inequalities, then there exists a control system state (x * , g * ) at which Theorem 1 is satisfied and which is therefore a fixed point.

Recurrent Networks and Integration
A recurrent excitatory network has been shown to operate as an integrator when neuronal excitability and connection strength are appropriately tuned [9,10]. Such a network can maintain a range of different firing rates indefinitely by providing excitatory feedback that perfectly counteracts the natural decay of the population firing rate. When input causes the network to transition from one level of activity to another, the firing rate of the network represents the cumulative effect of this input over history. Thus, the input is "integrated." Below, we show that the parameter values that make such a network an integrator can be stably maintained through dual homeostasis as described before. Importantly, we also show that an integrator network made stable by dual homeostasis is robust to variations in control system parameters and (as in the previous examples) unaffected by changes in input mean and variance. In this section, we build intuition for this phenomenon by investigating a simple example network consisting of one self-excitatory firing rate unit, which may be taken to represent the activity of a homogeneous recurrent network. In Appendix 5, we perform similar analysis for N rate-model neurons with heterogeneous parameters. In this case, we do not prove stability, but we do demonstrate that if any neuron's characteristic variance is sufficiently high, then the network is arbitrarily close to an integrator at any fixed point of the control system.
We consider a single firing rate unit described by the equatioṅ where η is the level of intrinsic noise, I (t) is a second-order stationary synaptic input with mean φ and autocovariance R(w), and ξ(t) is a white noise process with unit variance. (For simplicity, we have rescaled time to set the time constant τ r to 1.) Let m(t) denote the expected value of r at time t. Taking the expected values of both sides of the equation, we haveṁ Let μ denote the expected value of r once it has reached a stationary distribution. Settingṁ = 0, we calculate Let s(t) denote the deviation of r from m at time t: s(t) := r(t) − m(t). From the equations above we haveṡ If we set g = 1, then the s-dependence drops out of the right side, and we have In this extreme case, s acts as a perfect integrator of its noise and its input fluctuations, that is, as a noisy integrator. For g close to 1, the mean-reversion tendency of s is weak, so on short time scales, s acts like a noisy integrator. Next, we write a differential equation for the variance of r. From (18) we write Squaring both sides out to O(dt) and taking the expected value, we have where C is a positive constant depending on τ r and R(w) as in the previous section. Let ν := lim t→∞ var(r(t)) = lim t→∞ s(t) 2 denote the expected variance of r when it has reached a stationary distribution. At this stationary distribution, we have This relation between g and ν is plotted in Fig. 6. The right side of this equation is η 2 2 at g = 0 and increases with g until it asymptotes to infinity at g = 1. So, given ν * > η 2 2 , there exists exactly one g * at which a firing rate variance of ν * is achieved. The larger the characteristic variance, the closer g * will be to unity. As discussed before, the firing rate is a good integrator on short time scales if g is close to unity. So, given a sufficiently large characteristic variance ν * , the system's only potentially stable state in the range 0 < g * < 1 will allow it to act as a good integrator on short time scales. The larger ν * , the more widely ν * can be varied while still remaining sufficiently large to make g * close to unity. So, if target firing rates are chosen to make the characteristic variance ν * sufficiently large, then an integrator achieved  (19) is plotted with η 2 = 5 and C = 1. When synaptic strength is zero, all firing rate variance is due to noise, so ν = η 2 2 . As synaptic strength increases, firing rate variance increases. As synaptic strength approaches unity, recurrent excitation acts to reinforce variations in firing rate, and variance asymptotes to ∞. If target firing rates are set such that the characteristic firing rate variance ν * is large, then the synaptic strength g * at a control system fixed point must be close to unity, making the network an integrator of its inputs in this way is robust to variation in characteristic variance ν * (and unaffected by variation in μ * ).
We can use (17) f x (μ * ) > 0 and target rates are chosen to create a sufficiently large characteristic variance ν * , then dual homeostasis of intrinsic excitability and synaptic strength stabilizes a recurrent excitatory network in a state such that the network mean firing rate acts as a near-perfect integrator of inputs shared by the population. The corollary to Theorem 1 tells us that, to first approximation, the characteristic variance is proportionate to the difference between the target rates, and a large characteristic variance is achieved by setting the homeostatic target rates far apart from each other. The integration behavior created in this way is robust to variation in the characteristic mean and variance, and therefore robust to the choice of target firing rates.
This effect can be intuitively understood by noting that as a network gets closer to being a perfect integrator (i.e., as g approaches 1), fluctuations in firing rate are reinforced by the resulting fluctuations in excitation. As a result, the tendency to revert to a mean rate grows weaker and the firing rate variance increases toward infinity. (In a perfect integrator, a perturbed firing rate never reverts to a mean, so the variance of the firing rate is effectively infinite.) Thus, the network can attain a large variance by tuning g to be sufficiently close to 1. If this large variance is the characteristic variance of the control system, then there is a fixed point of the dual control system at this value of g.
In some sense, this behavior is an artifact of the model used-perfect integration is only possible if the feedback perfectly counters the decay of the firing rate over a range of different rates, which is possible in this model because rate increases linearly with feedback and feedback increases linearly with rate. However, such a balance is also achievable with nonlinear rate/feedback relationships if they are locally linear over the relevant range of firing rates. In particular, if the firing rate is a sigmoidal function of input and the eigenvalue of firing rate dynamics near a fixed point is near zero, the upper and lower rate limits act to control runaway firing rates while the system acts as an integrator in the neighborhood of the fixed point. In [8], we show that a recurrent network of biophysically detailed neurons with sigmoidal activation curves can be robustly tuned by dual homeostasis to act as an integrator.
In Appendix 5, we show that integration behavior also occurs at set points in networks of heterogeneous dually homeostatic neurons if one or more of them have a sufficiently large characteristic variance. If only one neuron's characteristic variance is large, the afferent synapse strength to that neuron grows until that neuron gives itself enough feedback to act as a single-neuron integrator as described before. But if many characteristic variances are large, then all synapse strengths remain biophysically reasonable, and many neurons participate in integration, as might be expected in a true biological integrator network.
In Fig. 7, we show simulation results for homogeneous and heterogeneous recurrent networks with target firing rates set to create sufficiently large characteristic firing rate variances. In addition to corroborating our analytical results, these simulations provide empirical evidence that the fixed points of heterogeneous networks are stable under similar conditions to those guaranteeing the stability of the single self-excitatory rate unit discussed before.

Discussion
This mathematical work is motivated by the observation that the mean firing rates of neurons are restored after chronic changes in input statistics and that this firing rate regulation is mediated by multiple slow biophysical changes [3,5]. We explore the possibility that these changes represent the action of multiple independent slow negative feedback ("homeostatic") mechanisms, each with its own constant "target firing rate" at which it reaches equilibrium. Specifically, we focus on a model in which the firing of an unspecified model neuron is regulated by two slow homeostatic feedbacks, which may correspond to afferent synapse strength and intrinsic excitability or any two other neuronal parameters.
In a previous work [8], we showed in numerical simulations that a pair of homeostatic mechanisms regulating a single biophysically detailed neuron can stably maintain a characteristic firing rate mean and variance for that neuron. Here, we have analytically derived mathematical conditions sufficient for any model neuron exhibiting Fig. 7 Dual homeostasis creates integrators from a single recurrently excitatory neuron and a heterogeneous excitatory network. (A) Dual homeostasis tunes a single neuron with a recurrent excitatory connection to function as an integrator from two different initial conditions (orange and blue). First row: the target rates r x and r g are plotted as a point in (r x , r g ) space. Second row: the resulting x and g trajectories are plotted in phase space and over time. Third row: before dual homeostasis, the neuron is tested for integrator-like behavior by injecting pulsatile input I (t). The firing rate r returns to a baseline after each pulse. Fourth row: after dual homeostasis, firing rate r increases at each pulse and retains its approximate value from one pulse to the next. This neuron is an integrator: its firing rate at any time represents an integral of the pulse history. (B) Analogous plots for a heterogeneous recurrently excitatory network of 200 neurons initialized from two initial conditions (orange and blue). Top row: target firing rate pairs of all neurons are plotted in (r x , r g ) space. Note that r g is always chosen to be greater than r x . Second row: average values of g and x across the network are plotted in phase space and over time. In phase space, a representative trajectory of a single neuron in the network for each of the two simulation runs is also plotted in lighter colors. Note that these trajectories are significantly removed from the average trajectories. Third row: both times the network is initialized, the average network rate r does not act as an integrator for pulsatile input. Fourth row: after dual homeostasis, the average network rate acts as an integrator such dual homeostasis to exhibit this behavior. Importantly, the homeostatic system reaches a fixed point when the firing rate mean and variance reach characteristic values determined by homeostasis parameters, so the mean and variance at equilibrium are independent of the details of the neuron model, including stimulus statistics. Thus, this effect can restore a characteristic firing rate mean after major changes in the neuron's stimulus statistics, as has been observed in vivo, while at the same time restoring a characteristic firing rate variance.
In Theorem 1, we have provided expressions for the characteristic firing rate mean and variance established by a specific set of homeostatic parameters. They show that when the separation between the target rates r a and r b is appropriately small, the relative convexities of the functions f a and f b (by which the firing rate exerts its influence on the homeostatic variables) determine which target rate must be larger for a fixed point to exist. When a fixed point does exist, the characteristic firing rate variance at the fixed point is proportional to the difference between r b and r a .
In Theorem 2, we find that any fixed point of our dual homeostatic control system is stable if a specific expression is positive. This expression reflects the mutual influences of firing rate on the homeostatic control system and of control system on the firing rate mean and variance.
Both these theorems are proven under the simplifying assumption that f a and f b are constant. However, in Appendices 1 and 2, we drop this assumption and find that qualitatively similar results hold as long as these second derivatives do not vary too widely. In particular, stability is guaranteed if the expression in Theorem 2 exceeds a certain positive bound that is close to zero if f a and f b are nearly constant across most of the range of variation of the firing rate.
We have explored the implications of our results for a system with slow homeostatic regulation of intrinsic neuronal "excitability" x (an additive horizontal shift in the firing rate curve) and afferent synapse strength g. From the corollary to Theorem 1 we find that (to first approximation) stable firing rate regulation requires r g > r x . Using Theorem 2, we show that for rate-based neuron models and Poisson-spiking models, stable firing rate regulation is achieved when the f g is sufficiently concaveup relative to f x .
We predict that these conditions on relative concavity and relative target firing rates should be met by any neuron with independent intrinsic and synaptic mechanisms as its primary sources of homeostatic regulation. Experimental verification of these conditions would suggest that our analysis accurately describes the interaction of a pair of homeostatic mechanisms; experimental contradiction of these conditions would suggest that the control process regulating the neuron could not be accurately described by two independent homeostatic mechanisms providing simple negative feedback to firing rate.
Our results have special implications for neuronal integrators that maintain a steady firing rate through recurrent excitation. We have found that the precise additive and multiplicative tuning necessary to maintain the delicate balance of excitatory feedback can be performed by a dual control system within the framework we study here if the target firing rates are set far enough apart to achieve a large characteristic firing rate variance. The integrator maintained by dual homeostasis is robust to variation of its target firing rates. This occurs because the circuit is best able to achieve a large firing rate variance when it is tuned to eliminate any bias toward a particular firing rate, exactly the condition necessary for integration. This robust integrator tuning scheme should be considered in the ongoing experimental effort to understand processes of integration in the brain.

Declarations
Code mentioned in the manuscript will be made available on Figshare.

Appendix 2: Generalization of Theorem 2
Here we present the sense in which the results of Theorem 2 persist for functions f a and f b with nonconstant second derivatives.
The Jacobian of the averaged control system is continuous with respect to perturbations to f a and f b in C 2 . Therefore, if the determinant of the Jacobian is negative (i.e., a fixed point of the averaged system is stable) for functions f a and f b with constant second derivatives, then when these functions are perturbed in C 2 , the Jacobian is negative at the perturbed fixed point (see the previous appendix).
The persistence of fixed point stability for nonconstant f a and f b can be expressed as a generalization of Theorem 2. Rather than requiring an expression to be positive, this generalized theorem requires that the same expression exceed some nonnegative lower bound, which is close to zero if f a and f b are sufficiently close to constants over most of the support of the distribution of firing rates r.
Theorem 3 Let (a * , b * ) denote a fixed point of the averaged control system described before. We assume the following: 1. The functions μ and ν are differentiable at (a * , b * ). 2. ∂F a ∂a and ∂F b ∂b are negative at (a * , b * ), that is, on average, a and b provide negative feedback to r near (a * , b * ).
Let μ * = μ(a * , b * ) and ν * = ν(a * , b * ) denote the firing rate mean and variance at this fixed point. Below, all derivatives with respect to a and b are evaluated at (a * , b * ).
We define The fixed point (a * , b * ) of the averaged system is stable if Remark 7 R D a is the maximal deviation of f a (r) from f a (μ * ) on the ball of radius D. ν D (a, b) is the amount of variance attributable to probability mass on that ball. Note that ν(a, b) μ(a, b) as μ and ν(a, b)  , where c ∈ [μ, r] depends on r and μ. Integrating, we have We let We now define a simpler matrixJ that approximates J : where X ab := − ∂μ ∂b f a (μ * ) − 1 2 f a (μ * ) ∂ν ∂b , and likewise for the other three terms inJ . This matrix is identical to the matrix J under the assumption of constant f a and f b . Therefore, from (13) we have We will estimate how closely the determinant ofJ approximates that of J . From (21) we write out the expected value as an integral over a probability distribution: Given any radius D > 0, we can split the integral into two parts, one on a ball of radius D about μ * and one outside that ball: Applying the intermediate value theorem to both integrals, we have for some c 1 ∈ B(μ * , D) and some c 2 ∈ R. By the definition of ν D , We compare this quantity to X ab : This bound holds for all D > 0. Taking an infimum over all D, we have Using this bound (and the equivalent bound for the other three terms of J −J ), we set bounds on how closely the determinant ofJ approximates the determinant  Let ρ be a positive quantity initialized at zero at t = 0 that increases by δ at each event of a Poisson process with rate gI (t) + x and decays exponentially with time constant τ d between events. We can write and N(ω, ds) is a Poisson random measure with mean measure λ(ds) = (gI (s) + x) ds [16]. Note that this process is doubly stochastic: I (t) is a random process, and N(ω, ds) is a random measure depending on a specific instantiation of I (t). We will write · to signify the expected value over instantiations of I (t), and · I to signify the expected value given a specific instance of I (t). According to [17], the nth moment of this process can be written The second moment is given by

Appendix 5: Recurrent Networks with Heterogeneity
Here we show that some of our results on homogeneous recurrent excitatory networks generalize to heterogeneous networks in which each neuron has its own dual control system, target rates, and level of intrinsic noise. In particular, we show that if the neurons in such a network have target firing rates that set at least one of the characteristic variances sufficiently high, then when the combined control system reaches a fixed point, the network behaves like an integrator on short time scales. This condition is robust to variation in individual neuronal parameters, including target firing rates. Since the combined control systems of a heterogeneous network with N neurons is 2N -dimensional rather than two-dimensional, we do not attempt to prove the stability of fixed points here. However, we have verified the stability of the behavior described below in simulations (data not shown, code available online).
We consider a collection of N neurons with firing rates r [1] , . . . , r [N ] connected allto-all by excitatory synapses. Each neuron n has its own homeostatic variables x [n] and g [n] , its own extrinsic second-order stationary synaptic input I (t) [n] with mean φ [n] and autocovariance R [n] (w), its own intrinsic noise η [n] ξ [n] (t) (where ξ [n] (t) is a white noise process with unit variance), a second-order stationary synaptic input I all (t) with mean φ all and autocovariance R all (w) shared by all neurons, and excitatory synaptic input from all other neurons, which together determine the firing rate as in the rate models above. Each neuron n has its own target firing rates r x[n] and r g [n] . Each neuron n delivers an additional excitatory synaptic input 1 N r [n] to each other neuron.
We write the rate model for neuron n in the form Let g := 1 N n g [n] . For g < 1, s(t) is a mean-zero OU process. If we set g = 1, then the term −(1 − 1 N n g [n] )s(t) drops out, and we can solve for s(t): s(T ) = s(0) + In this case, the mean network firing rate acts as a perfect integrator of all of the network's input fluctuations and noise. The closer g is to 1, the more s(t) acts like an integrator on short time scales. For g < 1, the vector of rates r [n] (t) is a linearly filtered second-order stationary random process and therefore approaches a stationary mean and variance. Using (25), we write a recursive equation for the covariance of s [n] (t) and s [ where g is the vector of all g [m] . Suppose that target firing rates are set such that each neuron n has characteristic variance ν * [n] . Then at any fixed point (x, g) of the combined control system, we have ν * [n] = F [n] (g) for all n. It is easy to check that F [n] (g) is 1 2 η 2 [n] at g [n] = 0, that it increases with any with respect to any g [m] on the domain { m g [m] < N} ∩ {g [n] > 0}, and that it asymptotes to infinity as g = 1 N m g [m] approaches 1. Thus, for sufficiently large ν * [n] , the constraint ν * [n] = F [n] (g) can be solved on the domain, and any solution must be close to the hypersurface g = 1. As discussed before, the system acts as an integrator on short time scales when g is close to 1; thus, if any one or more of the characteristic variances ν * [n] are sufficiently large, then at a fixed point of the combined control system, the mean network firing rate acts as an integrator on short time scales.
We can intuitively see that if only one neuron, neuron n, has a high characteristic variance, then it acts like the single-neuron integrator as g [n] approaches N , and neuron n dominates the global feedback signal. This scenario is not particularly biophysically interesting. However, if many neurons have high characteristic variance, then no individual neuron will dominate the feedback signal. In the special case that all neuron parameters are identical, a symmetry argument shows that g [n] will be close to 1 for all n-no single synapse is strong enough to make a single neuron an integrator, but the fluctuations in network firing rate cause correlated fluctuations in the individual firing rates, which collectively reinforce the network rate fluctuations, causing the network as a whole act as an integrator.
By the corollary of Theorem 1, the characteristic variance of a neuron is proportionate to the difference between target firing rates as long as this difference remains small on the scale of the convexity of f a and f b ; so if f a and f b are small, then a large characteristic variance can be achieved by setting target firing rates sufficiently far apart, and perturbations to this difference cause proportionate perturbations in characteristic variance. If the characteristic variances ν * [n] are sufficiently large to make the network an integrator, then they will remain sufficiently large even in the face of such perturbations. We conclude that a network integrator stabilized by dual homeostasis is robust to perturbations in the target firing rates r x and r g .