Divisive Gain Modulation with Dynamic Stimuli in Integrate-and-Fire Neurons

The modulation of the sensitivity, or gain, of neural responses to input is an important component of neural computation. It has been shown that divisive gain modulation of neural responses can result from a stochastic shunting from balanced (mixed excitation and inhibition) background activity. This gain control scheme was developed and explored with static inputs, where the membrane and spike train statistics were stationary in time. However, input statistics, such as the firing rates of pre-synaptic neurons, are often dynamic, varying on timescales comparable to typical membrane time constants. Using a population density approach for integrate-and-fire neurons with dynamic and temporally rich inputs, we find that the same fluctuation-induced divisive gain modulation is operative for dynamic inputs driving nonequilibrium responses. Moreover, the degree of divisive scaling of the dynamic response is quantitatively the same as the steady-state responses—thus, gain modulation via balanced conductance fluctuations generalizes in a straight-forward way to a dynamic setting.


Introduction
Gain modulation (or gain control) is an adjustment of the inputoutput response of neurons, and is widely observed during neural processing [1]. Gaze direction sets the response gain in primary visual [2], posterior parietal cortex [3], and auditory brainstem [4]. In specific species, gain control mechanisms produce an invariance of receptive field properties [5] and orientation selectivity [6] to changes in overall stimulus contrast. Higher cognitive processes, such as attention, modulate the response gain of cells in primary visual cortex [7], as well as in V4 [8]. Finally, it has recently been shown that gain control schemes are needed to control behavior in invertebrates [9]. Despite the clear importance of gain control in a variety of neural computations, the biophysical mechanisms that support specific gain control mechanisms have been elusive [10][11][12][13][14][15][16][17][18][19][20].
Noise induced phenomena in nonlinear systems are a rich avenue of study [21], with recent interest on the impact of fluctuations on excitable systems, such as neurons [22]. Chance et al., Doiron et al.,and Hô & Destexhe [11,16,23] all report that an increase in the fluctuations of background conductance inputs results in a decrease of the overall gain of the transfer between a static driving input and the mean output firing rate. In particular, if the balance between background excitation and inhibition is carefully controlled [11], then the gain control is purely divisive (or multiplicative). This means that an increase in conductance fluctuations acts to scale the transfer function over a large range of input by a simple constant multiplier (,1). Related work has further explored the impact of fluctuations on spike response [13][14][15]20,24,25], with a the manipulation of the neural transfer function by background fluctuations being a central focus.
These studies address the gain control of a transfer function where the signal is either static or statistically stationary and the neural output is the time averaged firing rate. However, many neural coding tasks involve the processing of time-varying, high frequnecy stimuli. In these situations neural response are often transient, and a quasi-static approximation of input-output transfer fails to capture the actual spike response. For example, in the rodent vibrissa sensory [26], auditory [27][28][29], and electrosensory systems [30] stimuli and responses modulate on the order of a few milliseconds, i.e., on the order of, or even faster, than typical membrane time constants of neurons. Even in the visual system, where the relevant timescales of natural scenes are much slower, the response precision of thalamic neurons is at the millisecond level, and standard static transfer function analysis fails to capture neural response [31], yet contrast induced gain control persists [32]. In this study, we address the question of whether the fluctuation induced gain control mechanism explored for static transfer [11,16,23] can be operative for dynamic stimuli as well.
Any theoretical treatment of this problem requires 1) a framework accurately capturing the time varying spike response owing to time varying input statistics (e.g. temporally inhomogeneous input and output statistics), and 2) sufficient biophysical detail to incorporate conductance based synaptic inputs within spike creation. A useful tool for incorporating these two features into neuron models is the population density method [33][34][35][36][37][38][39]. In particular, Nykamp & Tranchina [35] have developed a simple one-dimensional population density method of conductance based leaky integrate-and-fire models (LIF). The one-dimensional version of the population density method allow us to easily study the firing rate responses to dynamic stimuli in the conductance based formalism of Chance et al. [11]. Minor differences in our proposed model and their dynamic clamp experiments to mimic conductance based inputs are presented in the discussion section. We first show that divisive gain modulation of the steady-state responses only hold for low output firing rates, in particular, where neurons are in the classical subthreshold regime. Second, when restricted to this regime we find the transient responses to dynamic stimuli, which can differ greatly from the quasi-static equilibrium response, also exhibit divisive gain modulation via fluctuation background conductances with the same scaling factor as computed in the static case. Thus, the divisive gain modulation proposed by Chance et al., Doiron et al.,and Hô & Destexhe [11,16,23] generalizes to the dynamic situation in a very natural way.

