A Diffusive Homeostatic Signal Maintains Neural Heterogeneity and Responsiveness in Cortical Networks

Gaseous neurotransmitters such as nitric oxide (NO) provide a unique and often overlooked mechanism for neurons to communicate through diffusion within a network, independent of synaptic connectivity. NO provides homeostatic control of intrinsic excitability. Here we conduct a theoretical investigation of the distinguishing roles of NO-mediated diffusive homeostasis in comparison with canonical non-diffusive homeostasis in cortical networks. We find that both forms of homeostasis provide a robust mechanism for maintaining stable activity following perturbations. However, the resulting networks differ, with diffusive homeostasis maintaining substantial heterogeneity in activity levels of individual neurons, a feature disrupted in networks with non-diffusive homeostasis. This results in networks capable of representing input heterogeneity, and linearly responding over a broader range of inputs than those undergoing non-diffusive homeostasis. We further show that these properties are preserved when homeostatic and Hebbian plasticity are combined. These results suggest a mechanism for dynamically maintaining neural heterogeneity, and expose computational advantages of non-local homeostatic processes.


Introduction
Nitric oxide (NO) is a diffusive neurotransmitter which is widely synthesized in the central nervous system, from the retina to the hippocampus [1,2]. Its properties as a small nonpolar gas molecule allows rapid and unconstrained diffusion across cell membranes, a phenomenon often called volume transmission [3]. An important role of NO signaling is to regulate neural excitability through the modulation of potassium conductances in an activity-dependent manner, effectively mediating a form of homeostatic intrinsic plasticity (HIP). Experiments characterizing this effect also demonstrated that NO-synthesizing neurons can induce changes in the excitability of neurons located up to 100 μm away [4,5]. These findings are corroborated by a recent study demonstrating neurovascular coupling mediated through activity-dependent NO diffusion [6]. We build upon these observations, postulating a general form of HIP mediated by a diffusive neurotransmitter such as NO which we will refer to as diffusive homeostasis. This contrasts with canonical models of HIP, here referred to as non-diffusive homeostasis, which assume that each neuron has access to only its own activity [7].
Theoretical studies of HIP have generally focused on its role in maintaining stable network dynamics [8,9]. It has also been recently demonstrated that HIP can improve the computational performance of recurrent networks by increasing the complexity of network dynamics [10]. However, little is known about the effects of HIP on the heterogeneity typically observed in cortical networks; in particular, a growing body of evidence supports the finding that even neurons of the same type have a broad and heavy-tailed distribution of firing rates [11]. Rather than an epiphenomenon of biological noise, neural heterogeneity has been proposed to improve stimulus encoding by broadening the range of population responses [12,13]. However, this form of heterogeneity is difficult to reconcile with canonical models of HIP, which generally suppress cell-to-cell variability [14]. While some degree of heterogeneity in populations of the same type of neuron may emerge naturally [15], we found that such independent sources of variability will generally limit the responsiveness of a network through neuronal saturation.
Using network models and dynamic mean field analysis, here we show that networks with HIP mediated by diffusive neurotransmission exhibit a very different and unexpected behavior. Firstly, we report that diffusive homeostasis provides a natural substrate for flexibly maintaining substantial heterogeneity across a network. Secondly, the resulting population heterogeneity enables linear network responses over a wide range of inputs. This not only improves population coding, but enables a good use of available resources by ensuring that all neurons remain functionally responsive to changes in network dynamics. Finally, we demonstrate that these effects are preserved in networks whose recurrent synaptic inputs undergo Hebbian plasticity.

