extfieldsofstudy: Biology Medicine

text:
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 .


s2fieldsofstudy: Biology Medicine

text:
Metabolic heterogeneity signature of primary treatment-naïve prostate cancer

To avoid over- or under-treatment of primary prostate tumours, there is a critical need for molecular signatures to discriminate indolent from aggressive, lethal disease. Reprogrammed energy metabolism is an important hallmark of cancer, and abnormal metabolic characteristics of cancers have been implicated as potential diagnostic/prognostic signatures. While genomic and transcriptomic heterogeneity of prostate cancer is well documented and associated with tumour progression, less is known about metabolic heterogeneity of the disease. Using a panel of high fidelity patient-derived xenograft (PDX) models derived from hormone-naïve prostate cancer, we demonstrated heterogeneity of expression of genes involved in cellular energetics and macromolecular biosynthesis. Such heterogeneity was also observed in clinical, treatment-naïve prostate cancers by analyzing the transcriptome sequencing data. Importantly, a metabolic gene signature of increased one-carbon metabolism or decreased proline degradation was identified to be associated with significantly decreased biochemical disease-free patient survival. These results suggest that metabolic heterogeneity of hormone-naïve prostate cancer is of biological and clinical importance and motivate further studies to determine the heterogeneity in metabolic flux in the disease that may lead to identification of new signatures for tumour/patient stratification and the development of new strategies and targets for therapy of prostate cancer.

INTRODUCTION
Prostate cancer is a clinically heterogeneous disease covering a wide variety of phenotypes, ranging from indolent disease to aggressive lethal disease. At present, there is a lack of methods for reliably identifying such phenotypes, rendering use of first-line androgen ablation therapy a challenge with a risk of over-or under-treatment [1]. As such, there is an urgent need for signatures to stratify indolent from aggressive prostate cancer which could be developed through identification of differences in energy metabolisms of prostate cancer subtypes [2].
Reprogramming energy metabolism is an important hallmark of cancer [3,4] and abnormal metabolic characteristics of cancers have been implicated as potential diagnostic/prognostic biomarkers [5][6][7][8][9]. Thus in contrast to normal resting cells whose energy production is mainly based on mitochondrial oxidative

Research Paper
Oncotarget 25929 www.impactjournals.com/oncotarget phosphorylation, the majority of cancer cells rely on aerobic glycolysis coupled to increased glucose uptake and lactic acid production/secretion (the Warburg effect) [3,7]. The augmented glucose consumption can be observed using ( 18 F)-fluoro-2-deoxyglucose positron emission tomography (FDG-PET) [10,11]. However, prostate cancers are characterized by a lack of increased glycolysis rendering FDG-PET less effective in imaging this malignancy [12], (although a small group of tumours do utilize aerobic glycolysis [9]). Instead, prostate cancers have abnormal metabolisms for citrate and choline [13][14][15] and are characterized by increased levels of cholinecontaining compounds (ChoCCs). These compounds have been linked to increased cell proliferation and survival [16]. Tracing metabolites with ( 11 C)-and ( 18 F)labeled choline derivative-based PET and PET/CT scans is increasingly used in the clinic. Together, these glucose and choline-based imaging studies demonstrated the metabolic heterogeneity of clinical prostate cancer [17,18]. However, it is not clear whether prostate cancer subtypes are metabolically different, or whether metabolic heterogeneity (intra-tumoural and inter-tumoural) can affect tumour response to treatment. In view of this, a better understanding of metabolic heterogeneity of prostate cancer at a molecular level will provide mechanistic insights into disease progression and may lead to identification of novel diagnostic/prognostic signatures for stratification of patients' prostate cancers.
To focus on the metabolic heterogeneity of prostate cancer cells and avoid complication from the stromal compartment, we used our unique panel of transplantable, patient-derived xenograft (PDX) models, derived from hormone-naïve prostate cancer tissues (www. livingtumorlab.com) [19,20]. Using these hormonenaïve prostate cancer PDX models and clinical samples, we, for the first time, demonstrated the heterogeneity of expression of genes involved in cellular energetics and macromolecule biosynthesis and identified an inter-tumour metabolic heterogeneity signature in clinical samples that is associated with patients' prognosis.

