article: The standard architecture of neocortex is a network with excitation and inhibition in closely maintained balance. These networks respond fast and with high precision to their inputs and they allow selective amplification of patterned signals. The stability of such networks is known to depend on balancing the strengths of positive and negative feedback. We here show that a second condition is required for stability which depends on the relative strengths and time courses of fast (AMPA) and slow (NMDA) currents in the excitatory projections. This condition also determines the response time of the network. We show that networks which respond quickly to an input are necessarily close to an oscillatory instability which resonates in the delta range. This instability explains the existence of neocortical delta oscillations and the emergence of absence epilepsy. Although cortical delta oscillations are a network-level phenomenon, we show that in non-pathological networks, individual neurons receive sufficient information to keep the network in the fast-response regime without sliding into the instability.
  It is generally accepted now that the model of the brain operating as a feedforward system is incorrect. Instead, circuitry in cortex and other brain areas constitutes a finely balanced network of strongly interacting excitatory and inhibitory neuronal populations [1–3]. These networks can respond to their input with high temporal precision [4–6], selectively amplify patterned input signals [7], transmit multiple signals simultaneously between neural assemblies embedded in large networks [8], and maintain activity on a broad range of time constants [9] including those of short-term memory [10]. However, the existence of positive feedback in these networks requires careful maintenance of stability. Previous work [11–13] has shown that stability requires a balance between overall excitatory and inhibitory feedback, see Eq (5) below. Other studies have examined the impact of N-methyl-D-aspartate, (NMDA) on stability in the context of working memory [14,15] and its effect on homeostasis [16]. However, the impact of fast α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid, (AMPA) and slow NMDA glutamatergic currents on the relative speed of positive and negative feedback has not been previously examined. Here, we show that the fact that excitatory currents have two main components with vastly different dynamics requires an additional stability condition, the “temporal balance condition”, which describes the relative speed of positive and negative feedback, Eq (6) below. Even when the strengths of excitatory and inhibitory connections are perfectly balanced, violation of the temporal balance condition makes the network unstable. As we show, deviations from the temporal balance condition by a few percent are sufficient to move the system into a state where it oscillates in the delta range and, subsequently, becomes unstable. This instability may cause absence epilepsy (petit mal) seizures [17–19]. Conversely, as the network approaches this instability it is also able to respond more quickly to changes in input. In fact, we show that near this instability the addition of slow NMDA receptors can make the network respond more quickly than when recurrent feedback is modulated by AMPA receptors alone. To understand stability and response times in networks with AMPA and NMDA projections, we will study populations with a balance of excitatory and inhibitory input. To gain an intuitive understanding of such a network, we will begin with a simplified single population having both excitatory and inhibitory recurrent connections, Fig 1A, whose basic dynamics are those of a simple second-order system (dampened spring), Fig 1C. We then examine a network with both excitatory and inhibitory populations, Fig 2A, and describe how the long time constant NMDA receptors impact network oscillations and response times, Fig 3. Finally, we show that the described behavior is naturally produced in networks with short-term depression (STD), Fig 4, and that it occurs not only in mean-rate approximations but also in spiking networks, Figs 2C and 4D. Since, as we show, small changes in the AMPA/NMDA ratio can affect both a network’s stability and its response time, homeostatic mechanisms that maintain the balance of the time constants of EE and IE projections are required to keep the system stable and, at the same time, its dynamics in the physiological range. Current theories of homeostasis use as control parameter either the average firing rate of a neuron or the activity history of individual synapses, and corrections are implemented in terms of spike timing dependent plasticity mechanisms [20]. Neither of these provides the information needed for maintaining the AMPA/NMDA ratio within the stability range of the temporal balance condition. What is needed is, instead, information about the frequency response of the network. We show that individual neurons in a spiking network have sufficient information about the network frequency response to allow them to identify, and if necessary counteract, instability in the local network, Fig 5. Failures of homeostatic control may allow the system to transition through the oscillatory regime into the instability, leading to absence epilepsy.
  A network is considered to be balanced if its excitatory and inhibitory inputs cancel in the long term, even though transient changes can still have a significant impact. Fig 1A shows a network with a single population, P. This network has one excitatory and one inhibitory recurrent projection, each with a combination of fast and slow receptors, and both having exactly the same strength over the long term. Therefore, if the ratios of fast to slow connections are equal on both projections, excitatory and inhibitory feedback exactly cancel at all points in time. However, if the ratio of fast to slow inputs on the excitatory projection is increased compared to the inhibitory projection, small increases in the firing rate of the population will cause a transient excess of excitatory feedback followed by a smaller but longer period of inhibition, Fig 1B. On the other hand, when the ratio of fast to slow inputs on the inhibitory feedback is greater, a small increase in firing rate will cause a transient inhibitory feedback followed by a smaller but longer period of excitation. Since the overall strength of the connections are the same, the total recurrent excitation and inhibition still ends up balanced. These transient but balanced changes act mathematically like a first derivative which either slows down or speeds up the network, similar to the “derivative networks” that have been previously described [10]. For physical systems the first derivative often acts like a friction term, or viscosity in a fluid. The network in Fig 1A can be approximated by a spring in a viscous fluid with the equation, d 2 R (t) d t 2 + 2 ζ ω 0 d R (t) d t + ω 0 2 R (t) = I (t). (1) where R (t) is the firing rate of the population P, ζ is a damping coefficient, ω0 is the frequency at which the spring oscillates when ζ = 0, and I (t) is some time dependent input. A full derivation of this representation, showing how it is approximated from the balanced network described in Eqs (2) – (4) and shown in Fig 2A is given in the Supplement. In the spring approximation, changes in the relative timing of excitatory and inhibitory feedback change the damping coefficient, ζ, effectively acting like changes in the viscosity of a fluid surrounding the spring. The diagrams in Fig 1C depict four different types of dynamical responses that a spring can exhibit as a function of ζ. The first type of response is where the fluid has no impact at all. In this case the spring oscillates with a constant amplitude and at the system’s natural frequency ω0, about 3. 5 Hz for realistic cell and synaptic parameters in the spring approximation where ω 0 = 1 / τ e τ n m d a. We note that the units of ω0 in this formula are rad/s. We define τe as the membrane time constant and τnmda as the decay time of the NMDA receptors (see Supplement). This behavior, shown in the first column in Fig 1D, is called an undamped response. The increasing frequency of ω0 that can be observed as the full system is reduced to the harmonic oscillator, Fig 1D bottom row compared to the top row, is due to a reduction of the number of time constants in the system. For example, setting the value of τampa = 0 removes two terms in the denominator of the equation for ω0, compare the analytical solutions in Eqs (S5) and (S11), causing a commensurate increase in ω0. Although we do not have an analytical solution for the full network, the addition of τi and τgaba, where GABA is gamma-Aminobutyric acid, could explain the lower value of ω0 for the full network. For even smaller values, i. e. ζ < 0, the network is unstable and oscillations will increase in an unbounded manner. On the other hand, for increased viscosity a step input generates a transient oscillation which decays over time toward constant steady state activity, second column in Fig 1C and 1D. This is an underdamped system, 0 < ζ < 1. As the viscosity continues to increase it reaches a point where oscillations are no longer possible and the system responds with an exponential decay to the steady state, third column in Fig 1C and 1D. This system is called critically damped, ζ = 1. A critically damped response is the fastest possible non-oscillating response for that system. Any further increases in viscosity continue to slow the exponential approach of the network of its steady state, fourth column in Fig 1C and 1D. Such networks are called overdamped, ζ > 1. So far we have discussed the behavior of the idealized spring system, Eq (1), shown in the top row of Fig 1D. The middle and bottom rows show, respectively, the equivalent results for the reduced network depicted in Fig 1A, and the full network, Fig 2A. Clearly, the spring model approximates well the responses of both of these networks. The behavioral repertoire of all three systems comprises an unstable regime, undamped and damped oscillations (first and second column), fast responses to input close to critical damping (third column), and a slower asymptotic approach to steady state for larger damping (fourth column). For both network models (Figs 1A and 2A), these different behaviors are obtained by only changing the ratio of fast and slow inputs on the excitatory-to-excitatory (EE) or excitatory-to-inhibitory (IE) projections. In the next section we will consider in more detail how this activity comes about in the full network of Fig 2. We now consider a recurrent rate-based network with AMPA, NMDA and GABA synapses, shown schematically in Fig 2A and described by τedRedt=−Re+JeeampaSeeampa+JeenmdaSeenmda−JeigabaSeigaba+I (t) (2) τidRidt=−Ri+JieampaSieampa+JienmdaSienmda−JiigabaSiigaba (3) τ m n l d S m n l d t = - S m n l + R n (4) Variables Re and Ri represent the firing rates of the excitatory and inhibitory populations, with intrinsic time constants τe and τi. J m n l is the synaptic strength of the projection from population n to population m of synaptic type l which is GABA for the inhibitory and either AMPA or NMDA for the excitatory projections. S m n l is the synaptic activation level of the projection from population n to population m, with synaptic time constant τ m n l. Time varying input to the excitatory population is denoted by I (t). We use two primary concepts to describe the speed at which a network responds to inputs: network time constant and rise time. We use the term network time constant, τn, to represent the time constant associated with the dominant eigenvalue, λ, of the network where τn = −1/λ. However, in many cases dominant eigenvalue approximations are not appropriate. Therefore, in order to compare across all networks, we also use the term rise time which is the time it takes the network to go from 10% to 90% of its steady state value. When the dominant eigenvalue approximation is appropriate, consider a system described by a single decaying exponential with time constant τn, then the rise time is approximately equal to ln (9) × τn. Projections from the excitatory population have AMPA and NMDA components, each carrying part of the total synaptic strength. The stability conditions and time constants for the linear network can be derived from its eigenvalues, all of which need to have negative real parts for the system to be stable. Assuming that the synaptic strength J is large compared to the neuronal leak current and that all projections are O (J), we develop approximate conditions for both stability and the network time constant, τn by approximating the coefficients of the characteristic polynomial in the highest order of J. Requiring the coefficients of the characteristic polynomial to be all positive is a necessary condition for stability while the ratio of the first two coefficients determines τn. This leads to two conditions for the network which are derived in the supplementary material, J e i (J i e a m p a + J i e n m d a) - J i i (J e e a m p a + J e e n m d a) > 0 (5) J e i (J i e a m p a τ i e n m d a + J i e n m d a τ i e a m p a) - J i i (J e e a m p a τ e e n m d a + J e e n m d a τ e e a m p a) > 0 (6) For ease of notation, we omit the superscript gaba for the inhibitory connections in this equation. In both equations, the first term characterizes negative feedback and the second positive feedback. Eq (5) is the previously described condition for the balance between strengths of inhibitory and excitatory connections [10,21]. When the left-hand side (LHS) of this “balance condition” equation becomes negative, the network is unstable since recurrent negative feedback (the first term) is smaller than recurrent positive feedback (second term). If the LHS is positive, the network is stable and dominated by inhibition. When Eq (5) is fulfilled, the novel “temporal balance condition, ” Eq (6), describes instabilities due to the relative timing of negative and positive feedback. If negative feedback weighted by synaptic strength (first term) is too slow relative to the weighted positive feedback (second term), i. e. when the temporal balance condition becomes negative, it cannot balance the excitatory feedback, thus making the network unstable. This is true even if total negative feedback (first term in Eq (5) ) is stronger than the total positive feedback (second term in Eq (5) ), i. e. when the steady state network is dominated by inhibition. On the other hand, positive and increasing values of the LHS of Eq (6) lead to increases in τn which can vary over a large range. When the LHS of Eq (6) approaches O (J2), then τn increases up to the range of several seconds. The network is then similar to the negative derivative feedback network introduced as a model for working memory [10]. While that study defined the stability conditions in terms of a single time constant on both excitatory projections, we focus on the more biophysically realistic implementation in terms of AMPA and NMDA receptors. Since the proportion of AMPA and NMDA receptors on each projection controls both the value of τn and the stability of the network there is a direct trade-off between speed of input response and stability. When the network responds quickly it is close to an AMPA dominated instability. As the NMDA component increases, the network moves further from this instability and τn increases. In order to analyze the stability of this network a change of parameters is useful. Let w be the base synaptic strength, k the inhibitory to excitatory ratio, q the proportion of synaptic strength carried by NMDA receptors (on all excitatory projections), and Δq a relative shift in the proportion of synaptic strength through NMDA and AMPA receptors on the EE projections only. With q ∈ [0,1] and Δq ∈ [−q, 1 − q], we have J e e a m p a = (1 - q - Δ q) w, J e e n m d a = (q + Δ q) w, J i e a m p a = (1 - q) w, J i e n m d a = q w and J i i g a b a = J e i g a b a = k w. Clearly, k parametrizes the impact of inhibitory vs. excitatory synapses and q the relative strength of NMDA vs. AMPA. Varying Δq represents changes in the temporal balance, note that for this parameterization Eq (6) is zero when Δq = 0. We chose to modulate the relative difference in time constants between EE and IE projections by adding Δq to EE; subtracting Δq from IE yields the same network dynamics. The simulation in Fig 2B shows the response of a linear network with q = 0. 3, k = 1. 2 and w = 30 to a step input at t = 0. Instability occurs for relatively small deviations from temporal balance. Increasing the AMPA contribution by setting Δq = −0. 02 makes the network unstable. The network resonates with a low frequency (in the delta range) as it approaches instability. S1F and S1G Fig shows the resonant frequency due to the AMPA dominated instability as a function of w and k. The frequency is fairly robust, varying between 1. 4 Hz and 2. 8 Hz over a large range of parameter values. The network is, however, very sensitive to changes in Δq: increasing Δq by 0. 1 towards NMDA currents increases the rise time of the network to approximately 5 seconds. For even larger Δq (slightly below 0. 15) another instability occurs, to be discussed below. A similar set of behaviors occurs for networks of leaky integrate and fire (LIF) neurons. Fig 2C shows the response of a network of 3,200 excitatory and 800 inhibitory neurons to a step input at t = 0. In our simulations, each LIF neuron has a connection probability of p = 0. 2 to all other neurons. The average synaptic strength from each population onto a downstream neuron is w = 5. 0, and the relative strength of inhibition is k = 0. 65. Results are qualitatively similar to those of the rate model, Fig 2B. However, there is an increase in the (absolute) values of Δq at which the network approaches either one of the instabilities (compare figure legends). The overall behavior is the same as in the rate network, a transition from instability and low frequency (delta) oscillations for negative Δq to long time constant integration for large and positive Δq. Fig 2D illustrates the trade-off between a network’s stability and its rise time for a step input. For the network to respond quickly to an input, it must maintain small negative values of Δq close to the dashed black line representing a rise time of 100 ms. Any increase in Δq will cause a large concurrent increase in τn slowing down the network response. Conversely, decreasing Δq leads to the emergence of a bifurcation with poles in the transfer function which move away from the real axis (S1H Fig). The right edge of the orange portion of the line in Fig 2D represents this bifurcation. As Δq becomes more negative, the resulting oscillations in the delta range (1–4 Hz, blue lines in Fig 2B, 2C and 2E) continue to increase in frequency until the poles cross the imaginary axis (S1H Fig) and the network becomes unstable, left red square in Fig 2D. To maintain fast response times, the system has to stay close to this instability. Therefore, balanced networks with short response times sit on the edge of stability: relatively small uncompensated changes in AMPA strength on the EE projections yield either an unstable network or a much slower stimulus response, both highly undesirable in sensory cortex. For large Δq the network rise time increases beyond 5 seconds at which point gamma-range oscillations appear (Fig 2E, green). Subsequently, the poles of the transfer function cross the imaginary axis, S1H Fig, and the system becomes unstable, right red square in Fig 2D. Thus, our model predicts that increasing τn by increasing the proportion of NMDA over AMPA receptors results in gamma oscillations and eventual instability. This type of gamma oscillations has been described previously [22]. The rise times spanned in Fig 2D extend all the way from physiologically realistic responses to sensory stimuli (tens of ms) to persistent activity suitable for working memory (seconds). These large shifts in the response dynamics of the network occur for a change in Δq less than 0. 2. The emergence of delta oscillations in balanced networks is consistent with the Stargazer model of absence epilepsy. Stargazer mice lack the stargazin protein which is expressed in inhibitory interneurons and is involved in AMPA receptor trafficking [17–19]. Animals lacking this protein are prone to seizures with increased EEG power in the delta and low theta range. The cause of these seizures has previously been ascribed to a reduction in the strength of inhibitory feedback from the loss of AMPA receptors on inhibitory interneurons [19,23] but this does not explain the occurrence of delta oscillations. While decreased inhibitory feedback may also play a role, our model suggests that a relative slowing of the IE projection through a reduction in the proportion of AMPA receptors breaks the temporal balance condition. Since, as discussed, stability depends on the relative timing of the EE and IE projections rather than their absolute values, a reduction of AMPA in the IE projection is equivalent to an increase in AMPA on the EE projection and will also induce a delta oscillatory instability. This view is strongly supported by recent work showing a compensatory mechanism which increases the strength of NMDA currents in stargazin deficient mice [18]. The addition of NMDA receptors compensates for the disinhibitory effect caused by the loss of AMPA receptors, allowing the system to fulfill the balance condition, Eq (5). However, the increase of NMDA also slows the negative feedback which moves the system even closer to the left instability in Fig 2B, thereby violating the temporal balance condition, Eq (6). This explains the occurrence of delta oscillations consistently observed in the stargazer model of absence epilepsy. This hypothesis can be tested directly by optogenetic interneuron and pyramidal neuron activation controlled by closed-loop, real-time recordings of pyramidal cell activity, a technique that is becoming well-established [24–27]. The prediction is that relative excess of slow (NMDA-like) optogenetically generated inward currents into interneurons over fast optogenetically generated (AMPA-like) currents into pyramidal neurons will result in delta oscillations and, when the imbalance is increased further, to seizure-like pathologies. Importantly, the model predicts that this occurs even though the balance of excitatory to inhibitory strengths is maintained. The dynamics observed in Figs 1 and 2 depend on the addition of NMDA receptors to the excitatory connections. A number of previous studies have considered response times in similar networks but without NMDA. They showed that in spiking networks response times on the order of or faster than the membrane time constant were possible and did not observe a delta oscillatory instability [4,5]. In Fig 3 we show how the dynamics in Fig 2 depend upon the amount of NMDA on the excitatory projections, q, and the time constant of decay for the NMDA receptors, τnmda. We consider a broad range of values for τnmda as has been previously observed [28–31] as well as values of q ranging from 0 as in previous computational work that did not consider NMDA receptors to NMDA/AMPA ratios that have been observed in cortex [31,32]. Qualitative characteristics of the network dynamics are maintained for almost the entire range of both q and τnmda. Fig 3A shows that the delta oscillatory instability is maintained for all parameter values except for a small region near q = 0. For example, when τnmda = 100 ms the network has no delta oscillatory instability for q ∈ [0,0. 01]. This is consistent with the stability observed in previous work but also shows that even small amounts of NMDA require consideration of the temporal balance condition. The oscillatory instability is also still in the delta range across the whole parameter regime, Fig 3B. Fig 3C shows the rise time for each parameter set when Δq is chosen such that the network is critically damped which we define as the point just before the poles involved in the delta oscillations separate from the real line. When the recurrent connections in the network have a higher percentage of AMPA receptors (smaller q) then the network generally responds more quickly to changes in input as would be expected from a network with faster synaptic responses. However, for some non-zero values of q the network actually responds more quickly than the network with no NMDA receptors, i. e. when q = 0. In these cases, each recurrent projection is slower then one with only AMPA receptors yet the negative value of Δq at critical damping drives the system causing an overall reduction in rise time. Interestingly, this driving of the network for negative Δq produces a counterintuitive interaction between rise time and τnmda: longer decay times for the NMDA receptors produce a faster rise time in the critically damped network. In fact, the fastest rise time we observe in our simulations of linear networks, 52. 5 ms, occurs when q = 0. 004, Δq = −0. 003 and τnmda = 400 ms. Rise times on this order have been previously observed in visual cortex [33]. Faster rise time for all our networks are possible if shorter membrane time constants and faster τampa and τgaba are used. Observation of an inverse relationship between rise time and τnmda in cortex would be an important experimental confirmation of our model. Although fast responses are possible for a subset of values of Δq, most of the stable values of Δq sit in the overdamped, large rise time regime. This regime shows effectively linear increases in rise time as a function of Δq as can be seen in the portion of the blue line in Fig 2D where Δq > 0. In Fig 3D we show how the slope of this line depends upon q and τnmda. Contrary to the critically damped network, in the overdamped network larger τnmda produces longer rise times for each positive increment in Δq as would be predicted from slower synapses. However, the value of q has almost no impact. The fact that the value of q does not change the relationship between rise time and Δq implies that even when q is zero the addition of NMDA receptors on EE connections can cause increases in the rise time of the network. Therefore, even in networks which start with no NMDA receptors, small additions of such receptors could still slow the response of the network significantly, making our analysis important even in such corner cases. In the previous section we showed that the stability, rise time and oscillatory activity of a balanced network depend on the relative strengths of AMPA and NMDA receptors that form the EE and IE projections. Our analysis assumed that the synaptic strengths are static. However, most synapses experience changes in their effective strength through short term plasticity which, on excitatory projections, is dominated by STD [34,35]. Since STD has a range of observed values, synapses with different AMPA/NMDA ratios may have STD with different strengths. Such a combined distribution of STD strength and AMPA/NMDA ratio will cause the effective value of q on each projection and therefore Δq to change over time. We therefore study the balanced network model with non-uniform STD on excitatory projections, as described in Materials and Methods, Eqs (9) – (12). We consider the evolution of the synaptic strength on excitatory projections due to STD by parameterizing synaptic dynamics around the base usage rate u = 0. 2, defining uampa = u − Δu and unmda = u + Δu. We also use a recovery time constant of τr = 500 ms. The base STD parameters are taken from within distributions found in rat visual cortex (we are not aware of equivalent data in primates) [36]. The parameter Δu models the same effect as Δq but in a time dependent manner. As with Δq, Δu = 0 maintains the temporal balance condition while positive values of Δu increase network damping and negative values yield an underdamped network. Fig 4A shows the behavior of the rate network for three different values of Δu. The parameters of the network meet both the balance and the temporal balance conditions at t = 0, and as the network evolves it exhibits dynamics similar to the linear network (without STD) studied above. Negative Δu causes the temporal balance condition to move towards the AMPA dominated unstable regime and the network begins to oscillate in the delta range. Positive Δu yields an effective increase in the proportion of NMDA currents resulting in a greater network time constant, in which case the network moves further from the AMPA dominated instability. Parameters of STD vary over a broad range in cortex [36]. Within the range of values of u and τr observed in that study, the qualitative dynamics of the network are the same as reported here, viz delta oscillations followed by instability for negative Δu, and long time constant responses for positive Δu, except for some deviations occurring for very strong STD. When u > 0. 35 and τr > 1. 5 s, the STD associated with a large initial spike dampens the subsequent response, thereby reducing the initial oscillatory activity seen in Fig 4A. However, the slow response for positive values of Δu remains the same for strong STD. In addition, increasing values of u and τr require larger values of Δu in order to produce similar dynamics due to weaker recurrent connections caused by the increasing strength of STD. Since the strength of the connections between populations are constantly changing due to STD, the networks in Fig 4A have different response characteristics as a function of time. In order to visualize how STD impacts these networks we plotted the rise time and parameters of linear networks with the same instantaneous synaptic strengths as the nonlinear networks at each time point, Fig 4B and 4C. We define instantaneous synaptic strength as the values of J in Eqs (10) – (12). Fig 4B shows how the rise time of the linear networks change over time while Fig 4C shows how the decreasing synaptic strengths due to STD project onto two parameters of the linear network, the strength of the EE and the IE projections We (see figure caption) and Δq. The regions of instability and the rise times plotted in Fig 4C are computed from the linear system using the synaptic strengths as defined in Materials and Methods. Although the network with STD is nonlinear, the approximate linear systems show qualitatively similar dynamics. When negative Δu yields faster rise time the linear system at each time point also shows faster rise times, Fig 4B. The same is true for slow responses with positive values of Δu. Additionally, as the linear networks cross into the delta oscillatory instability the nonlinear networks begin to oscillate as well, Fig 4A and 4C. This implies that the impact of STD on the network dynamics can be broadly understood using the temporal balance condition. Results transfer to spiking neurons. Fig 4D shows how networks of LIF neurons respond to similar STD parameters as the rate network in Fig 4A. As in the network without STD, larger changes in the underlying parameters are required to produce similar dynamics. We have shown that differences in the distribution of STD across EE and IE projections in a balanced network can produce significant changes in both network stability and response time (τn). In addition to maintaining the base distributions of NMDA/AMPA receptors and synaptic strength, the network must have homeostatic mechanisms adjusting the joint distribution of parameters controlling STD and NMDA/AMPA ratios. Concurrent modulation of synaptic strength, AMPA/NMDA ratio and short term plasticity parameters has been shown to occur in cortical networks [37,38]. Cortical networks have also been shown to carefully control the NMDA/AMPA ratio both at individual synapses and across cortical areas [31,32]. In the next section we address the question of whether the necessary information is available to the network to implement such mechanisms. We have shown that many neocortical networks operate close to an oscillatory instability associated with the relative ratios of AMPA to NMDA currents in the EE and IE projections. The ratio is determined by distributions of synaptic variables across the population but individual neurons cannot have direct access to information about the distributions. However, since each neuron in a local network receives both a sampling of external inputs and of the recurrent output of the population, it has implicit access to the network state through the frequency response, defined as its output as a function of its input. As seen from an individual neuron, the network input is simply the external input, and the network output is the recurrent input to the neuron from the network. Given a large enough sample of these inputs and outputs, an individual neuron can thus obtain an approximation of the frequency response of the network. Since increasing delta oscillations indicate that the network is approaching instability, analysis of the frequency response allows corrective action to be taken to maintain network stability. We tested whether individual neurons in a network of randomly connected LIF neurons can detect local network oscillations as Δq is altered. We used a constant Poisson input to drive the network. An approximation of the frequency response was then obtained for each excitatory neuron using the discrete Fourier transform of its recurrent glutamatergic inputs, i. e. the network outputs (note that at this point, our interest is only in the availability of the information, not detailed biophysical mechanisms). The magnitude of the average frequency response across ten independent runs shows a clear peak in the delta range as the network approaches instability (Fig 5A). The mean area under the transfer function in the delta range [0. 5Hz, 5. 5Hz] shows a clear increase as Δq becomes more negative (Fig 5C). A homeostatic mechanism based on this information could thus cause neurons to shift Δq towards stability. Similar responses, although in a slightly higher frequency range, were observed when the instabilities were due to changing synaptic strengths caused by STD (Fig 5B and 5C), allowing regulation of Δu.
  Networks with balanced excitation and inhibition are ubiquitous in cortex and other brain areas. It is well-known that the relative strength of excitation and inhibition in these networks needs to be controlled to keep them stable and to maintain their functionality. We show that in addition to this constraint on overall synaptic strength, a second condition is required to maintain network stability. This “temporal balance condition” specifies the relative weight of (fast) AMPA and (slow) NMDA receptors in the balance of positive and negative feedback. Temporal balance provides a highly sensitive parameter for setting the response time of the network which can be adapted over a range from tens of milliseconds to seconds by small adjustments of the AMDA/NMDA ratio. The range of stability is bracketed by two different oscillatory singularities. Networks approaching the NMDA-dominated instability express gamma-range oscillations, of a type described previously. A novel finding is that AMPA-dominated networks close to instability oscillate in the delta range, a possible source of cortical delta waves and a potential cause of absence epilepsy. Finally, we show that individual neurons in the network have access to information that allows them to homeostatically tune the set point for temporal balance to the optimal range. While previous work has shown that neurons in the asynchronous state could be driven to respond more quickly than their membrane time constant we show that different time constants on excitatory connections may also be able to speed up the network response [4,5]. In Fig 3C the network responds more quickly to input when the excitatory connections have a small amount of long time constant NMDA receptors than if the excitatory connections only have AMPA receptors. As we note in the main text, this makes the counterintuitive prediction that longer NMDA decay times may cause cortical networks to respond more quickly. Although a connection between fast network rise times and low frequency oscillations seems counterintuitive, recent work has made a connection between fast reaction times and delta oscillatory phase entrainment [45,46]. The connection between delta oscillations and fast reaction times has been explained as increased excitability for neuronal populations when the phase of delta oscillations and stimulus onset are appropriately aligned [47]. Our work indicates that power in the delta range may also correlate with faster rise times in the underlying cortical networks. The connection between delta oscillations and rise time also implies an interesting trade-off. A network which responds most quickly to a sensory input may also have long periods of relative quiescence during the trough of each delta oscillation. Therefore, such a trade-off may only be useful when the timing of the sensory input is predictable [48]. In future work it would be interesting to examine the interaction between excitability, rise time and delta phase as a function of AMPA/NMDA ratios. While the temporal balance condition is important in maintaining the stability of the underlying network it also allows for significant changes in the rise time. Previously it has been shown that, unlike positive feedback networks, the network in Fig 2A can change its gain without significantly altering its rise time [7]. Here we show that the converse is also true, such a network is able to alter its rise time without changing its gain, Fig 1C. Therefore, it may allow a network to change the temporal aspect of its neural code independently of its magnitude. The model we use to examine the impact of STD on the temporal balance condition assumes that STD impacts AMPA and NMDA receptors differently. This could be viewed as having two populations of synapses with either AMPA or NMDA synapses and different values for STD. Although physiological evidence does not support such a binary model there is evidence to support a spatial distribution of each type of receptor across a neuron and a broad range of values for STD parameters [36,50]. Additionally, LTP appears to induce transient changes in the AMPA/NMDA ratio which could be impacted by a distribution values of STD parameters [51]. How the AMPA/NMDA ratio is co- distributed with different STD parameters would determine the effect on the temporal balance condition. Future work on this subject should examine how different distributions of AMPA/NMDA ratios and STD parameters impact the rise time and stability of cortical networks. Maintenance of such a network requires homeostatic mechanisms which ensure that both balance conditions are met. Much of the current work on homeostatic mechanisms in cortex has focused on global synaptic scaling which has been shown to maintain the balance between excitatory and inhibitory projections [52]. These global mechanisms are exactly what would be required for maintenance of the balance condition in derivative feedback networks. Some experimental paradigms have also shown that the strength of AMPA and NMDA currents scale proportionately or are co-regulated during scaling [32,53,54]. Such mechanisms may be able to act quickly to regulate individual synapses [55]. Other experimental work has shown that STD and the ratio of AMPA to NMDA receptors are concurrently regulated during LTP through both pre-synaptic and post-synaptic processes [56,57]. These experiments show that at least in principle the mechanisms required for network level modulation of synaptic strength, STD and AMPA/NMDA ratio exist in neocortical synapses. Although we use the frequency response to examine the stability of our network, a full Fourier transform is not necessary. Knowledge of the change in amplitude across a set of relevant frequencies is sufficient, Fig 5C. Such a homeostatic mechanism would only require appropriate bandpass filtering of the incoming signal. Many studies have shown that frequency selective calcium signalling is important in intracellular homeostasis and control [39–41]. Since, the oscillatory activity in our model is in part driven by calcium permeable NMDA receptors similar frequency selective mechanisms could be used to drive the homeostatic response. Additionally, frequency selectivity is a relatively general characteristic of chemical systems implying that other signalling pathways are also possible [42]. The scaling of AMPA/NMDA ratios implied by our mechanism is determined by the parameters of the network as a whole rather than of the individual neuron. Therefore, the homeostatic mechanism should operate on the synapses associated with independent subnetworks to which the neuron is connected rather then scaling across all synaptic connections for the neuron. For example, local recurrent connections within a cortical column could be one appropriate subnetwork. If that network shifts towards instability then only synapses connected to other neurons in that network should be impacted. If an individual neuron is connected to multiple networks then the stability of each network should be approached separately. Although an examination of homeostatic mechanisms as a function of network connectivity is not available, synapse specific homeostatic scaling has been observed in previous work [43]. Additionally, target-specific short-term plasticity has been shown to exist in cortex implying that information about network level connectivity may be available to individual neurons [44]. It is generally assumed that synaptic strengths form an essential part of long-term memory. It therefore is natural to ask how the homeostatic mechanisms described here interact with the potential storage of memory contents. If, as is usually assumed, the information contained in each synapse is primarily related to its steady-state strength, then our mechanism should have minimal impact on memory contents and coding efficacy. The proposed frequency-based homeostatic mechanism while altering the AMPA/NMDA ratio does not alter the steady state value of each synapse for a given input. It only changes the transient activity of the synapse on the time scale of the NMDA receptors. Therefore, learning rules, read-out mechanisms etc. , as long as they are defined in terms of steady-state synaptic strengths, will be unaffected by the homeostatic adjustments. In this work, our model of NMDA does not include voltage dependence. A rate based approximation to voltage dependence has been developed and would provide a more physiologically plausible model [49]. However, when there are different amounts of NMDA on the EE and IE connections, the rate dependence of NMDA receptors tends to have a greater impact on the balance between excitation and inhibition than the relative timing of the feedback. Although maintaining such a balance must be an ubiquitous requirement in cortex and is not unique to our model, studying the excitatory/inhibitory balance condition was not the focus of our work. The implications of this additional nonlinearity would be an interesting topic for a future study.
  The reduced rate-based network depicted in Fig 1A is defined by five ordinary differential equations, a simplified version of the full network defined by Eqs (2) – (4): τedRdt=−R+w ( (1−q−Δq) S+ampa+ (q+Δq) S+nmda) −w ( (1−q) S−ampa+qS−nmda) +I (t) (7) τ m l d S m l d t = - S m l + R (8) where R represents the firing rate of the population, with intrinsic time constants τe. There are two recurrent projections, one excitatory and one inhibitory, each with a total synaptic weight w. S m l represents the synaptic activation of these two projections. m is the projection type, + for excitatory and − for inhibitory, and l is the synapse type, either AMPA or NMDA. In the diagram in Fig 1A, fast and slow are respectively equivalent to AMPA and NMDA in this formulation of the model. The parameter q is the proportion of each projection that is carried by the NMDA receptors. Δq represents a change in the ratio of AMPA and NMDA receptors on the excitatory projection relative to the inhibitory projection. Time varying input to the neuronal population is denoted by I (t). In the simulations in Fig 1C, middle row, we use w = 30, q = 0. 3, τe = 20 ms, τampa = 5 ms and τnmda = 100 ms. For the plots from left to right Δq is −0. 0425, −0. 0340, −0. 0095 and 0. 125. The value of Δq for the undamped oscillator was computed analytically from the characteristic polynomial, see Supplement. Parameters for critical damping were found by looking for the emergence of two complex roots near the point of instability. The values for the underdamped and overdamped systems were chosen to make the dynamics similar to the simulations from the spring approximation. The dynamics of the network are described by eight ordinary differential equations, Eqs (2) – (4). We consider time varying input to the excitatory population I (t) given by a step input of amplitude 5 Hz for the linear network and 10 Hz for the network with STD filtered by a 100ms Gaussian kernel. The time constants are: τe = 20 ms, τi = 10 ms, τampa = 5 ms, τnmda = 100 ms and τgaba = 10 ms. The values of the firing rate variables, Re and Ri, are constrained to always be positive except in Fig 1. The synaptic strengths are redefined in terms of w, k, q and Δq (see section on Rise Time and Stability in the Linear Network) to reduce the dimensionality of the parameter space and to automatically enforce the balance conditions. For all networks without STD we set q = 0. 3. In the simulations in Fig 1C bottom row we use w = 30, k = 1. 5 and from left to right Δq is −0. 0226, −0. 0196, −0. 0068 and 0. 05. The value of Δq for the undamped oscillator was found by searching for a root with zero real part as a function of Δq. Values for critical, overdamped and underdamped systems were found as for the reduced networks, see section on the Reduced Rate-Based Network. We use the STD mechanism described in ref. [58] for excitatory synapses (there is no STD on inhibitory projections): d x l d t = 1 - x l τ r l - u l x l R (t) (9) where xl (t) ∈ [0,1] represents the depressed synaptic efficacy due to STD, τ r l is the recovery rate of the synapse. The instantaneous rate of the input to the synapse is R (t) and ul is the usage rate which is a constant in our formulation. l is either AMPA or NMDA. Therefore, we have one Eq (9) for xl = xampa and another for xl = xnmda. We redefine the synaptic strengths of excitatory connections in Eqs (2) and (3) to include STD, J e e a m p a = 1 2 x a m p a w (10) J e e n m d a = 1 2 x n m d a w (11) J i e a m p a = J i e n m d a = x a m p a + x n m d a 4 w (12) Evolution of xampa and xnmda is determined by Eq (9) with parameters uampa, τ r a m p a, unmda and τ r n m d a. We define Δu = (unmda − uampa) /2. Since we are predominantly interested in the effect of STD on the temporal balance condition, we set STD on the IE projections to exactly balance the EE projections without changing the effective time constant. The balance condition is always met if q = 0. 5 and STD on the IE projections is the average of xampa and xnmda, Eqs (10) – (12). This also ensures that q remains constant and changes in the parameters of xampa and xnmda only change Δq. All our LIF networks have Ne = 3,200 excitatory and Ni = 800 inhibitory neurons. Each neuron is represented by the standard LIF equation and synapses have exponentially decaying activation, τ d V m d t = - (V m - E l) + ∑ n, l J m n l S m n l + I (t) + η (t) (13) τ m n l d S m n l d t = - S m n l + ∑ n, k δ (t - t m n α) (14) where τ is the membrane time constant and takes the values τe for excitatory and τi for inhibitory neurons as before. El = −60 mV is the reversal potential, Vm is the membrane voltage for neuron m with firing threshold −40 mV and a reset potential after firing of −52 mV. J m n l is the strength of the synapse from neuron n to neuron m on receptor type l, taking the value w/ (Nep) for excitatory synapses and w/ (Nip) for inhibitory synapses. The probability of making a connection from any neuron to any other is p = 0. 2 and the proportion of AMPA and NMDA receptors is redefined by q and Δq as in the rate model. S m n l is the synaptic activation, t m n α is the time of the α-th incoming action potential at synapse mn and τ m n l is the decay time for type l of synapse where τgaba = 10 ms and τnmda = 100 ms as before. τampa was increased to 10 ms to improve the stability of the network. For the simulations in Fig 2, each neuron receives Poisson-distributed background input η (t) from 1,000 neurons with event strength 0. 2 mV. The background input to excitatory neurons has firing rate 1. 05 Hz and to inhibitory neurons 1. 0 Hz; these values were chosen to maintain the same baseline firing rate for both populations. In addition, at t = 0 s the voltage of the excitatory population receives a step input of I (t) = 3 mV filtered by a 100 ms Gaussian kernel. All spiking simulations were run using the Brian spiking neural network simulator [59]. For LIF neurons with STD, R (t) in Eq (9) is replaced by ∑ n, k δ (t - t m n k). The p = 0. 2 connection probability is evenly split between an AMPA and an NMDA projection, effectively giving q = 0. 5. The usage rates, uampa and unmda on the EE projections differ by Δu as in the rate model. On the IE projections, the sign of Δu is reversed such that, Δuie = −Δu. For the simulations in Fig 4, neurons receive the same background stimulation as for the LIF network without STD. At t = 0 the voltage of the excitatory population receives a step input of I (t) = 6 mV filtered by a 100 ms Gaussian kernel ending at t = 1 s. The frequency response in the LIF network with and without STD is computed from a 4 second window of spiking activity in both networks. Each neuron receives the same background activity as in previous simulations. The network without STD receives a 2 mV step input and runs for 4 seconds prior to the window. The network with STD receives a 12 mV step input and runs for 1. 5 seconds prior to the window.
summary: Many networks in the brain are finely balanced, with equal contributions from excitation and inhibition. Deviations from this balance, if for instance the total amount of excitation exceeds that of inhibition, lead to potentially devastating instabilities. Unlike previous work we consider the interaction between fast and slow excitatory connections. We show that not only the amount of excitation needs to be controlled to achieve network stability but also the ratio of slow to fast excitation. Furthermore, optimally fast network performance requires that networks approach instability. However, networks very close to this instability develop oscillations in the delta range (1–4Hz) which potentially cause absence epilepsy. We show that a normal (non-pathological) network can auto-regulate its activity to avoid the instability.
keywords: medicine and health sciences
  neural networks
  nervous system
  electrophysiology
  neuroscience
  homeostatic mechanisms
  physiological processes
  viscosity
  systems science
  mathematics
  homeostasis
  materials science
  materials physics
  chemical properties
  physical chemistry
  epilepsy
  computer and information sciences
  system instability
  animal cells
  chemistry
  physics
  cellular neuroscience
  cell biology
  anatomy
  synapses
  frequency response
  physiology
  neurons
  neurology
  biology and life sciences
  cellular types
  physical sciences
  neurophysiology
title: The Edge of Stability: Response Times and Delta Oscillations in Balanced Networks