Depolarization induced suppression of excitation and the emergence of ultra-slow rhythms in neural networks

Ultra-slow fluctuations (0.01–0.1 Hz) are a feature of intrinsic brain activity of as yet unclear origin. We propose a candidate mechanism based on retrograde endocannabinoid signaling in a synaptically coupled network of excitatory neurons. This is known to cause depolarization-induced suppression of excitation (DISE), which we model phenomenologically. We construct emergent network oscillations in a globally coupled network and show that for strong synaptic coupling DISE can lead to a synchronized population burst at the frequencies of resting brain rhythms.

Endogenous cannabinoids (CBs) represent a fundamentally new class of retrograde messenger [1], which are released postsynaptically and bind to presynaptic CB receptors.CB synthesis is stimulated when levels of calcium rise inside the neuron or when certain G-proteincoupled receptors are activated.One function of endogenous CBs is to regulate neurotransmitter release via activation of presynaptic CB receptors, allowing neurones to regulate, via feedback, their upstream neuronal inputs [2].This suppression of upstream presynaptic release of GABA or glutamate is termed depolarization-induced suppression of inhibition (DISI) or excitation (DISE) respectively [3,4].Cannabinoid receptors are ubiquitous within the brain and CB1, the most abundant CB receptor, can be found in different areas such as the hippocampus, neocortex, amygdala, basal ganglia and hypothalamus [1].They have already been implicated in the temporal coordination of cell assemblies and the modulation of certain brain rhythms [1,5].The mammalian cortex is known to show oscillatory electroencephalogram (EEG) activity in a wide range of frequencies from approximately 0.05 Hz to 500 Hz [6].Recently the ultra-slow (0.01-0.1 Hz) fluctuations in functional magnetic resonance imaging (fMRI) signals (which have corresponding EEG correlates), as well as spatial patterns of their coherence, have gained the attention of the neuroscience community.This is mainly due to their robust and reproducible manifestation across subjects and their disturbance in pathological states including depression and Alzheimer's disease [7].So far, the mechanisms underlying such ultra-slow cortical rhythms, as observed in fMRI or as a modulation of the fast oscillations in electrophysiological signals [8], remain unknown.As regards the former Ghosh et al. [9] have developed a model using anatomical connectivity data (linking brain modules) that incorporates cortico-cortical communication delays and noise.However, by focusing on a non-synaptic form of coupling this work is more relevant to issues pertain-ing to the spatial coherence of fMRI signals than it is to probing physiologically based slow emergent rhythms.In this article we introduce a synaptically coupled network with a phenomenological form of retrograde second messenger signaling that can support DISE.Importantly we uncover a mechanism for the emergence of nested fast and ultra-slow oscillations.When linked to other modules in a larger network the latter would be reflected as an ultra-slow component of the macroscopic network dynamics and could therefore underlie those seen in the resting brain state.
Our single neuron model of choice is the Morris-Lecar (ML) neuron model [10].This is a two dimensional conductance based model, often used as an idealized fastspiking pyramidal neuron, written in the form v = f (v, w) + I + s(t), ẇ = g(v, w). ( Here v plays the role of a voltage variable, w that of a gating variable, I is a fixed input and s(t) represents a time varying synaptic input.Rather than list the details of the functional forms for f (v, w) and g(v, w), which we take from [11] (with time measured in ms), it is more natural to show a plot of the phase-plane and nullclines as in Fig. 1 for s = 0.The model can support either one or three fixed points, dependent on the choice of I.
The bifurcation diagram inset in Fig. 1 shows that the largest of the fixed points undergoes a sub-critical Hopf instability.Beyond this bifurcation there is a window of parameter space where a fixed point is bistable with a periodic orbit.Indexing each neuron in the network with i = 1, . . ., N the synaptic drive to the i-th neuron is given by where T m j is the m-th firing time of the j-th neuron, v s the synaptic reversal potential and W ij the connection 0 0.2 0.4 -0.5 -0.strength between neurons i and j with a global conductance scaling g s .The function η(t) captures the shape of a conductance change in response to the arrival of an action potential.Here we choose an alpha function and write , where H is a Heaviside step function.The firing times are specified in terms of a threshold h according to We focus on the case of an excitatory globally coupled network and set W ij = 1/N and v s > 0 with respect to the resting state.To implement a phenomenological model of DISE (or DISI if v s < 0) we first introduce a spatio-temporal average level of depolarization v e (t): where K(t) = 0 for t < 0. Since CBs are free to diffuse through neural tissue they allow for volume signaling, which we will consider as a form of global feedback.Their production is directly linked to depolarization, which we shall take to be the average tissue level as defined by (4).Here we choose K(t) = λe −λt H(t), where λ −1 is an indirect measure of the long time-scale for cannabinoid dynamics, which is on the order of tens of seconds to minutes [4].As a minimal model of DISE we will imagine that if the global CB level is sufficiently high then all excitatory synapses are blocked.In this case the network becomes uncoupled in the sense that excitatory synaptic currents drop to zero.By noting that this is equivalent to the suppression of firing events in (2), we may implement this model of DISE by letting the firing threshold adjust in response to v e (t) according to The threshold v th e controls whether the level of CB is sufficient to trigger DISE.In essence the model ( 5) means that synaptic interaction is curtailed if the mean level of depolarization becomes too large.
To probe the effects of DISE on network dynamics we first focus on the most symmetric states expected to exist in a globally coupled system -namely synchronous and asynchronous solutions [12].In the synchronous state all neurons have identical T -periodic trajectories with firing times given by T m i = mT for all i.In this case the synaptic drive to every neuron takes the identical form s(t) = g s (v s − v(t))P (t), where P (t) = m∈Z η(t − mT ) can be calculated as with P (t) periodically extended outside [0, T ).Equation (1) may then be solved as a periodic boundary value problem (PBVP) for the periodic orbit (v(t), w(t)) = (v(t + T ), w(t + T )) with v(0) = v th .This describes the synchronous orbit given that the corresponding mean depolarization does not trigger the DISE mechanism.This PBVP can be solved numerically, say using xppaut [13].
Although solutions exist for small g s a weakly coupled oscillator description, using standard techniques reviewed in [14], can be used to establish that such solutions are unstable.For stronger coupling there is a further window of g s where solutions exist, and by comparison with direct numerical simulations (implemented in Python using the Brian module [15]) are found to be stable.For an asynchronous splay state the firing times are given by T m j = mT + jT /N .In the limit N → ∞ network averages may be replaced by time averages due to: for any T -periodic function F (t) = F (t + T ).Hence a splay state in which all neurons fire is given by v i (t) = v(t + iT /N ), where v(t) is a T -periodic solution of (1) with s(t) = g s (v s −v(t))P 0 and P 0 = T 0 P (t)dt/T = 1/T .Note that for the splay state v e (t) takes on the constant value v 0 = T 0 v(t)dt/T .For small g s the splay state is found to have a similar period to that of the synchronous solution, though once again a weak-coupling analysis shows that this solution is unstable.Interestingly a clustered solution can occur for a wide range of g s .To see this consider two clusters of neurons, one in a splay state and the other sitting at rest.This can be described using the differential-algebraic system where r is the fraction of firing neurons and (v(t + iT /M ), w(t + iT /M )) with M = N r, and (v, w) describe neurons in the splay and resting cluster respectively.In this case v e = rv 0 + (1 − r)v.For v e < v th e the parameter region of existence for such a solution is illustrated in the inset of Fig. 2, where a pair of splay states (with r = 1) only coexists with a rest state for rg s ∈ [L, H].
Here the splay state is annihilated in a saddle-node bifurcation at rg s = H, while at rg s = L the lower branch solution seizes to cross the firing threshold.As g s is increased it is possible that v e could grow until it reaches v th e and activate the DISE mechanism.The border in the (r, g s ) parameter plane where v e = v th e for a cluster state is shown in Fig. 2, and we see that it defines a critical curve marking the onset of DISE which we can write in the form g s = g c (r).The line r = H/g s shows that in the absence of DISE cluster states would exist for a greater area of parameter space.For g s < g c (r) direct numerical simulations do indeed show cluster states with properties in excellent agreement with the solution of (8) (with v(0) = v th ) up to small fluctuations.An example is shown in Fig. 3.For a given value of g s the fraction of neurons r in the firing state is a function of initial data, as expected.Importantly, after transients, the mean depolarization signal is flat (no oscillations) and the period of oscillation of a firing neuron is of the same order of mag- showing that the network has split into two clusters (one with a common periodic orbit shown in blue with a period T ∼ 6 and a rest state in purple).v th = 0.05 (green line), v th e = 0.02 (red line).
nitude as a single isolated neuron.In the region where g s > g c (r) and DISE precludes the existence of a cluster state we expect more exotic non-periodic network states to emerge.In illustration of this point we note that stable synchronous oscillations are possible with increasing g s and that the average depolarization for these rhythms is an increasing function of g s .Hence there is also a critical value of g s at which the DISE mechanism will also preclude the existence of a periodic synchronous state.In this parameter window a synchronous (or near synchronous) solution can lead to a strong level of average depolarization for which v e (t) > v th e .This activates the DISE mechanism, precluding further synaptic input and subsequently leading to a drop in network firing activity and hence a drop in v e (t).Once v e (t) drops sufficiently to cross v th e from above then excitatory synaptic currents can once again drive the network leading to an increase in v e (t) so that the process may repeat over.In this case the emergent time scale of the network rhythm is set by the duration of v e (t) above v th e .Even for a synchronous solution this will depend on initial data, so that network oscillations would not generically be periodic.To quantify the value of possible inter-spike intervals (ISIs) we focus on synchronous rhythms with (v(0), w(0)) = (v th , w 0 ) for some given w 0 and solve the BVP v e (0) = v th e = v e (∆) with s(t) = g s (v s − v(t))P (t).The growth of the ISI, ∆, as a function of g s is shown in Fig. 4, together with results from direct simulations.The numerical spread of ISIs for low g s can be ascribed to fast multi-spike bursts.With higher g s a single spike response is more common and the period of the network state is accurately predicted by the theory.Note that the spike times considered here are only those that contribute to synaptic currents, while the neurons do in fact spike on a fast time scale during the synaptically silent period.Hence, the network as a whole shows nested oscillations with a slow variation of synaptic currents superimposed on fast oscillations of the instantaneous average network voltage (see Fig. 5 bottom left).To understand how decreasing v th e can lead to a rapidly increasing ∆, as shown in the inset of Fig. 4, it is useful to develop the correspondence of the evolution of the network (fixed parameters) with that of a single neuron with varying background drive I. Referring to the inset of Fig. 1 the network can leave point A, corresponding to a synchronous firing state with average voltage v 2 , when v e (t) drops below v th e .The subsequent large increase in synaptic drive causes a transition to the right of the saddle-node of periodics, where firing is not possible, and so synaptic currents fall which causes the transition to point B. This unstable fixed point, with voltage v 1 , generates orbits which spiral outward for a time T 1 = T 1 (g s ) generating a signal with v e (t) > v th e (so that synaptic currents are suppressed).These transition to full blown nonlinear oscillations, with average voltage v 2 and ve (t) < 0, and complete the path to point A so that the process may repeat over.Making the convenient (and obviously not accurate) assumption that v 1,2 are constant then the BVP may be solved by hand for λ = 0 to give ∆ = T 1 (g s )(v 1 − v 2 )/(v th e − v 2 ), explaining the dependence of ∆ on v th e seen in Fig. 4. Importantly, without any parameter fine-tuning, we see the emergence of very large ISIs for large values of g s , which are largely independent of the network size.Moreover, in contrast to other network models of slow oscillations (< 1 Hz) [16] we do not require a mixture of excitation and inhibition, and as shown in the inset of ISIs on the order of tens of seconds.Thus DISE in the strong coupling regime is a candidate mechanism for the generation of ultra-slow rhythms.

FIG. 2 :
FIG.2:(Color online).Fraction of firing neurons r as a function of the synaptic coupling strength gs.The inset shows the bifurcation diagram (in maximum amplitude) for a splay state with ve < v th e as a function of the product rgs.Here the rest state v is less than the threshold v th and the two branches of splay solution coalesce in a saddle-node of periodics at rgs = H.In the main figure the unlabeled dotted curved line shows the solution branch for which ve = v th e = 0.02.

FIG. 3 :
FIG. 3: (Color online).A cluster state for N = 100, g = 0.03, v th e = 0.02, λ = 10 −5 .Top left: a plot of the average signal ve(t), showing that after transients the emergent state lies below the threshold to activate DISE (red line).Top right: A raster plot of spike times, illustrating the drop-out of some neurons and the emergence of a splay state with the fraction of firing neurons r = 0.38.Bottom left: The average network potential b v = P N i=1 vi/N oscillates around the predicted value rv0 + (1 − r)v (magenta line) for r = 0.38.Bottom right: Phase plane dynamics for the network (dropping transients)showing that the network has split into two clusters (one with a common periodic orbit shown in blue with a period T ∼ 6 and a rest state in purple).v th = 0.05 (green line), v th e = 0.02 (red line).
FIG.4:(Color online).The predicted synchronous population ISI (in ms) as a function of gs, for w0 = 0.121 (green line), fits the ISIs seen in direct numerical simulations with N = 100 (red dots).Other parameters as in Fig.3.The inset shows the increase in ISI with decreasing v th e for gs = 1.

Fig. 4 FIG. 5 :
FIG. 5: (Color online).A similar plot to Fig.3showing the emergence of slow synchronized firing patterns in the strong coupling regime with gs = 0.5.Other parameters as in Fig.3.Bottom right shows voltage traces of 5 neurons (arbitrary offset for better display).Any variability due to heterogeneous initial conditions does not affect the emergence of ultra-slow near-synchronous oscillations.