Results
We investigated the effects of diffusive homeostasis in a recurrent neural network (Fig 1A, see Methods) with sparse and random connectivity, based on conventional models of cortical networks giving rise to asynchronous irregular spiking activity [16]. Each neuron received external input with a rate randomly drawn from a normal distribution. Ca 2+ influx during a spike triggered NO synthesis through nNOS activation (Fig 1B, see Methods). Periodic boundary conditions were used to simulate NO diffusion, as we are simulating a subsection of a cortical network embedded in a larger section with similar activity. This resulted in a 2D toroidal surface on which diffusion was simulated. Boundary conditions with concentration gradients fixed at zero were also tested, with no appreciable differences in the outcome. Each neuron's firing threshold θ i underwent modulation through negative feedback mediated by the concentration of NO (Eq 6 in Methods).
The effect of a diffusive neurotransmitter mediating HIP within the network was investigated by comparing two cases: first where NO was allowed to diffuse freely across cell membranes as observed experimentally [4], and second without diffusion such that intracellular NO (B) Intracellular homeostatic signals in a model neuron. Each spike triggers calcium influx, which leads to nNOS activation and NO synthesis. (C) Mean population firing rates for networks with diffusive or non-diffusive homeostasis after an increase in external input (red triangle). Spatial distribution of NO concentrations at different times across the network with diffusion are shown below. (D-E) Distributions of firing rates and log firing rates (insets) after homeostasis from network simulations (D) and mean-field analysis (E), both receiving independent Poisson inputs drawn from a Gaussian distribution. (F) Distributions of firing rates in the mean-field analysis for low and high covariance σ of threshold and input rate. concentration was affected only by a neuron's own recent activity (Fig 2). The latter corresponds to a canonical model of HIP as investigated before [8,9].
Diffusive homeostasis enables a broad firing rate distribution Fig 1C illustrates that both forms of homeostasis stabilized network activity following an increase in input. There was however a crucial difference in how the neurons reacted to this change. While for non-diffusive homeostasis each neuron simply returned to its target firing rate, diffusive homeostasis caused each neuron to sense a mixture of its own activity level and that of the rest of the network. This can be seen in the spatial concentration profiles in Fig 1C. It is important to note that the spatial position of each neuron was random and independent of its connections, meaning that there was no explicitly defined structure in the NO concentrations.
As a result, these networks exhibited a very different steady state behavior. The firing rate distribution was narrow as expected for non-diffusive homeostasis, but broad and heavy-tailed for diffusive homeostasis (Fig 1D). The latter is consistent with recent experimental results indicating that firing rate distributions in cortex are generally heavy-tailed, approximating lognormal distributions [11]. There were no noticeable differences in inter-spike interval statistics between networks with diffusive and non-diffusive homeostasis (not illustrated).
We investigated the difference in firing rate distributions by modeling the relation between activity read-out and homeostatic compensation in these two cases using a dynamic mean-field model (see Methods). This approach considered an unconnected population of neurons with random inputs, where each of the two scenarios was simulated by using an appropriate activity read-out. HIP was implemented as in the full spiking model, but the degree of diffusive signaling was now controlled by a single parameter, α (Eq 11 in Methods), which determined the balance between local and global activity read-out. If small, neurons used primarily their own activity to modulate their firing threshold, while increasing α caused the firing threshold to depend more strongly on the average population activity. Setting, for instance, α = 0.8 led to a broad and heavy-tailed rate distribution similar to the full model, while α = 0 yielded a narrow distribution as in the non-diffusive case (Fig 1E).
This model provides a simple and intuitive explanation for this effect. For a non-interacting population, non-diffusive homeostasis can be thought of as precisely matching a neuron's input μ i and its threshold θ i to maintain the target firing rate. We can imitate this by introducing a covariance σ(μ, θ) between μ i and θ i , such that a high input rate implies a high firing threshold and a low input rate a low firing threshold. Since setting α > 0 (analogous to diffusive homeostasis) introduces a correlation between a neuron's threshold θ i and the average population threshold y, this effectively results in a decorrelation of μ i and θ i in comparison with setting α = 0 (analogous to non-diffusive homeostasis). In line with the previous results, populations with for instance σ(μ, θ) = 0.6 yielded a broader and more heavy-tailed distribution of firing rates than populations with σ(μ, θ) = 0.99 (Fig 1F).
Since non-diffusive homeostasis directly relates the firing threshold of a neuron to its input, we observed a wider distribution of firing thresholds, which in turn ensured that all neurons assumed similar firing rates. Diffusive homeostasis, on the other hand, yielded similar firing thresholds across the population (Fig 3A and 3B). When combined with the nonlinear inputoutput relation of neurons [17], this gave rise to the broad firing rate distributions we observed (see also Discussion). This result was robust to changes in the rate of NO diffusion. While decreasing the rate of diffusion, D, did result in slightly narrower firing rate distributions, they were broader than in networks with non-diffusive homeostasis across a wide range of values ( Fig 3C). A similar trend was observed when varying the width of the external input rate distribution. While decreasing this width led to a decrease in the width of the firing rate distribution, they were consistently broader in networks with diffusive homeostasis (Fig 3D).
Since one may argue that diffusive homeostasis is merely adding variability to each neuron's homeostatic signal due to the influence of neighboring neurons' activity, we now ask whether it is possible to achieve broad firing rate distributions with non-diffusive homeostasis. Indeed, by introducing variability of homeostatic targets (see Methods), we could produce a distribution of firing rates similar to that observed with diffusive homeostasis (Fig 1D, red histogram). However, as we will show next, the effect of diffusive homeostasis is quite distinct from that of activity-independent, 'quenched' heterogeneity arising from randomly distributed homeostatic targets.