Integrate-and-fire neuron
We consider a leaky integrate-and-fire neuron (LIF) driven by a pre-synaptic population of excitatory (e) and inhibitory (i) cells. The neuron's voltage change is given by a random differential equation: Dividing by the leakage conductance g l yields: G e t ð Þ~X ke A ke d t{T ke ð Þ , where t m~C =g l is the membrane time constant, A k e=i is the random size (for simplicity, chosen from the same distribution) of the excitatory/inhibitory k th synaptic event. The arrival times T k of both excitatory and inhibitory synaptic inputs are governed by modulated Poisson processes with mean rates n e t ð Þ and n i t ð Þ, respectively. Throughout, e r is the resting membrane voltage, e e is the excitatory while e i is the inhibitory reversal potential. When the neuron's voltage crosses u th , a spike is recorded and the neuron enters a refractory period for a fixed time of t ref , after which, its voltage is reset to u reset . Consequently, the neuron's voltage u varies between e i and u th (e i ƒuƒu th ). Throughout this paper, we will set u reset~er , e r~{ 70 mV, e i~{ 80 mV, e e~0 mV, u th~{ 55 mV, t ref~2 ms, and t m~2 0 ms in accordance with estimates from experimental measures. We choose the average value of the random variables m Ae and m Ai so that the neuron's voltage change (from u {~e r ) is 60.5 mV [40]. The random variable A e=i has a parabolic distribution function with finite The background synaptic activity is balanced [11,40], namely n 0 e and n 0 i are chosen so that, in the absence of the driver input (n d~0 ), the random target voltage will have mean equal to the resting potential: m V~er . This will be true if: Of interest are the output threshold crossing times, and we estimate response statistics by combining the responses from N trials where the arrival times T k are statistically independent across trials yet share the same generating intensities n e t ð Þ and n i t ð Þ. The instantaneous firing rate of the neuron is defined as where t ij is the j th threshold crossing recorded during trial i. Throughout the paper we are interested in the relationship between the driver n d t ð Þ and the response r t ð Þ, and specifically how the balanced activity n 0 e and n 0 i can modulate the relationship.

Author Summary
Many neural computations, including sensory and motor processing, require neurons to control their sensitivity (often termed 'gain') to stimuli. One common form of gain manipulation is divisive gain control, where the neural response to a specific stimulus is simply scaled by a constant. Most previous theoretical and experimental work on divisive gain control have assumed input statistics to be constant in time. However, realistic inputs can be highly time-varying, often with time-varying statistics, and divisive gain control remains to be extended to these cases. A widespread mechanism for divisive gain control for static inputs is through an increase in stimulus independent membrane fluctuations. We address the question of whether this divisive gain control scheme is indeed operative for time-varying inputs. Using simplified spiking neuron models, we employ accurate theoretical methods to estimate the dynamic neural response. We find that gain control via membrane fluctuations does indeed extend to the time-varying regime, and moreover, the degree of divisive scaling does not depend on the timescales of the driving input. This significantly increases the relevance of this form of divisive gain control for neural computations where input statistics change in time, as expected during normal sensory and motor behavior.
1847.40 s 21 , respectively), which we respectively label low (black), medium (red), and high (blue). Chance et al. [11] modified the background level by various rate factors and labeled the regimes 1X, 2X, etc., which is slightly different than our convention of low, medium, and high. However, the resulting steady-state input/output curves ( Fig. 2) are similar to those in Chance et al. [11]. Also, our results below hold equally well for many other sets of balanced background activity. For a particular background intensity (low in this case) with random excitatory drive n d t ð Þ, the output is random (see spike raster plots). As the balanced background activity is increased, the variability in the voltage also increases (Fig. 1B).