A panel of PDX models of hormone-naïve prostate cancers demonstrates heterogeneity of expression of genes involved in energy metabolism and recapitulates metabolic heterogeneity of clinical samples
A panel of 11 well-characterized LTL (Living Tumour Laboratory) prostate cancer PDX-models that are derived from hormone-naïve prostate cancer tissues, presents a unique opportunity for focusing on metabolic alterations in prostate cancer cells by analyzing human genes (while stromal genes are of mouse origin) in vivo [19]. To investigate whether these hormone-naïve PDX models present with molecular metabolic heterogeneity, transcriptome microarray data were generated and cluster analysis was performed using mean pathway scores (derived from per gene z-scores normalized to cancer mean and standard deviation, Supplementary Table 2). As shown in Figure 2, the metabolic pathway scores demonstrated variation among these models. For example, four tumours had increased pathway scores for glycolysis, while seven other tumours had decreased scores for this pathway. Heterogeneity was evident for all pathways.
To investigate whether such heterogeneity of metabolic gene expression observed in the PDX models represents the clinical situation, we analyzed transcriptome sequencing data of treatment-naïve, primary prostate cancer specimens of the VPC cohort (n = 14). The VPC prostate cancer per-gene z-scores were normalized to cancer mean and standard deviation, to generate mean pathway scores. The cluster analysis combining PDX models and VPC clinical samples showed that (in Figure  3), all the PDX models were clustered into different arms as observed for the clinical samples, suggesting that this set of PDX models recapitulates the heterogeneous metabolic nature of the clinical samples. In addition, the matched PDX/parental tissue pairs LTL-331/VPCT24 and LTL418/VPCT23, clustered closely together, as seen in dendrogram arms A and B respectively, illustrating the fidelity of the PDX models (arrow brackets). Interestingly, the LTL-313 xenograft lines, which had been derived from needle biopsies at five different foci of one patient's primary prostate cancer (asterisk), show a diverse heterogeneity signature in each focus. The observed differences between the five cores from which these models originate suggest the existence of diverging events that are continually associated with the development of intra-tumoural metabolic heterogeneity of prostate cancer. In summary, we documented molecular evidence of cancer metabolic heterogeneity in PDX models and clinical hormone-naïve prostate cancer samples. Importantly, the metabolic heterogeneity discovered in PDX models was validated in clinical prostate cancer samples.

Change in expression of genes involved in metabolism is consistent with the metabolomic changes in prostate cancer
To determine whether the gene expression involved in metabolic pathways is consistent with the metabolic properties of prostate cancer cells, we first analyzed the transcriptomic differences, between benign and treatmentnaïve prostate tumours, of genes involved in metabolic pathways and we then explored the correlation with clinical metabolomics. A comparison of the MSKCC dataset [21] consisting of 28 benign prostate tissues and Oncotarget 25930 www.impactjournals.com/oncotarget 112 prostate adenocarcinomas from treatment-naïve patients, showed significantly increased expression of genes involved in choline/phosphocholine generation (p < 0.0002) and in proline synthesis (p < 0.0029) ( Table  2) in tumours compared to benign tissue. These results reflect the findings of several recent clinical metabolomics studies reporting increased choline metabolites and increased proline [22][23][24][25] in prostate cancer compared to benign tissue.
Interestingly, expression of genes involved in pyruvate to acetyl-CoA formation was significantly decreased (p < 0.0092). As expected, there were no significant differences in the expression of glycolysis pathway genes (i.e. upstream of the PDH complex) between benign tissue and cancer (Table 2), which is consistent with previous reports and reflects the clinical finding that FDG-PET is less effective in imaging prostate malignancy [15] than other cancers.