Diffusive homeostasis retains input heterogeneity
To investigate the functional consequences of heterogeneity caused by a diffusive homeostatic process, we next simulated specific changes in external input. First, we stimulated small random groups of neurons at higher input rates of 5 Hz and 10 Hz (versus a baseline of 2.5 Hz), as illustrated in Fig 4. Such inputs may, for instance, reflect developmental or other plastic changes that lead to a long-lasting change in network input. In these simulations, the average network firing rate was reliably brought back to the original target firing rate by both forms of homeostasis ( Fig 4A-4C, black traces). As above, in networks with non-diffusive homeostasis this was achieved by returning the rate of each neuron to the target firing rate regardless of their external input (Fig 4A, colored traces). In contrast, for networks with diffusive homeostasis, we found that the separability of firing rates of individual groups are maintained according to their input, while the firing rates of all groups were simultaneously reduced so that the average network firing rate again reached the target (Fig 4B, colored traces). Introducing variability in homeostatic targets for the non-diffusive case, as described previously, did not maintain separability of individual groups as in the diffusive case. Instead, the different groups returned to their mean firing rates that existed before inputs were elevated ( Fig 4C).
The distribution of final firing thresholds explains these differences (Fig 4D-4F). For nondiffusive homeostasis, neurons in the group receiving 10 Hz input had the highest thresholds since they needed to reduce their firing rate the most, followed by the 5 Hz and 2.5 Hz groups respectively. This led to the final threshold of each neuron reflecting its input. Note that the distribution of firing thresholds is broader in this setup than in (Fig 3A and 3B), as a broader range of inputs is given to the network. For a diffusive signal, a neuron's firing threshold is modulated by the activity of nearby neurons. Since group membership of a neuron is   independent of its position, this effect again introduced a correlation between each neuron's threshold and the mean threshold of the entire network, resulting in a distribution of final thresholds which are less segregated according to their input compared to a network with nondiffusive homeostasis. Thus, firing thresholds in neurons undergoing diffusive homeostasis were more weakly related to their external input. This in turn preserves local firing rate differences in input groups while maintaining constant average network activity. Introducing variable targets for non-diffusive homeostasis caused the thresholds to depend more strongly on their external input, similar to the original non-diffusive case.
We could broadly reproduce the distinctions between diffusive and non-diffusive homeostasis in the dynamic mean-field approach by varying α. For α = 0, modeling non-diffusive homeostasis, we obtained identical firing rates in input groups, as in the recurrent network ( Fig  4G). Note that changing the input of groups of neurons in the recurrent network also affects the activity of neurons with fixed input (Fig 4B, red traces) due to recurrent connections, an effect that is obviously absent in the dynamic mean-field description. Increasing α led to local firing rate differences persisting for longer periods of time. However, these differences eventually decay very slowly, only remaining stable for the case where α = 1 (Fig 4H and 4I). The reason this occurs is, even after the population activity has quickly reached its homeostatic target, the deviations of the input groups still exert a small homeostatic force when α < 1. For example, if α = 0.95, there will be a relatively fast change in thresholds as the population activity reaches its target, followed by much slower changes, at 1 − α = 0.05 times the speed ( Fig 4I). This does not happen to the same extent in the spiking network simulations with diffusive homeostasis, as diffusion of NO ensure that deviations from the population activity are directly compensated for by neighboring neurons. Differences persist for 3245 ± 440 s, compared with 115 ± 6 s and 140 ± 40 s for non-diffusive homeostasis with uniform and variable homeostatic targets, respectively (± symbol denotes standard error of the mean of 6 independent network realizations in each case, see Methods). Since we have increased the speed of homeostasis in order to reduce simulation time (see Methods), a more realistic time course of 15 minutes for NO modulation would cause input differences to persist in networks with diffusive homeostasis for many hours to days [5]. Taken together, this shows that diffusive homeostasis can retain input heterogeneity due to the influence of neighboring neurons' activity on an individual neuron's firing threshold.

Population heterogeneity during diffusive homeostasis enables linear network responses
In the simulations shown so far, each neuron received a static input throughout since we were interested in the final network states. We now investigate how these networks respond to fast changes in input; specifically how faithfully each neuron represents its change in input. Since networks with diffusive homeostasis simultaneously maintain constant average network activity and firing rate heterogeneity, we expected that this should allow input modulations to be followed more precisely due to a greater representational capability.
After the network reached steady state under an initial distribution of external inputs, we froze homeostasis so as to simulate fast changes in activity, since we assume that homeostasis is not active over these time scales. We then regenerated the external inputs to each neuron from the same distribution presented during homeostasis. This can be thought of as a re-configuration of inputs due to external fluctuations. To best represent such changes in a simple population coding paradigm, each neuron should respond linearly to a change in input; non-linear transformations may lead to an information loss and hence affect neural computations, although this may indeed be desirable in some brain regions. We interpreted the range of changes in input over which this response is linear, or non-saturating, as the range over which homeostasis does not interfere with the network response. Fig 5A-5C show the change in input rate versus change in output rate of each neuron. A highly nonlinear response was observed in networks with non-diffusive homeostasis, with rectification for large decreases in input and superlinear responses for large increases in input. This effect was quantified by an R 2 value of 0.57 from a linear regression. Conversely, networks with diffusive homeostasis exhibited a linear response across the entire range of input changes, with an R 2 value of 0.85. Population heterogeneity can also be achieved, as discussed before, by introducing target variability during non-diffusive homeostasis. This yielded a similar non-linear response as in the non-diffusive network with homogeneous targets, with an R 2 value of 0.38.
A consequence of the asymmetry in responses to input changes for networks with non-diffusive homeostasis was that the population rate increases upon regenerating inputs, despite the fact that mean input to the network remained unchanged (Fig 5D). This did not occur for networks with diffusive homeostasis, suggesting that these networks are more adept at maintaining a target level of activity in conditions where external inputs are dynamic and fast-changing. Crucially, the benefits of a diffusive homeostatic signal can be achieved by a relatively broad range of values for the rate of diffusion, D, indicating that the effects we describe are robust to precise parameter choices ( Fig 5E). Increasing the rate of NO decay, λ (see Methods), has a similar effect to decreasing D (S1 Fig). This difference in responses to input changes could again be reproduced in the dynamic mean-field approach. This allowed us to characterize population responses across different effective diffusive ranges, using the R 2 value from a linear regression as a measure of response linearity. Fig 5F shows R 2 values across a range of different input distribution widths, δ, as α is varied to model different diffusion coefficients (see Methods). This revealed a dependency on δ: While values of α * 1 exhibited the best response for smaller δ, hence cases where the inputs are rather narrow, the optimal α decreased as δ increased, as well as the overall response linearity. This dependence on input width can be explained by considering the manner in which a population of neurons with a distribution of dynamic ranges span a range of inputs. If this range of inputs is small, then all neurons will span it regardless of their dynamic range (determined by their firing threshold), hence the high values of R 2 for δ = 0.1. For an intermediate range of inputs, neurons whose dynamic range is best adapted to the average input are most responsive. This is achieved by increasing α. If the range of inputs is very large (δ = 1.0), R 2 values are low since the dynamic ranges of the population cannot span the inputs. This effect is stronger at high α, as firing thresholds are more correlated, and the dynamics range of most neurons cannot capture the full input variance.
Since connection probability falls off with spatial distance in cortical networks [18], we additionally simulated recurrent networks featuring such connectivity profiles. These networks exhibited qualitatively similar behavior under diffusive and non-diffusive homeostasis compared to networks without any spatial dependence in connectivity (S2 Fig).
Up until this point we have presented diffusive and non-diffusive homeostatic mechanisms as dichotomies, which has enabled a clear investigation of their distinct effects on network properties. However, it is more biologically relevant to investigate networks in which both mechanisms are simultaneously active. Fig 6 shows that the increased neural heterogeneity and response linearity observed in networks with diffusive homeostasis are also present in networks with both diffusive and non-diffusive homeostasis, and that the degrees of neural heterogeneity ( Fig 6A) and response linearity (Fig 6B) are determined by the relative timescales of these mechanisms (see Methods). As the ratio of timescales of non-diffusive homeostasis to diffusive homeostasis is increased (i.e. as non-diffusive homeostasis becomes slower than diffusive homeostasis), the network goes from narrow steady state firing rate distributions to broad, and exhibits an increase in response linearity, thus becoming more similar to networks with only diffusive homeostasis.
Overall, these results suggest that networks undergoing diffusive homeostasis are better suited to linearly represent a range of inputs. We investigated this by presenting the networks with time-varying inputs after freezing homeostasis. Groups of excitatory neurons received additional inputs which were randomly and independently generated after fixed time intervals (see Methods). Fig 7A shows the representation of such a time-varying input pattern (dotted black line) for each network (colored lines). Networks which have undergone diffusive homeostasis were capable of tracking this input significantly better than their non-diffusive counterparts, as characterized by the RMS error between the network response and input pattern (0.12 for diffusive homeostasis; 0.23 and 0.19 for non-diffusive homeostasis with uniform and variable targets, respectively; Fig 7B).
We can explore these differences further by constructing a simplified task in which a population of orientation-selective neurons respond to the orientation of a stimulus (see Fig 7C-7F, Methods). This is not intended to represent circuits which perform this task in the brain, but to serve purely as a demonstration of the relative merits of linear and non-linear network responses.
Neurons in the network are randomly assigned a preferred stimulus orientation, independent of their spatial position. A stimulus of a certain orientation can then be presented to the network by varying the external input rates of each neuron, with neurons whose preferred orientation is closest to the stimulus orientation receiving the highest input rate. The stimulus orientation can be decoded from the network by taking the vector average of the stimulus response across all neurons. The orientation of this vector average, or population vector, is the decoded stimulus orientation. Networks with linear responses perform better than those with non-linear responses in decoding stimulus orientation, as measured by the standard deviation of errors in the orientation of the population vector compared to the stimulus orientation (41°, 63°and 72°for diffusive homeostasis, non-diffusive, and non-diffusive with variable targets respectively, Fig 7G).