Population density approach
A Monte Carlo simulation of Equation (2) would be computationally expensive to ensure an accurate result. In many studies only qualitative effects are reported, and thus quantitative accuracy is not at a premium. However, in our study the accuracy demands are large, as we will quantitatively compare the time dependent r t ð Þ for various levels of background intensities. To overcome the errors inherent in finite data from Monte Carlo simulations we use population density methods [35], known to give very accurate estimates of r t ð Þ (Fig. 1E) for the idealized neural models described by Equations (2)- (5).
In the population density method, neurons with similar biophysical properties are grouped together, and the evolution of a density function r u,t ð Þ is considered. In brief, r u,t ð Þ describes the voltage probability density over many statistically independent neurons. Integrating the density over a region in state space gives the probability that a neuron randomly chosen from the population will be in that region of state space: Let J V u,t ð Þ denote the probability current; a signed quantity with the convention that positive/negative J V is the probability per unit time of crossing u from below/above. The evolution of r u,t ð Þ is governed by a continuity equation [35]: We separate the probability current into three distinct terms: The first term, J l u,t ð Þ represents the deterministic leak to rest in the absence of synaptic events. The second and third terms, J e u,t ð Þ and J i u,t ð Þ, model the excitatory and inhibitory synaptic input driving the population. Mathematically we have: whereF F C Ã e=i is the complementary cumulative distribution function: With the chosen distribution for A e=i , the functions above are (setting x~{t m log 1{y ð Þ):F Finally, the instantaneous firing rate, r t ð Þ, is the flux of probability current through u th from below: The firing dynamics are implemented by an absorbing boundary condition at spike threshold, r u th ,t ð Þ~0, and the source term in Equation (7), modeling membrane reset after a refractory delay. The population average firing rate r t ð Þ by the population density method (see Eqs (6)-(8)) is a computationally efficient way of capturing r t ð Þ, compared to computationally expensive Monte Carlo simulations.

Divisive gain modulation
Gain modulation is typically studied in the equilibrium regime [10,11,16], where the driver input is constant in time n d t ð Þ~n d and the response r eq n d ð Þ~lim t?? r n d ,t ð Þ denotes the equilibrium firing rate as a function of input. A divisive gain modulation for where r eq j n d ð Þ is the response of the population with some background activity 'j' (for our purposes j = 1 is low, j = 2 is medium, and j = 3 is high background activity) and c j is a scalar. To measure the divisive gain modulation of a response r eq j v d ð Þ we fit the scaled response curves to r eq 1 u d ð Þ (low) by minimizing the mean-squared error E eq j : Finding the c j that minimizes E eq j is easily obtained by orthogonally projecting onto the subspace spanned by r j n d ð Þ: Because we want the largest possible range for divisive gain modulation, we steadily increase the maximum of n d n m ð Þ until the error in (10) becomes significant (i.e., the scaled curves no longer lie on top of each other). Let that maximum n d value be n max .
The nonequilbirum response r t ð Þ to a time-varying input n d t ð Þ is given by Equation (8). We extend divisive gain modulation to the The response curve r eq n d ð Þ for different levels of balanced background activity (n 0 e is determined once n 0 i is specified, see Equation (4)). For high output firing rates outside of the boxed region, the slopes of all of the response curves are almost equal. (B) Top Panel: zoomed region of (A). Bottom Panel: a least squares fit was used to find the scalar c eq j (see Equation (11)), and we plot the scaled response c eq j r eq j . Here, c eq 2~1 :72 and c eq 3~3 :2. doi:10.1371/journal.pcbi.1000365.g002 nonequilbrium setting, with the analogous description: The degree of divisive gain modulation in the nonequilibrium setting is determined by how well the time-varying responses scale with one another over a range of t. Thus, the scalar c j is now computed from integrals over t, compared to integrals over n for the equilibrium case (11).

Results
Divisive gain modulation in the steady state for subthreshold firing rates We compute the equilibrium input/output relationship, r eq n d ð Þ, using the same framework as Chance et al. [11] (section: Integrateand-fire Neuron), in hopes of first reproducing their results. Using balanced excitation and inhibition to mimic background synaptic activity, we compute r eq n d ð Þ for different fixed values of the driving input n d (Fig. 2). To obtain computationally accurate results in reasonable time we employ population density methods (section: Population Density Approach) to estimate the response r eq . Divisive gain modulation via increases of the background rates n 0 e and n 0 i occurs in the low firing rate region (boxed region of Fig. 2A). In this regime the neuron response is dominated by input fluctuations rather than any intrinsic spike rhythmicity, thereby replicating the high variability observed in the spike responses in cortical networks [40]. Throughout we refer to this as the fluctuation driven regime. In the fluctuation dominated regime the responses can be scaled, in the sense of Equations (10)- (11), to nearly quantitatively match one another (Fig. 2B). This rescaling of the response by background fluctuations qualitatively matches the results presented in [11].
In contrast, for very high output firing rates divisive gain modulation does not occur. The responses n eq n d ð Þ are nearly linear with very similar slopes ( Fig. 2A), showing only a background activity induced translation of the response (often termed subtractive gain modulation [10]). This region corresponds to a regime where input fluctuations have limited impact and the neuron response is predominately determined by the mean value of the input rates, and we refer to this as drift dominated regime. Fluctuation induced divisive gain control restricted to low firing rates is consistent with [16], where simulations of a large-scale compartmental neuron model were used. The insensitively of n d to input fluctuations at large n d has also been recorded in pyramidal cells and fast-spiking interneurons [41,42]. However, the exact n d where gain manipulation changes from divisive to subtractive (as n d increases) is difficult to compute and is often model specific [14]. Indeed, there are neurons where the influence of noise persists at high firing rates, such as in layer 5 of rat medial prefrontal cortex [24], however, the biophysical mechanisms that support this effect are absent in the standard LIF model.
In summary, population density methods (section: Population Density Approach) can replicate fluctuation-induced divisive gain modulation of the equilibrium response at low firing rates, previously observed in: simple integrate-and-fire models [14][15][16], simulations of biophysical realistic cell models [16,20], as well as simulated conductance experiments in vitro [11].

Divisive gain modulation with dynamic stimuli
We study the influence of background fluctuations on the nonequlibrium response to a highly time-varying excitatory drive. We choose an input rate n d t ð Þ consisting of sums of sinusoids with various amplitudes, phases, and frequencies to mimic 'rich' time varying stimuli (for an example see Fig. 3A). This produces an inhomogeneous Poisson process driver input n d t ð Þ, resulting in a non-stationary in time stochastic driving current. The response r t ð Þ inherits the non-staionarities of n d t ð Þ and is temporally modulated (Fig. 3B, black curve). Even though the stimulus results in a rather narrow range of response firing rates r t ð Þ, it has adequately rich temporal modulation to produce output firing rates that are different than the quasi-static response (Fig. 3B, brown curve), obtained by setting r t ð Þ~r eq n d t ð Þ ð Þ. The main result of our study is that fluctuation-induced divisive gain modulation is robust for low to moderate output firing rates in response to dynamic stimuli, despite the complicated dynamics of the leaky integrate-and-fire neuron in the nonequilibrium regime (i.e L Lt r u,t ð Þ=0). To demonstrate we compute the nonequilibrium responses (r 1 t ð Þ, r 2 t ð Þ, and r 3 t ð Þ), for the three levels of background activity used for the equilibrium case (low, medium, and high). For larger background activity the overall response is reduced, observed here since r 3 t ð Þƒr 2 t ð Þƒr 1 t ð Þ for all t (Fig. 3C). We compute the dynamical analogue of c eq j , c j (see Equations (11) and (14)) and the scaled response c j r j t ð Þ, which quantitatively matches the base response r 1 t ð Þ (Fig. 3D). This mimics the results for the equilibrium case (compare Figs. 2B and 3C and 3D). It is, a priori, unexpected that the dynamic response r t ð Þ (with timescale t m and refractory period t ref ) should scale in the same way as the equilibrium response r eq . Previously, Holt & Koch [10] showed that an increase in membrane shunting without a change in input fluctuation causes a translational shift, rather than division of the equilibrium response curves, which was also verified by Chance et al. [11]. To verify that a pure shunting change cannot result in divisive gain modulation of the nonequilibrium responses (Fig. 3D), we fix the background fluctuation level and driver n d t ð Þ, but increase the deterministic leakage conductance g l (Equation 1) to mimic different background synaptic activity (conductance) levels. Equivalently, t m is replaced with a scaled version: t m 1zn 0 e m Ae zn 0 i m Ai À Á , which has the same mean conductance in the absence of driving input as (n 0 e ,n 0 i ). In simulations where the background activity is set by deterministic leak rather than by synaptic conductance fluctuations, the neurons had negligible firing rates because they were unlikely to fire by random chance and did not scale in a divisive manner (not shown). For exposition, we set all of the background fluctuation levels to that of low and vary g l , and hence t m , to mimic deterministic effects of changing background activity so that there are less fluctuations, but still some amount to induce background firing. The unscaled responses (Fig. 4A) were scaled via a least squares fit (Equation (14)). Not surprisingly, the responses do not scale in a divisive manner (Fig. 4B). Thus, divisive gain modulation in the nonequilibrium regime critically depends on changing the background fluctuation levels. We remark that Chance et al. were in the high conductance state when they verified this whereas our regime has less overall conductance.
When the dynamic stimuli are increased so that resulting output firing rates are larger, the neurons no longer exhibit divisive gain modulation. Increasing the overall intensity of the driving input n d t ð Þ (compare Fig. 5A with Fig. 3A) yields firing rates r t ð Þ that are an order of magnitude larger (compare Fig. 5B with Fig. 3C). Increasing the overall background activity reduces the overall response magnitude (Fig. 5B), similar to what is observed in both the equilibrium and nonequlibrium regimes. However, when the response curves are scaled by c j computed for the low rate case pure divisive gain modulation is not observed for the high rate response. There is no trivial (a time independent c j ) or natural way to scale the output firing rate curves so that they lie on top of each other. Since divisive gain modulation does not hold in the equilibrium setting for high output firing rates (drift dominated regime), one would expect that it does not hold in the nonequilibrium state. However, both the equilibrium and nonequilibrium states are quite different and we present the failure of fluctuation induced division for the sake of completeness. It is interesting to note that for periods of time when the output firing rates are low, divisive gain modulation appears evident, likely owing to a transient excursion into the fluctuation driven regime.
In our model, when the driving input rate n d t ð Þ is low, the population of neurons rarely fire action potentials (i.e., low spontaneous activity). The firing rates in our simulations in this state range from nearly 0 to 3 s 21 , depending on the background level of activity. Although extracellular recordings in the cortex suggest the neurons can fire spontaneously at rates larger than 2 s 21 [43], such experiments are usually biased towards active neurons. Extracellular recordings by [44] that were unbiased towards responsive neurons suggest that many neurons have low spontaneous firing rates and that only a small fraction of neurons respond 'well' to stimuli in unanesthetized animals; this fact was also discussed in [43]. Moreover, calcium imaging experiments of awake and anesthetized rats in layer 2/3 of the cortex show that many neurons have resting firing rates less than 1 s 21 [45]. The actual firing rate of neurons in the resting state is a contentious issue, but our results hold for many parameter regimes (see Fig. 6).
Divisive gain modulation with dynamic stimuli is robust in the subthreshold regime (Fig. 6). To illustrate this point, the response r 1 t ð Þ with low background level to a time-varying driver input and the response r j t ð Þ to the same driver input with a second level of background activity are computed. We plot the logarithm of the area (or error E j , see Equation (13)) between the time-dependent response scaled by the equilibrium scale factor c eq j r j t ð Þ and the A logarthmic scale was used (Fig. 6) to better highlight the variety of values. The driver input is scaled as follows: where n t ð Þ is the driver input used previously (see Fig 3A), a is the scaling parameter, and m is a parameter that insures n d t ð Þ in positive and not too small. Notice a~1 corresponds to the driver input in Fig. 3A, and a~1:8 with m~1550 s {1 corresponds to the driver input in Fig. 5A. The vertical dot-dashed line in black corresponds to 0 error because it is the reference background curve for a given a. The two points marked by stars (*) in Fig. 6A at a~1 and n 0 i~1 400, 1900 s 21 correspond to the difference in area between the curves in Fig. 3D, which is quite small. In fact, for a large region of parameter space, divisive gain modulation holds (any patch that is orange to blue in Fig. 6A). The two black circles (N) in Fig. 6A at a~1:8 and n 0 i~1 400, 1900 s 21 correspond to the difference in area between the curves in Fig. 5B. With larger a values, the neurons are in the drift dominated regime, and divisive gain modulation no longer holds, as expected (red regions in Fig. 6A).
The average (unscaled) time-dependent response of the neurons with the same parameters and driver inputs as Fig. 6A are plotted in Fig. 6B on a logarithmic scale: The response with the same driver input as before (Fig. 3A), but with increased leakage conductance g l to mimic various background synaptic activity (conductance) levels in a deterministic way. The balanced background fluctuations levels (n 0 e ,n 0 i )are the same (low) in all 3 curves. The black curve has t m~2 0 ms (same curve in Fig. 3), the red curve has t m~1 8 ms, and the blue curve has t m~1 6 ms.  The three stars (*) at a~1 are the average firing rates of the unscaled responses in Fig. 3C, and the three black circles (N) at a~1:8 are the average firing rates of the unscaled responses in Fig. 5A (bottom panel). The average firing rate gives a qualitative idea of how large or small the response is as a and n 0 are varied. For example, the average firing rate at the low level is about 2 s 21 but the firing rate response can be quite low and as high as 15 s 21 (see Fig. 3B). Thus, divisive gain modulation holds for many parameters in a variety of subthreshold regimes.
Comparison of optimal scaling factor of equilibrium and nonequlibrium responses A gain control scheme will be effective in unpredictable environments if it is quantitatively insensitive to the timescales of the input, or in other words the degree to which the response is scaled should not depend on the spectral content of the signal. For fluctuation induced gain control we then require that the scaling factor c j associated with a specific background state would need to be independent of the temporal frequencies in the driver input n d t ð Þ. To test this we compare the optimal scaling factor c between two dynamical responses (each with a distinct balanced background state) where the synaptic driving input is: Notice the specified n d t ð Þ varies from 0 to n max , so that synaptic input rates are non-negative. Let us denote c by c v ð Þ for two given background rates driven by sinusoidal input with frequency v in Hz (here v is not the conventional radian frequency). When n max is in a low range the differences between c v ð Þ and c 0 ð Þ are negligible over a wide range of v (Fig. 7). This result is robust for a range of background states (Fig. 7A-D). The quantitative match between the divisive scaling of equilibrium and nonequlibrium responses extends to more complicated temporal modulations of the driving input (Fig. 3A). Specifically, we find that c j &c eq j for the results shown previously (c 2 &c eq 2 &1:72 and c 3 &c eq 3 &3:2 in Figs. 2 and 3). Thus, fluctuation induced gain control is quantitatively insensitive to the timescales of the driving input.

Frequency response to weakly time varying inputs
To better describe the mechanism underlying fluctuationinduced divisive gain control in the nonequlibrium, we focus on a weak time modulation of the input drive and compute the linear frequency response [34,46,47]. The frequency response function gives the first order temporal modulation of output firing rate assuming the synaptic driving input consists of a large constant component and a small time-varying component: Here we have set n eq d to be some fixed driver synaptic input rate, making the overall time independent excitatory input n e~n eq d zn 0 e , while the inhibitory input is still n i t ð Þ~n 0 i . In total we then have the equilibrium state defined by the triplet B~n eq d ,n 0 e ,n 0 i È É , and the time dependent component of the driver simply ee i2pvt . Assuming e%1 we approximate: The time modulation of the response is characterized by Z v,B ð Þ, indicating how large or small the first order response is to time- Figure 6. Parameter space where divisive gain modulation extends to the nonequilibrium regime. (A) The logarithm of the area (or error) between the time-dependent response curve scaled by the equilibrium scale factor c eq j r j t ð Þ and the r t ð Þ curve for many background levels and many driver inputs: The a values on the vertical axis corresponds to a scaling of the driver input used in Figure 3A (see text for details). The two points at a~1 and n 0 i~1 400, 1900 s 21 marked by stars (*) correspond to the difference in area between the curves in Figure 3D, and the two black circles (N) at a~1:8 and n 0 i~1 400, 1900 s 21 correspond to the difference in area between the curves in Figure 5B. Any region with colors in the range of orange to blue correspond to parameters where divisive gain modulation persists. (B) The average (unscaled) time-dependent response of the neurons with the same parameters and driver inputs as (A) are plotted on a logarithmic scale. The three stars (*) at a~1 are the average firing rates of the unscaled responses in Figure 3C, and the three black circles (N) at a~1:8 are the average firing rates of the unscaled responses in Figure 5A  Our earlier results (Fig. 3) show that for the same driver input u d t ð Þ and different background inputs that r j t ð Þ&r 1 t ð Þ c j for some scaling factor c j . However, we know that in limit e?0 the equilibrium response also satisfies r eq j &r eq 1 c j (Fig. 2). Combining these two results, and neglecting the O e 2 À Á terms in Equation (15), where B j is a background state. Satisfyingly, when neurons are in the fluctuation-dominated regime Z v,B ð Þ j jdoes indeed multiplicatively scale in the same quantitative manner for different levels of balanced background synaptic input scale ( Fig. 8A and 8B). The phase component W v,B ð Þ is the same for all B and v tested (insert Fig. 8A) and hence can not change the response r t ð Þ for different B. Thus from the quantitative scaling match of both r eq and Z v,B ð Þ j jwe expect fluctuation-induced divisive gain control to extend to weak inputs.
We remark that the near exact scaling of Z v,B ð Þ j jin the high frequency range (vw10 Hz) is important; if the multiplicative scaling was only true in the flat region of Z v,B ð Þ j j(vv10 Hz) then fluctuation-induced divisive gain control would fail in the nonequlibrium, i.e. when the quasi-static approximation fails. To be more specific, lim v?0 Z 0,B ð Þ j jdr eq =dn d , if we neglect the fluctuations given by the driver Poisson process. Thus the scaling of Z v,B ð Þ j jfor vv10 Hz is completely explained by the scaling of the gain of r eq . However, multiplicative scaling for vw10 Hz ensures that fluctuation induced gain control will extend into the nonequilibrium regime, even though the quasi-static approximation fails.
When the neurons are in the drift dominated regime, the frequency responses does not scale in a multiplicative manner because there are resonant peaks at integer multiples of the steadystate firing rate [34,48,49], and these resonant peaks occur at different frequencies for various balanced background synaptic activity (see Fig. 5B). Thus, divisive gain modulation with dynamic stimuli cannot possibly occur. The frequency responses in the drift dominated regime are scalar multiples of each other up to 10 Hz, where there appears to be divisive gain modulation with the same equilibrium scaling factors (see Fig. 5B). As explained in the previous paragraph, frequency response for 0vvv10 Hz is equal to the frequency response for v~0. However, the multiplicative scaling breaks down for the same v range where the quasi-static approximation breaks down, meaning that for drift dominated responses any fluctuation induced gain control in the equilibrium regime will not transfer to the nonequilibrium response.

Discussion
Chance et al. [11] described a mechanism by which divisive gain modulation results from a balanced, fluctuating background synaptic activity which both shunts and linearizes the membrane to spike transfer. The response r eq j is a scaled version of a baseline condition r eq 1 , and the dividing factor c j is independent of the driver intensity u eq d . However, many stimuli induce input and output statistics which vary on the timescale of neural integration [26,27,30,31]. Extending fluctuation induced gain control to accurately divide the response to these inputs is not automatic, as the spike-reset and refractory dynamics significantly shape the response in the nonequilibrium regime to be significantly different than the quasi-static approximation. However, our results show that the fluctuation induced gain control does extend to the nonequilibrium regime, increasing the potential utility of this form of gain control in neural processing. Furthermore, establishing the independence of the scaling term c j from the timescale of the driver greatly simplifies the circuitry required to implement gain control. In its simplest scenario, the gain of the response is set by the background rates n 0 e and n 0 i which maintain their scaling effect despite processing unpredictable environments where inputs statistics can vary dramatically.
The analysis of the time dependent response for weak signals showed how a scaling of r t ð Þ is inherited from an equivalent scaling of r eq and the response function Z v,B ð Þ j jby fluctuating background conductances. The response to an input s t ð Þ of arbitrary strength and spectrum can be written using the Volterra expansion [50]: whereZ Z 1 t,B ð Þ is the inverse Fourier transform of the response function Z v,B ð Þ described in Equation (15). Fluctuation induced gain control extends well into the nonlinear regime, evidenced by the empirical agreement in regimes where r t ð Þ varies significantly about r eq (Fig. 3D). In this case, the influence of the higher order terms in the Volterra expansion are likely important. We conjecture that, within the fluctuation dominated regime, each response functionZ Z i t, . . . ,t i ð Þ ,B j À Á &Z Z i t, . . . ,t i ð Þ ,B 1 À Á c j , meaning that the multiplicative scaling extends, response function-byresponse function, analogously into the nonlinear regime. This scenario is opposed to the one where each term exhibits scaling with distinct terms, yet the sum of terms somehow scales with c j , forcing agreement with our results where s t ð Þ has large temporal variance (Fig. 3D). In principle computingZ Z i t, . . . ,t i ð Þ ,B j À Á is quite difficult, however, if this scaling is correct then the influence of the stochastic background onZ Z i t, . . . ,t i ð Þ ,B j À Á , in the fluctuation driven regime, becomes straightforward.
Divisive gain control is a central tool in many neural computations [1], yet robust biophysical mechanisms that produce gain control are elusive [10][11][12][13]16]. Our work gives further evidence that using background fluctuations as a mechanism to scale responses is a surprisingly stable mechanism operable for a variety of input statistics. Fluctuation induced effects on the equilibrium state transfer different from divisive gain control have been reported [24,25]. Notably, [24] have shown that the firing response of pyramidal cells in layer 5 is sensitive to fluctuations at high rates, where the mean current no longer determines the spike rate. The mechanisms responsible are not present in the standard LIF model, however, modifications could possibly be made to model these effects and a population density equation could, in principle, be derived. These models would require more state variables and/or equations and in general are not computationally tractable without some reduction or approximation. Sophisticated methods for other neuron models have been developed [51][52][53]. Extending gain control to nonequilibrium responses to a larger class of models is currently an open avenue of research.
The LIF model we have used is an approximation to the dynamic clamp experiments of Chance et al. [11]. One difference is that our model does not have temporal correlations in the synaptic conductances, while there are temporal correlations in the experiments even though Chance et al. average over time (and trials) to obtain the firing rate. Also, we are using a simple yet biophysical spiking neuron model, where the level of background activity determines the variance of the background voltage (see Fig. 1B), consistent with the observations that membrane potential variability changes with the internal brain state [54]. In Chance et al. the variance of background voltage was the same for all background fluctuation levels, ensuring that the variability of the output firing rate is constant. Despite these differences, our results suggest that fluctuation induced divisive gain modulation is viable with dynamic stimuli.
The population density equations (6)-(8) that characterize the LIF model contain a partial differential-integral equation that is difficult to analyze. Our model is more general than white noise models that have an advection/diffusion density equation (e.g, Fokker-Planck equation) because it allows for large voltage changes upon receiving synaptic input events. However, the simulations shown in this paper are in the regime where the diffusion approximation is good. If the voltage change upon receiving synaptic events (excitatory or inhibitory) is assumed to be small, a good diffusion approximation of (6)-(8) is obtained by replacing r u',t ð Þ with r u,t ð Þz Lr u,t ð Þ Lu u'{u ð Þ in the integrals in the probability current terms J e u,t ð Þ and J i u,t ð Þ (see Text S1 and Figure S1). With large voltage changes, a similar approximation can be obtained by a re-scaling of the equation around the (deterministic) mean. However, a direct comparison to the Fokker-Planck equation with white noise conductances still must be done numerically because of the conductance-based input (Text S1 outlines the Fokker-Planck approximation to the full density equation and Figure S1 shows the magnitude of the advection/ diffusion coefficients). Moreover, the analytic formulas obtained with advection/diffusion equations are often computed numerically and usually assume at least a quasi-static approximation. With Poisson current injection however, a closed form Fokker-Planck approximation is obtained with drift and diffusion coefficients that can be written exactly in terms of voltage, input rates, and the statistics of C Ã e=i . An analytical explanation of the robust scaling of the firing rate responses that is observed remains elusive yet is conceivable because of the many analytical results obtained for density equations of a variety of neuron models [51,52]. However, we remark that even in the equilibrium regime an analytic explanation of divisive gain modulation via conductance fluctuations is difficult to obtain [14].