Gene expression of a metabolic signature demonstrates heterogeneity in clinical samples and clinical outcome of prostate cancer patients
To further validate our finding of metabolic heterogeneity of treatment-naïve prostate tumours and explore its clinical significance, we analyzed the correlation of our heterogeneous metabolic signature  Z-scores were normalized to all tumours (mean/standard deviation) for each gene within a pathway. Average z-scores were then used to generate pathway scores. Red, upregulated; blue, downregulated; white, no change from tumour average. A, B, and b1 represent various arms of the dendogram.
Oncotarget 25932 www.impactjournals.com/oncotarget with clinical outcome in the MSKCC dataset [21]. As shown in Figure 4, heterogeneity of metabolic pathways between tumours is evident. Dendrogram arms I and II indicate two broadly different groups of tumours. Group II shows further heterogeneity as indicated by dendrogram arms (i and ii). The MSKCC validation cohort showed significant heterogeneity (Wilcoxon Rank test, p < 0.05) of metabolic pathways involved in cholesterogenesis, choline metabolism, gluconeogenesis, glutaminolysis, glycolysis, ketogenesis, lipogenesis, oxidative phosphorylation involving complexes I, IV and V, the TCA cycle and onecarbon metabolism (Supplementary Table 4). Pathways with less heterogeneity also occurred, including fatty acid activation and ketolysis. Interestingly, proline synthesis and proline degradation pathways showed opposite directions of expression. The proline synthesis pathway was consistently upregulated and proline degradation pathway was downregulated compared to benign tissue controls, which together supports metabolomics results [17] of increased proline content in aggressive prostate cancer. This is also consistent with results in the previous section (Table 2) where we showed that the proline synthesis pathway score was significantly increased in cancer compared to benign tissue (p < 0.0029).
To determine the clinical relevance of the observed metabolic heterogeneity we carried out survival curve analyses (5 year-overall survival as well as biochemical recurrence-free survival) for the MSKCC cohort (Supplementary Table 5). Patients with elevated expression levels of genes involved in one- The prostate cancer cohort of the Vancouver Prostate Centre (VPC) consisted of 14 untreated primary tumours and 3 matched benign prostate tissues [44]. Abbreviations: T, tumour; B, Benign; Dx, diagnosis; PSA, prostate specific antigen; nrn, never reached nadir; N, no; Y, yes; LTFU, lost to follow-up. There were no differences in stromal content across the cohort based on comparisons of desmin and vimentin expression (fibroblasts) [44]. www.impactjournals.com/oncotarget carbon metabolism, compared to controls ( > 2 Standard Deviations, Supplementary Table 6), showed significantly poorer biochemical-free survival compared to patients with lower expression of these genes (p < 0.01) ( Figure  5A-5B). Additionally, patients with low expression levels of genes involved in proline degradation, compared to controls (Supplementary Table 6), had significantly poorer biochemical-free survival compared to patients with higher expression levels ( Figure 5C-5D). Therefore, increased one-carbon metabolism or decreased proline degradation in these patients was associated with significantly decreased survival, suggesting these metabolic changes form clinically important metabolic signatures. There was a likely-correlation (Spearman) of one carbon metabolism and proline degradation pathway scores with Pathology stage (p = 0.06 and p = 0.08 respectively). This suggests that increased one carbon metabolism or decreased proline degradation (therefore increased proline levels) is observed in more aggressive tumours compared to less aggressive tumours. Such increased one carbon metabolism or decreased proline degradation signature was also observed in 50% or 75% (respectively) of PDX models with metastatic ability suggesting a predictive value of aggressiveness.