Properties of diffusive homeostasis are conserved in networks with Hebbian plasticity
In the networks described so far, we have used static and uniform synaptic weights for recurrent connections. We next considered whether the observed properties of diffusive homeostasis are altered by the presence of plastic synaptic weights, in particular when Hebbian spike-timing-dependent plasticity (STDP) is introduced (see Methods). Using a standard model of STDP with additive depression and potentiation for all recurrent excitatory synapses, we simulated networks with both STDP and homeostasis active until synaptic weight and firing rate distributions reached a steady state [19]. As before, firing rate distributions were broader in networks with diffusive homeostasis (Fig 8Bi). Broad distributions could also be achieved by introducing variability in homeostatic targets. Spiking activity remained asynchronous after STDP, as shown by the distribution of inter-spike intervals and the spike autocorrelograms, although STDP caused weakly synchronous activity in networks without any form of homeostasis (Fig 8Bii-8Biv) [20]. The additive STDP rule led to a bimodal distribution of synaptic weights (Fig 8A), as previously reported [19].
STDP amplified the differences in response linearity that were observed between homeostatic cases. Inputs to each neuron were regenerated from the same distribution presented during plasticity, and the corresponding change in output rate was compared to the change in input rate, as in Fig 5A-5C. While the response linearity, given by the mean R 2 value, for networks with diffusive homeostasis was 0.16, networks with non-diffusive homeostasis exhibited much lower mean values of 0.01 and 0.02, for uniform and variable homeostatic targets respectively (Fig 8C). Networks without any homeostasis had a mean value of 0.1. These R 2 values were lower than those from networks without any STDP (Fig 5E), which was likely due to a combination of weaker external input given to these networks (g ext of 40 ns compared to 80 ns, see Methods), and stronger recurrent excitation received by neurons in these networks due to potentiation during STDP. This was tested by assessing response linearity in networks with static weight matrices obtained by shuffling the steady-state weight matrix of a network which had undergone STDP without any homeostasis (see Methods). These networks were run for each different homeostatic case until a steady state was reached, and inputs to each neuron were regenerated in order to measure response linearity. While shuffling synaptic weights did increase R 2 values across all networks (Fig 8C, crosshatched bars), indicating that STDP plays a role in decreasing response linearity, they remained lower than in networks from Fig 5E, confirming that the reduced influence of external input compared with recurrent input is largely responsible for this difference. We observed qualitatively similar retention of broad firing rate distributions and response linearity with diffusive homeostasis when a weight-dependent update rule was used (not illustrated), which has been argued to lead to more realistic weight distributions [21].

Discussion
Neural homeostasis is commonly thought of as a local process, where neurons individually sense their activity levels and respond with a compensatory change if activity changes. Here we investigated a complementary mechanism, where homeostasis is mediated by a diffusive molecule such as NO that acts as a non-local signal. Using a generic recurrent network model, we show that this form of homeostasis can have unexpected consequences. First, we found that it enables and maintains substantial population heterogeneity in firing rates, similar to that observed experimentally in intact circuits [11], and that input heterogeneities can be preserved in the population activity. Second, the specific form of neural heterogeneity brought about by diffusive homeostasis is particularly suited to support linear network responses over a broad range of inputs. It is important to note that this behavior differs from that of networks where heterogeneity is simply introduced by randomly assigning a different target to each neuron. These results predict that disrupting neural diffusive NO signaling can affect perceptual and cognitive abilities through changes of neural population responses. While other non-diffusive homeostatic mechanisms would continue to stabilize neural activity, the lack of a signal related to the average population activity may disrupt the flexible maintenance of firing rate heterogeneity, and as a result the ability to represent network inputs.
Mean-field analysis revealed that these differences are essentially due to the diffusive messenger providing each neuron with a combination of the average network activity and its own activity as the homeostatic signal. Diffusion of the signal from highly active neurons causes a reduction in the activity of their neighbors, such that firing rates of highly active neurons do not have to be completely reduced in order for the population to achieve a target rate. As a consequence, diffusive homeostasis furnishes a network with an efficient way of flexibly maintaining heterogeneity of firing rates. These effects can also be understood by considering the neural transfer functions, as illustrated in Fig 9A and 9B, which provides an intuitive explanation for the differences in firing rate distributions observed under diffusive homeostasis [17,22]. For non-diffusive homeostasis the transfer function of each neuron is brought to center around its input, leading to a narrowing of the firing rate distribution. Diffusive homeostasis decorrelates the input and threshold of individual neurons, resulting in a population of neurons residing along the entire transfer function. This preserves the non-linear shape of the transfer function, causing broad and heavy-tailed firing rate distributions.
Narrow firing rate distributions are an obvious consequence of local homeostatic processes, as for instance shown recently with homeostasis implemented as local synaptic metaplasticity [14]. This is in apparent conflict with the growing body of experiments documenting broad and heavy-tailed distributions of firing rates in cortex [11]. One could argue that a straightforward explanation is a process, for example genetic or developmental, which randomly assigns neurons heterogeneous homeostatic targets. While we show here that this can result in broader firing rate distributions, we also found that this generally leads to networks with a mismatch between the neural dynamic ranges and input statistics, which in turn limits the responsiveness of the network.
A striking feature of diffusive homeostasis is the lack of requirement for any such distribution of homeostatic targets, as the diffusive signal can be effectively exploited through providing a context for heterogeneity-neurons which maintain a significantly higher firing rate than the rest of the network also synthesize a higher level of the diffusive signal, thus ensuring that their deviation from the average firing rate is counterbalanced by lowering neighboring neurons' firing rates. This mechanism essentially allows neurons to differ in activity from the population as long as the population as a whole provides some compensation for these deviations. Moreover, this mechanism is compatible with the recent finding that a minority of cells were found to consistently be the most highly active and informative across brain states [23]. While non-diffusive homeostasis would have a disruptive effect on such a 'preserved minority' of neurons by reducing their activity towards those of the less active majority, diffusive homeostasis provides a substrate for maintaining their differentiated activity.
A significant distinction between the effects of diffusive and non-diffusive homeostasis appears when network responses to rapidly changing input are considered (Fig 5). We show that networks with diffusive homeostasis represent input changes more faithfully than those with non-diffusive homeostasis. Saturation of neurons' responses to large changes are observed in networks without diffusion-this effect is further illustrated in Fig 9C and 9D.
Across a spatially homogeneous network, diffusing signals act to effectively shift the transfer function of each neuron towards the average network input, ensuring that neurons are responsive across the entire range of inputs presented to a network. This is in contrast to networks with non-diffusive homeostasis, in which individual neurons are only responsive in a range around their current input. Moreover, the asymmetric response of networks with non-diffusive homeostasis causes the average network activity to increase after fast input changes, while it is constant for a network with diffusive homeostasis (Fig 5D). The latter case is consistent with observations that mean population firing rates are preserved across novel and familiar environments and across different episodes of slow-wave sleep [24,25] Networks with diffusive homeostasis have an improved ability to accurately track time varying inputs (Fig 7A and 7B) as a direct consequence of their linear responses. Beneficial effects of neural heterogeneity for population coding have been suggested before [13,26], but here we find that the broad linear response regime maintained by diffusive homeostasis further improves network performance. This improvement in network performance is also observed in a simplified stimulus orientation decoding task (Fig 7C-7G). Networks with diffusive homeostasis perform better than those with non-diffusive homeostasis when a population vector is constructed from the neural responses in order to decode stimulus orientation (Fig 7G). Although there exist alternative methods for decoding stimuli, the population vector has been shown to exhibit performance close to the optimal maximum likelihood procedure for broad tuning, as was used in our example [27].
These distinctions between diffusive and non-diffusive homeostasis are conserved in networks with STDP (Fig 8). This demonstrates that the limitations of non-diffusive homeostasis in maintaining neural heterogeneity and responsiveness extend beyond the case of static inputs, towards more realistic situations in which neurons receive ongoing and diverse perturbations. Indeed, networks with non-diffusive homeostasis lost almost all sensitivity to external inputs after STDP, while networks with diffusive homeostasis retained this sensitivity (Fig 8C).
The consequences of diffusive and non-diffusive homeostasis coexisting were also explored, by implementing these mechanisms simultaneously in a single network (Fig 6). Stable activity could be robustly maintained, with the resulting network behavior depending on the relative timescales of the non-diffusive and diffusive mechanisms. If non-diffusive homeostasis acted faster than diffusive homeostasis, the network exhibited a narrow rate distribution and a low responsiveness to input changes. Conversely, if diffusive homeostasis acted faster than non-diffusive homeostasis, the network exhibited broad firing rate distributions and linear responses to input changes. This is a plausible scenario, as NO modulation of ion channels occurs over a timescale of 15 minutes [5], while other homeostatic processes which require transcriptional changes occur over a timescale of hours to days [28]. These results reflect what is observed as α is varied in the dynamic mean-field analysis, as local and global homeostatic mechanisms are simultaneously active for values in the range 0 < α < 1.
It is important to note that modeling HIP as a force acting on the threshold of an integrateand-fire neuron in order to achieve a target firing rate is a significant simplification. More physiologically realistic descriptions of homeostatic processes reveal the complex relationship between ion channel concentrations and the regulation of a wide range of neural activity [28]. Moreover, a number of previous studies have explored the effects of volume transmission on network dynamics, including its potential in implementing a winner-take-all function [29], the ability of a diffusive signal to reflect the average activity of a group of neurons [30], and the role of another diffusive neurotransmitter, TNFα, in epileptogenesis [31]. Here, we add a functional interpretation by exploring its effects on neural heterogeneity and responsiveness within a network.
While NO is involved in a wide variety of neural processes throughout development and learning [32][33][34], these were ignored throughout for the sake of simplicity and tractability. Nonetheless, the impaired performance of nNOS knock-out mice in cognitive tasks [35] and the prevalence of epilepsy following nNOS inhibition [36] could be linked to diminished homeostatic control of neural excitability. Finally, the outcome of this study is not necessarily confined to NO, and could equally apply to other diffusive neurotransmitters observed in the brain such as hydrogen sulfide and carbon monoxide [37]. Indeed, we conclude that it demonstrates the potential role of diffusive neurotransmitters as an economical and reliable signal of activity across a population of neurons.