DISCUSSION
Whereas the genomic and transcriptomic heterogeneity of prostate cancer is well documented [26], less is known about metabolic heterogeneity of this malignancy and its relevance to malignant progression. Studies of the heterogeneity of prostate cancer, in particular hormone-naïve cancers, have been impeded by a lack of clinically relevant experimental models. Commonly used prostate cancer xenograft models based on cultured human prostate cancer cells do not accurately recapitulate the tumour heterogeneity, nor the cancer-stroma architecture of the original cancer specimens from which the cell lines were derived [19,20]. Importantly, the majority of prostate cancer cell lines and xenograft models used were derived from advanced metastatic castrationresistant cancer, rather than untreated hormone-naïve cancer. In this study, we demonstrated the heterogeneous expression of cancer cell genes involved in energy metabolisms in 11 patient-derived prostate cancer PDX models derived from hormone-naïve tumours. These models belong to a previously produced panel of prostate cancer PDX models [19,20] and to our knowledge currently form the largest collection of hormone-naïve prostate cancer PDX models in the field. Furthermore, the metabolic heterogeneity signature in these models suggests, for the first time, that a panel of PDX models can recapitulate heterogeneous metabolic signatures observed in clinical samples. The metabolic program that correlated with poor prognosis in the MSKCC dataset was associated with metastasis in PDX models. As such, these models provide a valuable platform for studying metabolic heterogeneity.
Using RNA sequencing strategies and microarrays to investigate heterogeneity of expression of genes involved in cellular energetics and macromolecular biosynthesis pathways, the present study documents intertumour metabolic expression heterogeneity of several pathways/sub-pathways in PDX models and in two clinical cohorts of, treatment-naïve prostate cancers as summarized in Figure 6. Consistent with these findings, heterogeneity of several metabolic pathways has been Figure 6: Summary of heterogeneity of expression of genes related to core metabolic pathways involved in cellular energetics and macromolecule biosynthesis in hormone-naive prostate cancer. The one-carbon metabolism pathway score was significantly increased (red) and the proline degradation pathway score was decreased (blue) compared to benign and both were associated with decreased biochemical-free survival in the MSKCC subset. It has been reported that the one-carbon metabolism and proline metabolism pathways interconnect with glycolysis and pentose phosphate pathway to provide ATP, and NADH and NADPH (mediators of electron transfer for energy, reductive biosynthesis and redox defense) [8,9,42] Oncotarget 25938 www.impactjournals.com/oncotarget observed in a meta-analysis of datasets of 22 different tumour types [27]. Importantly, we identified metabolic gene signatures of increased one-carbon metabolism or decreased proline degradation associated with significantly decreased biochemical disease-free survival of patients. One limitation of this study is the lack of metabolomics validation. Further integrated analysis of transcriptomic and metabolomics data and functional studies will provide more insights regarding the significance of metabolic heterogeneity in prostate cancer.
It is widely reported that cancer cells can direct glucose metabolism toward aerobic glycolysis [28]. However, consistent with our own results, emerging evidence suggests that prostate cancer cells can also direct glucose towards another ATP-generating pathway involving 3-phosphoglycerate conversion to serine through phosphoglycerate dehydrogenase (PHGDH), resulting in activation of one-carbon metabolism [29][30][31]. One-carbon metabolism facilitates NADPH generation [32] and the anabolic synthesis of amino acids, proteins, nucleotides and phospholipids, and has a role in the maintenance of redox balance and methylation reactions involved in post-translational modifications [33,34] and so may be a potential driver of oncogenesis. Indeed, methotrexate, an inhibitor of the one-carbon cycle, is used as an anticancer agent and preclinical studies are evaluating several other enzymes in this pathway [33]. Importantly, however, methotrexate has recently been associated with accelerated progression of a previously indolent prostate cancer [35]. An alternative therapeutic, metformin, has also been shown to indirectly inhibit one-carbon metabolism [36] and is being tested in preclinical trials for prostate cancer [37][38][39]. Several agents targeting metabolic enzymes exist for other metabolic syndromes that could be re-purposed for prostate cancer [37].
Proline metabolism promotes cancer cell proliferation and energy production [40]. The cycling of proline catabolism and proline synthesis acts as a redox shuttle between the mitochondria and cytosol to maintain redox homeostasis [41,42]. Furthermore, increased proline biosynthesis contributes to cancer cell growth by providing NAD+/NADP+ molecules to the pentose phosphate pathway and to glycolysis [40]. Here we show that a decreased pathway score for proline degradation is associated with significantly decreased biochemical disease-free survival of prostate cancer patients. Concomitantly, the expression of the proline synthesis pathway genes was consistently increased in prostate tumours implying a strong propensity of these tumours to activate the proline synthesis cycle. In support of our findings, LC/GC-MS metabolomic profiling identified increased proline levels with increased progression of prostate cancer [22][23][24][25].
Increased proline synthesis has also been reported as the major metabolic shift in metastatic breast cancer cells, in melanoma cell lines compared to melanocytes and in ovarian cancer cells (OVCAR3) compared to ovarian cancer stem cells [40]. Furthermore, a recent study of 133 metabolic genes identified five genes whose siRNAinduced knockdown resulted in strong growth-inhibitory effects on breast cancer cells. Included were PYCR1 involved in proline synthesis and the serine biosynthesis gene, PDGDH, involved in one-carbon metabolism [43].
In conclusion, we have, for the first time, demonstrated heterogeneity of expression of genes related to metabolic pathways in both hormone-naïve prostate cancer PDX models and clinical samples and provided insight into the metabolic heterogeneity of hormone-naïve prostate cancer associated with patient clinical outcome. These results not only suggest the potential of a metabolic heterogeneity signature as a biomarker for patient risk stratification, but also motivate further work to determine whether the observed differences in expression of genes involved in metabolic pathways translate into differences in metabolic flux. A better understanding of metabolic heterogeneity of hormone-naïve prostate cancer will lead to the identification of new signatures for patient stratification and the development of new strategies and targets for therapy of prostate cancer.

Patient-derived xenograft (PDX) modelsdiscovery cohort
Transplantable xenograft lines of patients' prostate cancer tissues were developed by grafting of fresh cancer tissue (obtained from needle biopsies), and serial transplantation, into the subrenal capsule (SRC) graft site of male NOD/SCID mice supplemented with testosterone as previously described [19]. Animal care and experiments were carried out in accordance with the guidelines of the Canadian Council on Animal Care.

Gene expression analysis of xenografts using RNA microarray analysis
RNA microarray analysis of PDX tissue was performed as previously described [19]. Data are available at GEO accession number GSE41193. Total RNA samples were prepared and processed using Agilent's One-Color Microarray-Based Gene Expression Analysis Low Input Quick Amp Labelling v6.0. An input of 100 ng of total RNA was used to generate cyanine-3-labeled cRNA. Samples were hybridized on Agilent SurePrint G3 Human GE 8X60K Microarray (Design ID, 028004). Arrays were scanned with the Agilent DNA Microarray Scanner at 3-μm scan resolution and data were processed with Agilent Feature Extraction 11.0.1.1. The processed signal was quantile normalized with Agilent GeneSpring 12.0. www.impactjournals.com/oncotarget

Validation cohorts
The prostate cancer cohort of the Vancouver Prostate Centre (VPC) consisted of 14 untreated primary tumours and 3 matched benign prostate tissues [44] as described in Table 1. Gleason scores ranged from 7-10 and samples consisted of ≥ 70% of malignant glands. The MSKCC cohort consisted of 112 untreated primary prostate tumours and 28 benign prostate tissues [21]. Data are available at GEO accession number GSE21032.

Gene expression analysis of tumor sections using RNA sequencing
VPC clinical tumour sections were processed and RNA sequencing was performed at the BCCA Michael Smith Genome Sciences Centre according to standard protocol as previously described [44]. Briefly, RNAseq data reads were first mapped onto the hg19 human reference genome and exon-exon junctions by spliceaware aligner STAR [45], using known gene model annotation from Ensembl Release 75. Reads with an unmapped mate or multi-mapped location were filtered out using Bam Tools [46]. Data are available at GEO accession number GSE55016.
Using aligned RNA-seq reads, gene expression profiles for each sample were calculated on the basis of gene annotation from Ensembl Release 75. Only reads that were unique to one gene and exactly corresponded with the structure of the gene were counted for the corresponding genes by using tool HTSeq [47]. In order to eliminate the variance of sequencing depth among samples, the raw read counts were normalized by R package DESeq [48].

Metabolic pathway score
Twenty-three pathways analyzed in this study include those of choline metabolism [13,14] as well as energy production and macromolecule biosynthesis, as previously curated for cancer cells [6,27,49,50] outlined in Figure 1 (and detailed in Supplementary Table 1). Assigning an 'upregulation' or 'downregulation' status to a metabolic pathway requires the contribution of all genes in the pathway. To this end, we selected well-established genes for each pathway, excluding upstream effectors.
In the PDX models, gene expression z-scores were derived from the Mean/SD of all 11 xenograft lines. The pathway score was then calculated from the average z-scores of genes within each pathway/sub-pathway. To compare the VPC cohort with the PDX models, z-scores for VPC gene expression were derived from the Mean/SD of the 14 VPC cancer samples. We obtained enough matched patient tissue for a paired comparison for 2 of the PDX models and original patient tissues (VPCT24:LTL331 and VPCT23:LTL418).
To validate metabolic heterogeneity in a larger cohort of primary untreated prostate cancer specimens (MSKCC), and to ascertain potential clinical predictive value, prostate cancer gene expression was normalized to benign gene expression (Mean/SD) to generate per gene z-scores for the genes expressed in 23 metabolic pathways/sub-pathways ( Figure 1).

Statistical analyses
The pathway scores calculated from the average z-scores of genes within a pathway provide information about the direction of the gene/pathway regulation. In the MSKCC cohort, Wilcoxon rank tests (p < 0.05) were applied to determine if pathway scores were significantly up-or down-regulated in tumours compared to benign tissue. Pathways containing fewer than 6 genes were not tested by the Wilcoxon rank test. The MSKCC cohort was also subjected to t-test (p < 0.05) using log2 gene expression (cancer vs benign) and survival curve analysis using pathway z-scores normalized to benign (Kaplan Meier plots) performed in Graphpad 6 Prism (significance determined using the Mantel-Cox test and the Gehan-Breslow test). Cox regression analysis and Spearman correlation were also undertaken to determine if metabolic pathway scores correlated with clinical parameters. Hierarchical clustering was used as previously described [51] using the R language (http://cran.r-project.org/).

CONFLICTS OF INTEREST
The authors declare that they have no conflicting interests.

GRANT SUPPORT
This study was supported in part by the Canadian

Author contributions
DL and SLE were involved in the study design, data analysis, data interpretation, literature search and generation of figures; SQ carried out data analysis; SQ, NN were involved in generation of figures; NN, SYCC and PWG performed literature searches; HX, AMH, SYCC carried out data collection; FM and RHB performed