Network model
We simulated a spiking network of leaky integrate-and-fire (LIF) model neurons with conductance-based synapses and injected Ornstein-Uhlenbeck noise, as described by where v is the membrane potential, τ m the membrane time constant, E l the leak conductance reversal potential, and σ OU the variance of the injected noise. η(t) is an Ornstein-Uhlenbeck process with zero mean, unit variance, and correlation time τ OU = 1 ms [38]. g e and g i are the excitatory and inhibitory synaptic currents respectively, given by Eq (2), where t k denotes the time of all k incoming spikes. The reversal potential of the synapses are denoted by E e and E i , the synaptic conductances by J e and J i , and the synaptic time constants by τ e and τ i . The external input conductance is given by J ext , and t ext denotes the arrival time of external input, modeled as an independent homogeneous Poisson process for each neuron i with rate μ i . A spike is emitted whenever the membrane potential v exceeds the firing threshold θ, and the membrane potential is then reset to the resting potential value, v r , after a refractory period, τ ref .
The network was made up of N neurons; 0.8N excitatory and 0.2N inhibitory, with excitatory and inhibitory synaptic conductances scaled so that the network was in a balanced state [16]. Recurrent connections were random and sparse, with connection probability ¼ C N independent of neuron type, where C defines the mean number of synapses per neuron. The balanced state was achieved in the network through scaling the inhibitory synaptic conductances by a factor of g, such that J i = gJ e .

NO synthesis and diffusion
We assumed that neuronal NO synthase (nNOS) is activated by Ca 2+ influx following a spike Eq (3) and describe the relationship through the Hill Equation, Eq (4), which results in a sigmoidal concentration dependence, where n and K are parameters of the Hill equation [39] and τ Ca 2+ and τ nNOS are the timescales of Ca 2+ decay and nNOS activation respectively.
Throughout this paper we considered the case where all neurons, inhibitory and excitatory, express nNOS. The 2D diffusion equation, Eq (5), was solved numerically using a spatial resolution ds, diffusion coefficient D and a decay term λ [40]. Neurons were represented by a point source according to their activated nNOS concentration. This is a reasonable simplification, in particular as the networks we study are assumed to be large compared to the dimensions of individual somata. Periodic boundary conditions were used, as we assume we are simulating a subsection of a cortical network embedded in a larger cortical network with similar network activity. This resulted in a 2D toroidal surface on which diffusion was simulated. Boundary conditions with fixed concentration or fixed zero gradient were also tested, with no appreciable differences in the spatial distribution of NO. However, boundary conditions with nonzero fixed gradients of NO, approximating a large NO sink or source, can affect the spatial distribution of NO, and could hinder the ability of NO to act as a homeostatic readout. For further discussions of the diffusive properties of nitric oxide in neural tissue, see [41].
The homeostatic effect of NO was represented in neuron i by an increase in θ i , the firing threshold, according to the relative difference in intracellular NO concentration [NO] and a target concentration [NO] 0 ; where τ HIP is the timescale of homeostasis.

Non-diffusive homeostasis
For simplicity, the implementation of non-diffusive homeostasis is almost identical to that of diffusive homeostasis, in that the putative non-diffusive neuromodulator [NO non-diffusive ] is synthesized through Eqs (3) and (4), and modulates firing thresholds through Eq (6). The only difference is that the diffusion term in Eq (5) is removed, so that [NO non-diffusive ] is entirely determined by the rate of intracellular synthesis and decay;

Dynamic mean-field analysis
For a detailed derivation of equations used in our dynamic mean-field analysis, see [16] and [17]. Briefly, under the assumptions that the network is in an asynchronous regime and that a single EPSP is sufficiently small compared to the voltage required to elicit a spike from resting membrane potential, we can extract the mean firing rate of an LIF neuron in a recurrent network by solving a pair of equations under the condition of self-consistency. The synaptic current for a neuron i in a time interval τ can be described by its mean μ i and standard deviation σ i as follows: where J is the synaptic efficacy, C the number of synapses per neuron and ν the average population firing rate. The expected mean firing rate ϕ i (μ i , σ i ) of an LIF neuron with this synaptic current is given by where erfc is the complementary error function. Since the firing rate described by Eq (9) is determined by the synaptic current parameters μ i and σ i , which are in turn determined by the population firing rate ν, self-consistency requires that the rate which determines the synaptic current parameters must also be equal to the rate which is produced by these parameters, that is: We simulated a non-interacting population of neurons described by the mean-field theory, in which all neurons are identical except for their threshold θ i . Although there is no recurrent excitation within the population, the synaptic current statistics are comparable to that which a neuron within a recurrent network would receive. This enabled us to consider the firing rate distributions arising from presenting single neurons with distributions of synaptic currents, similar to the approach by [17].
We assumed that a neuron embedded in a homogeneous network receiving a diffusive homeostatic signal is analogous to a neuron using a combination of its own firing rate and the average population firing rate as a signal. The network can then be reduced to a population in which the firing threshold θ i of each neuron i is modulated according to where ϕ 0 is the target firing rate and ϕ i and are the firing rate of the neuron i and the population respectively. α was varied between 0 and 1 and can be thought of as the proportion of NO which a neuron receives due to diffusion from other neurons, with α = 0 indicating that each neuron senses only its own activity and α = 1 indicating that all neurons share an identical population-wide signal.
In order to implement homeostasis in this setup, we iterated through Eqs (8) and (9) until Eq (10) is satisfied to a precision of 10 −4 Hz, where Eq (9) returns ϕ i for each neuron in the population, and ¼ P i N is used as ν in Eq (8). At each timestep the thresholds θ i of each neuron were modulated according to Eq (11), and rates ϕ i were subsequently recalculated from Eq (9).

Procedure for Fig 1
External input rates μ i for each neuron i were randomly drawn from a Gaussian distribution such that μ i * N(10, 10 2 ) Hz. Since the mean NO concentration takes time to reach a steady state in the recurrent network simulations, we ran the network for 100 s without homeostasis and with all neurons receiving 5 Hz input, defining the target NO concentration [NO] 0 to be the mean NO concentration across all neurons at 100 s.
For the dynamic mean field analysis, we chose parameters which roughly match the rate statistics of the recurrent network simulations. Inputs to each neuron were drawn from a Gaussian distribution such that μ i * N(5.7, δ 2 ), s i ¼ ffiffiffiffi m i p . δ = 0.4 is the width of the distribution of mean inputs to the population. Note that the parameter δ referred to here differs from the δ in Eqs (1)-(3).

Adding target variability
In order to match the distribution of effective targets observed during diffusive homeostasis for networks with non-diffusive homeostasis, we assigned each neuron in the non-diffusive network a different homeostatic target, [NO 0 ] i . A network without any homeostasis is presented with input statistics (μ i * N(2, 5 2 ) Hz), tuned such that the firing rate distribution match that of the network with diffusive homeostasis. [NO 0 ] i for each neuron i can then be drawn from the distribution of steady-state intracellular concentrations of NO for this network. This results in a broad and heavy-tailed distribution of homeostatic targets, as opposed to the single homeostatic target which is used in networks with diffusive homeostasis and the unmodified non-diffusive homeostasis.
A similar approach was adopted in the dynamic mean-field analysis, with each neuron assigned a target firing rate ϕ 0,i from the steady-state firing rate distribution of a network with α = 0.8.

Procedure for Fig 4
External input for each neuron i was μ i = 2.5 Hz (N = 2500). NO 0 was set as described previously, although with an input rate of 2.5 Hz. 2 groups of 250 excitatory neurons each were randomly chosen, independent of neuron position, and stimulated with μ green = 5 Hz and μ blue = 10 Hz, keeping the inputs to the remaining neurons unchanged. Firing rates plotted in Fig 4A-4C were smoothed with a uniform time window of 20 s. Persistence of input differences were calculated by measuring the length of time it took for the signal-to-noise ratio between the two groups receiving elevated inputs to fall below 0. The signal-to-noise ratio is defined as (μ 1 − μ 2 )/(σ 1 + σ 2 ), where μ i and σ i correspond to the mean and standard deviation of the firing rates of group i. Fig 5A-5D were generated using the same simulation setup as described previously. After the network has reached the homeostatic target firing rate, we froze homeostasis. Input rates to each neuron were then regenerated from the same original input distribution, such that m after i $ N ð10; 10 2 Þ Hz. Dm i ¼ m after i À m before i is the difference in input rate each neuron experiences upon this change, and Dn i ¼ n after i À n before i is the corresponding change in output rate for each neuron. The black lines in Fig 5A-5C are from least-squares linear regression, and the R 2 values given were derived from this fit. A similar approach was used in the dynamic mean-field analysis, while varying δ, the width of the input distribution.

Simultaneous diffusive and non-diffusive homeostasis
In networks with both diffusive and non-diffusive homeostasis active, Eq (6) is replaced with Eq (12), where τ diffusive and τ non-diffusive represent the timescales of diffusive and non-diffusive homeostasis respectively, and the dynamics of [NO] and [NO non-diffusive ] are as described in Eqs (5) and (7) respectively.
Fig 6A was generated using the same simulation parameters as in Fig 1, while varying τ non-diffusive from 1000 to 15000 ms and keeping τ diffusive fixed at 2500 ms. These networks were run for 350 s in order to ensure they reached a steady state. Likewise, Fig 6B was generated using the same simulation parameters as in Fig 5A-5D, while varying the ratio of timescales.

Time-varying input
In addition to the external input μ i previously described, the network was randomly separated into groups of 250 neurons. Each group j was stimulated with external Poisson input with a rate given by μ j,t * N(0, 25 2 ) Hz. These inputs were regenerated at each timestep t of length 1 s.
The time-varying input μ j,t was also presented during homeostasis. The dotted black line in Fig  7A shows the normalized μ j,t , while colored lines show the normalized rate deviation of a randomly chosen group j from the mean population firing rate.

Decoding stimulus orientation
Each excitatory neuron was randomly assigned a preferred orientation. After the network reached a steady state, homeostasis was frozen. For each trial, each neuron i with preferred orientation θ i was stimulated with external Poisson input at a rate given by μ b + N(θ s , σ s , θ i ), where μ b = 20Hz is the base input rate and N(θ s , σ s , θ i ) is the amplitude at θ i of a Gaussian tuning curve centered around the stimulus orientation θ s , with a width given by σ s = 90°and a peak amplitude of 2.5 Hz. The angle decoded using the population vector method is the angle of the vector sum of all neural responses.

Spike-timing-dependent plasticity
A spike-timing dependent plasticity rule, as described in [19], was implemented in each recurrent excitatory-excitatory synapse. Both potentiation and depression are additive in this rule, with no weight dependence. For each nearest neighbor pair of pre-and post-synaptic spikes separated by a time Δt, the synaptic weight is updated by a value Δw given by if Dt < 0: ( Above, τ + and τ − denote the timecourse over which potentiation and depression occur respectively, while A + and A − denote the relative strengths of potentiation and depression. Weights are bounded by 0 < w < g max , therefore leading to a bimodal weight distribution. Initial weights are generated from a Gaussian distribution, given by w 0 * N(7.5, 2.5 2 ) nS. Depression is set to be slightly stronger than potentiation (A + < A − ) in order to prevent runaway potentiation, thus ensuring that irregular firing is maintained within a reasonable range of rates. The external input in networks with STDP was given by μ i * N(10, 10 2 ) Hz, and J ext = 40 nS, C = 250. STDP and homeostasis were frozen after 2000 s in order to facilitate measurements of the steady state for Fig 8A and 8B, and to measure the network response linearities in  Fig 5A-5D. The R 2 values for networks with shuffled, static synaptic weights were obtained by taking the steady-state synaptic weight matrix from a network with STDP and without any homeostasis, randomly shuffling this final weight matrix, and using these weight matrices in new networks for each different case of homeostasis, without STDP active. Using shuffled weight matrices obtained from networks with other forms of homeostasis present did not qualitatively affect the results.

Model parameters
Unless explicitly defined, the parameters used throughout the paper are given in Table 1. For synthesis, diffusion, and decay of NO we have attempted to match data when available [40,42], although the dearth of experimental measurements does not permit for great precision [43,44]. Additionally, parameters were chosen such that the timescale of homeostasis is separated from that of firing rate fluctuations. This is a reasonable assumption, given that activity-dependent NO modulation likely acts within 10 minutes or slower [5], although NO diffusion occurs on the order of 10 seconds. τ HIP was chosen to be long enough so as to avoid oscillations but short enough so as to allow feasible large scale simulations. This is a common assumption in computational studies [10]. Larger simulations, up to N = 25000, were run with no discernible difference in results.
All numerical simulations were implemented using the Brian simulator, v1.4.1 [45], and the mean-field analysis was implemented using IPython Notebook [46]. The 2D diffusion equation was solved numerically using an explicit finite difference equation method, using the numpy python package [47]. Data analysis was performed with the numpy Python package and plotting with the matplotlib package and seaborn library [47,48]. Simulation code and IPython Notebooks which perform the data analysis and plotting are available at https://github.com/ yannaodh/sweeney-2015. These results qualitatively agree with those for the random networks used throughout the study. The simulations shown here were identical to those in Figs 1D and 5 of the main text, but the connection probability between neurons had a Gaussian shape. (A-C) ρ, the ratio of the connectivity range and the diffusive range, is varied across a wide range of values (ρ = 5.0, 1.0, 0.5). (D) Each point represents the R 2 value of a linear fit as in Fig 5A-5C, for one network. Spatial dependence in the connection probability between two neurons was introduced as follows:

Supporting Information
where d is the Euclidean distance between the neurons and s is a constant defining the connectivity range of the network. 2D positions on the torus are bounded such that x, y 2 (0, 1). Given a diffusive range of 0.1mm, values for s were therefore set as 0.05, 0.1, and 0.5. The ratio of s and the diffusive range was defined as ρ, which had values of 0.5, 1.0 and 5.0. (TIF)