WWW.JNEUROSCI.ORG
-
The Journal of Neuroscience MBF Bioscience Autoneuron
 QUICK SEARCH:   [advanced]


     
-


HOME
  |  
SEARCH  |   ARCHIVE  |   SUBSCRIBE  |   CONTACT  |   HELP

This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Submit an eLetter
Right arrow Alert me when this article is cited
Right arrow Alert me when eLetters are posted
Right arrow Alert me if a correction is posted
Right arrow Citation Map
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (64)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Tabak, J.
Right arrow Articles by Rinzel, J.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Tabak, J.
Right arrow Articles by Rinzel, J.

 Previous Article

The Journal of Neuroscience, April 15, 2000, 20(8):3041-3056

Modeling of Spontaneous Activity in Developing Spinal Cord Using Activity-Dependent Depression in an Excitatory Network

Joël Tabak1, Walter Senn2, Michael J. O'Donovan1, and John Rinzel3

1 Laboratory of Neural Control, National Institute of Neurological Diseases and Stroke/National Institutes of Health, Bethesda, Maryland 20892, 2 Physiologisches Institut, Universität Bern, CH-3012 Bern, Switzerland, and 3 Center for Neural Science and Courant Institute of Mathematical Sciences, New York University, New York, New York 10003


    ABSTRACT
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
APPENDIX
REFERENCES

Spontaneous episodic activity is a general feature of developing neural networks. In the chick spinal cord, the activity comprises episodes of rhythmic discharge (duration 5-90 sec; cycle rate 0.1-2 Hz) that recur every 2-30 min. The activity does not depend on specialized connectivity or intrinsic bursting neurons and is generated by a network of functionally excitatory connections. Here, we develop an idealized, qualitative model of a homogeneous, excitatory recurrent network that could account for the multiple time-scale spontaneous activity in the embryonic chick spinal cord. We show that cycling can arise from the interplay between excitatory connectivity and fast synaptic depression. The slow episodic behavior is attributable to a slow activity-dependent network depression that is modeled either as a modulation of cellular excitability or as synaptic depression. Although the two descriptions share many features, the model with a slow synaptic depression accounts better for the experimental observations during blockade of excitatory synapses.

Key words: spontaneous activity; oscillations; depression; developing spinal network; recurrent excitation; rate model


    INTRODUCTION
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
APPENDIX
REFERENCES

Spontaneous activity generated by networks of synaptic connections is a characteristic, but poorly understood, feature of the developing nervous system. It has been detected in many different regions, including the retina, the hippocampus, and the spinal cord. In each region, bouts of activity occur episodically, separated by periods of quiescence. Although the mechanisms responsible for this activity are incompletely understood, it is known that developing networks are hyperexcitable, in part, because the classically inhibitory neurotransmitters, GABA and glycine, are depolarizing (Cherubini et al., 1991; Sernagor et al., 1995).

We have been studying this activity in the isolated lumbosacral spinal cord of the chick embryo. This preparation generates spontaneous episodes of activity in which spinal neurons are cyclically activated (see Fig. 1). This activity is not produced at a specific location, nor by a specific set of connexions (Ho and O'Donovan, 1993), and is still present in networks for which glutamatergic transmission or glycinergic/GABAergic transmission has been blocked (Chub and O'Donovan, 1998a). It has also been shown that episodes of activity transiently depress network excitability, which recovers in the quiescent phase; this depression is manifest as decreased synaptic strength (Fedirchuk et al., 1999) and hyperpolarization of individual cells. These observations have led to the proposal that the occurrence of episodes is a population phenomenon that is controlled by activity-dependent network depression (O'Donovan and Chub, 1997; O'Donovan, 1999). There is evidence that spontaneous bursts or episodes in other parts of the nervous system are also regulated by a form of activity-dependent depression. A slow "refractoriness" is thought to restrict wave generation and propagation in the newborn ferret retina (Feller et al., 1996, 1997; Butts et al., 1999). In hippocampal networks, NMDA receptor desensitization (Traub et al., 1994) or presynaptic depletion of glutamate pools (Staley et al., 1998) have been proposed to regulate the occurrence and/or duration of synchronous bursts.

Less is understood about the mechanism of cycling within an episode, but it has been proposed that this could also be regulated by activity-dependent depression with much faster kinetics than the depression controlling episodes [O'Donovan and Chub (1997); see also Streit (1993)]. Recently, it has been shown in modeling studies that a purely excitatory network with random connections susceptible to fast depression is capable of generating oscillations similar to the cycling observed during episodes (Senn et al., 1996; J. Streit and W. Senn, unpublished results).

Using these ideas, we develop a simple and general model of network activity. This model describes only the average firing rate of the population, rather than the membrane potential and spikes of the individual neurons. It assumes a purely excitatory recurrent network that is susceptible to both short- and long-term activity-dependent depression (see Fig. 2A). The model is described by only three coupled nonlinear differential equations. Therefore, the analysis can be done graphically, step by step, which facilitates an intuitive understanding of the network dynamics.

We show that this simple model can account for the occurrence of spontaneous episodes, the rhythmic cycling within an episode, and the developmental changes in the duration of episodes and interepisode intervals. Most importantly, it shows that the recovery of spontaneous activity in the presence of excitatory amino acid blockers (Barry and O'Donovan 1987; Chub and O'Donovan 1998a) is a property of networks with a slow form of synaptic depression.

Parts of this work have been presented previously in abstract form (Rinzel and O'Donovan, 1997) and conference proceedings (O'Donovan et al., 1998; Tabak et al., 1999).


    MATERIALS AND METHODS
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
APPENDIX
REFERENCES

Experimental

All experiments were performed on the isolated spinal cord of White Leghorn chicken embryos maintained in a forced-draught incubator. The spinal cord was dissected out under cooled (12-15°C) oxygenated (95% O2, 5% CO2) Tyrode's solution containing (in mM): 139 NaCl, 12 glucose, 17 NaHCO3, 2.9 KCl, 1 MgCl2, 3 CaCl2. The preparation was transferred to a recording chamber (at room temperature) and superfused with oxygenated Tyrode's solution that was later heated to 26-28°C. The concentration of KCl was raised to 5-8 mM during some of the recordings to accelerate the recovery of activity after glutamatergic blockade (Chub and O'Donovan, 1998a).

Neural activity was recorded from ventral roots or muscle nerves. Recordings were made using tight-fitting suction electrodes and amplified (DC 3 kHz) with high gain DC amplifiers (DAM 50 and IsoDAM, World Precision Instruments). Amplified signals were directly digitized through a PCI board (National Instruments) and/or recorded on tape (Neurodata) for further analysis. Recordings were digitized at a low sampling rate (5-20 Hz) because the events of interest occur on a slow time scale (Fig. 1).



View larger version (20K):
[in this window]
[in a new window]
 
Figure 1.   Spontaneous episodes of activity recorded from ventral roots of chick embryo spinal cord in vitro at embryonic day (E) 7.5 and 10. Recordings were made in DC mode and digitized at 10 Hz, showing the synaptic depolarization recorded electrotonically from motoneurons. Periods of activity (episodes) are clearly distinct from silent phases. An interval is defined as the time between the beginning of two consecutives episodes. Episodes are composed of cycles, or "network spikes." With development, both the interval and episode duration lengthen, and cycle frequency and the level of depolarization (tonic component) also increase. Also, on older embryos episodes have a leading tonic phase before cycling begins. For the examples shown here, interval and episode durations are, respectively, E 7.5: 440 ± 29 and 20 ± 1 sec; E 10: 1600 ± 140 and 68 ± 3.5 sec (mean ± SD).

Theoretical

Model equations. The two models are described and analyzed in depth in Results. Here we present them and the associated parameter values more succinctly. Both models have three dynamical variables, including a and d: a is the population activity or mean firing rate, and d is the fast depression variable, representing the fraction of synapses not affected by the fast synaptic depression. These two variables generate the cycles within an episode. The slow modulation that underlies onset and termination of episodes, and the long silent phases, is described by theta  or s. theta  represents the threshold for cell firing. In the theta -model, it increases slowly during an episode, eventually terminating the episode and recovering between episodes. In the s-model, s represents the fraction of synapses not affected by the slow synaptic depression, decreasing slowly during an episode, terminating it, and then recovering between episodes. In the s-model, the fraction of nondepressed synapses at any given time is therefore the product s·d.

The equations for the models are:

theta -model:
&tgr;<SUB><UP>a</UP></SUB><A><AC>a</AC><AC>˙</AC></A>=a<SUB>∞</SUB>(n·d·a)−a, (1)

&tgr;<SUB><UP>d</UP></SUB><A><AC>d</AC><AC>˙</AC></A>=d<SUB>∞</SUB>(a)−d, (2)

&tgr;<SUB>&thgr;</SUB><A><AC>&thgr;</AC><AC>˙</AC></A>=&thgr;<SUB>∞</SUB>(a)−&thgr;, (3)
s-model:
&tgr;<SUB><UP>a</UP></SUB><A><AC>a</AC><AC>˙</AC></A>=a<SUB>∞</SUB>(n·s·d·a)−a, (4)

&tgr;<SUB><UP>d</UP></SUB><A><AC>d</AC><AC>˙</AC></A>=d<SUB>∞</SUB>(a)−d, (5)

&tgr;<SUB><UP>s</UP></SUB><A><AC>s</AC><AC>˙</AC></A>=s<SUB>∞</SUB>(a)−s. (6)

The parameter n measures the connectivity in the network (see Results). All of the xinfinity functions are sigmoidal (e.g., see Fig. 2C for a plot of ainfinity ):
a<SUB>∞</SUB>(I<SUB><UP>ef f</UP></SUB>)=1<FENCE>(1+e<SUP><UP>−</UP>(<UP>I<SUB>ef f</SUB>−&thgr;</UP>)/<UP>k</UP><SUB><UP>a</UP></SUB></SUP>)</FENCE>

d<SUB>∞</SUB>(a)=1<FENCE>(1+e<SUP>(<UP>a−&thgr;</UP><SUB><UP>d</UP></SUB>)/<UP>k</UP><SUB><UP>d</UP></SUB></SUP>)</FENCE>

&thgr;<SUB>∞</SUB>(a)=1<FENCE>(1+e<SUP><UP>−</UP>(<UP>a−&thgr;</UP><SUB><UP>&thgr;</UP></SUB>)/<UP>k</UP><SUB><UP>&thgr;</UP></SUB></SUP>)</FENCE>

s<SUB>∞</SUB>(a)=1<FENCE>(1+e<SUP>(<UP>a−&thgr;</UP><SUB><UP>s</UP></SUB>)/<UP>k</UP><SUB><UP>s</UP></SUB></SUP>)</FENCE>.
The parameter theta x defines the argument's value for which half-activation occurs: xinfinity (theta x) = 1/2. The value of kx (positive) determines the sigmoid's steepness, smaller values being steeper. Note that the presence (respectively absence) of a minus sign in front of (a - theta x) indicates that the sigmoid is increasing (respectively decreasing). The variables d and s are high for low activity (low depression) and low for high activity (high depression), so dinfinity (a) and sinfinity (a) are decreasing functions of a.

Although ainfinity is the cellular input-output function, it should be considered as an average over the population of neurons. Assuming these neurons are not identical, ainfinity inherits a sigmoid shape; theta  can be interpreted as the mean threshold and ka as the threshold variance, assuming a unimodal distribution of thresholds (Wilson and Cowan, 1972). Because the connectivity coefficients are included inside the ainfinity function, a should be interpreted as mean population firing rate averaged over a brief period of time, not as instantaneous voltage (Pinto et al., 1996); ainfinity (n·s·d·a) is the instantaneous firing rate. Note that because ainfinity (0) is not equal to zero, we are implicitly assuming that there is some background activity or fraction of spontaneously active neurons in the network. However, the model does not depend on this feature. Indeed, if ainfinity (0) = 0, similar episodes of activity are produced, except that a stays virtually at 0 from just after an episode to just before the next one (our unpublished results).

Unless stated otherwise, we use the parameter values given in Table 1.


                              
View this table:
[in this window]
[in a new window]
 
Table 1.   Values of the parameters used for each model, unless mentioned otherwise in text or figure

Simulations. The model equations were implemented within XPPAUT (freely available software by G. B. Ermentrout, http://www.pitt.edu/~phase/), a general purpose interactive package for numerically solving and analyzing differential equations. XPPAUT includes a tool for calculation of bifurcation diagrams (AUTO). Simulations were performed using the Runge-Kutta integration method with a time step of 0.2 (dimensionless units). We checked that the results were unchanged when the time step was halved. Simulations were run on Unix PCs. A version of the software WinPP runs under Win95 or Win98.


    RESULTS
TOP
ABSTRACT
INTRODUCTION
MATERIALS AND METHODS
RESULTS
DISCUSSION
APPENDIX
REFERENCES

We have argued elsewhere (Ritter et al., 1999) that the spinal network responsible for spontaneous activity is composed of recurrently connected, functionally excitatory synaptic connections. Furthermore, because the expression of network activity does not appear to depend on the details of network architecture or membrane properties and because detailed biophysical information on cell and coupling properties is unavailable, we are justified in using a very idealized model. Throughout, we use mean field models that describe only the averaged population behavior, not the detailed electrical activity of individual cells. The models are also temporally coarse-grained so that action potentials are not represented per se, only firing rate.

Our graphical analysis of the models proceeds by considering reduced models, beginning with the fastest process first and expanding the model by stepwise introduction of the next slower variables. We begin with the single Equation 1 that governs the dynamics of a recurrent excitatory network, without the negative feedback of depression. This reduced model illustrates the classical behavior of bistability in a recurrent population, of either a low or high activity state of steady firing. Next, we introduce a fast activity-dependent synaptic depression, which can lead to network oscillations (cycling) in this two-variable system. We analyze this subsystem with phase plane methods, finding a different type of bistability. Finally, we add slow activity-dependent depression in either of two forms and show that these three-variable models can generate episodes similar to those observed experimentally in the chick embryo spinal cord and other networks.

A recurrent excitatory network can be bistable

In this section, we consider only the dynamics of the average activity in a random excitatory network; the depression variables are frozen. We use a graphical method to illustrate how the steady activity states of the one-variable system are determined. These steady states depend on parameters such as the network connectivity, and we show that for certain values of connectivity the network can have two stable states.

Here, and throughout, we let a be a measure of the network activity, for instance the population's mean firing rate, relative to the maximum possible rate. Its time behavior, according to classical descriptions (Wilson and Cowan, 1972), satisfies tau a a = ainfinity (Ieff- a, where a is the time derivative of a and Ieff is the effective input to cells in the network. Accordingly, the activity of the network reaches the value ainfinity (Ieff) with a time constant tau a. ainfinity is called the activation function and describes how Ieff is converted into firing rate (Fig. 2B). For small inputs, ainfinity is near 0; very few neurons are firing. For large inputs, ainfinity saturates at 1, with all cells in the network being fully active. There is a rather sharp transition between these two limits (if ka is small) at the value Ieff = theta  for which ainfinity (theta ) = 1/2. theta  corresponds to the average threshold for cell firing. The mathematical expression of ainfinity is given in Materials and Methods.



View larger version (16K):
[in this window]
[in a new window]
 
Figure 2.   A, Schematic diagram of recurrent network with random connections. The output from the network is fed back to the cells. The amount of feedback depends on network connectivity (n). Activity (a) can also be modulated at the synaptic level by a fast depression (d) and by a slow depression (s). The effective input to cells is therefore n·s·d·a. The resulting postsynaptic firing rate, at steady state, is given by ainfinity (a·d·s). This postsynaptic firing rate can also be modulated by the cells' threshold (theta ). Note, our mathematical formulation does not distinguish presynaptic and postsynaptic effects. B, Cellular input-output relationship, averaged over the population. Ieff is the effective input to the network (see Results). The cellular threshold parameter theta  can also be viewed as a slowly modulated variable. ka measures the width of the transition range around threshold; a small value of ka implies a steep activation function.

Because of the recurrent connectivity, the input Ieff is produced by the network's own activity, in addition to any external inputs. Therefore, in the absence of external inputs Ieff = n·a, where we introduce the parameter n as a measure of network connectivity (in arbitrary units). The greater the connectivity, the greater the recurrent input being fed back to the network for a given amount of activity (n is a gain parameter, measuring how the average cellular output is transmitted to other cells, and thus depends on many physiological parameters such as the average number of contacts from a neuron, average synaptic efficacy, etc.). Thus, the dynamics of this network is described by:
&tgr;<SUB><UP>a</UP></SUB><A><AC>a</AC><AC>˙</AC></A>=a<SUB>∞</SUB>(n·a)−a. (1′)
What is the behavior of the system defined by Equation 1'? After any transient input, a will evolve to a steady state, for which a = 0. Its value at steady state must then satisfy a = ainfinity (n·a). The solutions to this equation are represented graphically in Figure 3A as the intersections of the curves f(a) a and g(a) = ainfinity (n·a) for different values of n.



View larger version (15K):
[in this window]
[in a new window]
 
Figure 3.   Analysis of the single-variable activity model for a recurrent excitatory network. A, The steady states satisfy ainfinity (n·a- a = 0. They are determined by the intersection(s) between the curves f(a) = a (thick line) and g(a) = ainfinity (n·a) (smooth curves), shown here for different values of n. At these intersections ainfinity (n·a- a = 0. B, Time course of the activity for n = 0.5. At t = 20 a brief stimulus is applied but fails to bring the network above threshold (i). Stronger stimulation applied at t = 50 excites the network above threshold, so the high activity steady state is reached (ii). The network will stay at this high level unless a large perturbation brings it back below threshold. C, Plot of the steady-state values of a for all possible values of n between 0 and 1. The dashed portion of the curve indicates unstable steady states. In A and C, the steady states for n = 0.5 are represented by filled (stable) or open (unstable) circles. The arrows indicate that the activity, if initially below threshold (unstable steady state), will fall back to the low stable state (i, a < 0), and if it is initially above network threshold it will increase until it reaches the high stable state (ii, a > 0), as illustrated in B. theta  = 0.18.

For small values of n, there is only one steady state, at a low value of a, with the network barely active because of its low, functional connectivity. Even if the network receives a strong transient external input, it will return to this quiescent state because the connectivity is too low to sustain the activity. Similarly, for large connectivity, there is only one possible steady state (the curves do not actually intersect at the low level), for which the network is very active (a near 1). A more interesting case occurs for intermediate values of the connectivity parameter: there are three steady states, at low, intermediate, and high values of a. However, the intermediate steady state is unstable. That is, a small perturbation from it will grow, and the network activity will evolve toward one of the two other steady states. Therefore, the network has two possible stable steady states; it is bistable. If the network is in the low steady state, it can be switched to the high steady state by a transient stimulus that brings a above the intermediate steady state, which is the effective network threshold. If the stimulation does not increase the activity above threshold, the network will fall back to low activity after the stimulus. This network threshold is not equal to the cell threshold theta ; for example, if ainfinity is very steep, the network threshold approximately equals theta /n; it decreases with increasing connectivity. Note that once the network has switched from the low to high (or high to low) state, it will stay in the new state even though the perturbation has stopped. It will not return to the previous state, unless a new perturbation of sufficient magnitude occurs (if neural noise were incorporated explicitly in the model there might be spontaneous transitions between the two states). This type of bistable behavior, illustrated in Figure 3B (see figure legend), is characteristic of recurrent excitatory networks (Wilson and Cowan, 1972). In the next sections, we will find that bistability underlies our models' episodic rhythms; slow modulation leads to transitions between the attracting pseudostates of cycling during an episode and quiescence between episodes.

Figure 3C summarizes the above results by showing the activity levels at steady state (a = 0) over a range of n values. Again, it illustrates that for low connectivity/efficacy (n small) there is only one steady state at low activity. For strong connectivity (n near 1 and above) only the high activity state is present. For an intermediate range of n values both stable states coexist, along with an unstable state of intermediate activity (dashed line). The states are continuously connected as you travel along the S-shaped curve. However, as n is varied over the range, one encounters discontinuous transitions between the stable states, at the S-curve's knees. This figure, showing the different possible activity states for all values of a control parameter, is called a bifurcation diagram [see Strogatz (1994); a bifurcation is a point for which there is a qualitative change in the dynamics of the system, like at the knees]. It will be useful to distinguish the three "branches" of such S-curves: upper, lower, and middle. Here, upper and lower ones correspond to the stable steady states, and the middle one corresponds to the unstable steady state. The middle branch is a separatrix between the upper and lower branches: any state of the network above it (i.e., on the right side of the S-curve) will evolve toward the high steady state, and any state below it (i.e., on the left side of the diagram) will fall to the low steady state. It is easy to understand why: for any point on the S-curve, by definition, a = 0. If the connectivity is increased (i.e., the point shifted to the right) even slightly, this increases ainfinity (n·a), and therefore a > 0. Conversely, if n is decreased this leads to a < 0.

A recurrent excitatory network with fast activity-dependent synaptic depression can produce oscillations

In the above section, we studied the behavior of a recurrent excitatory network as a function of its (fixed) connectivity. Here, we will analyze network behavior when the effective connectivity is changing because of synaptic depression. We will analyze the resulting two-variable system graphically in the phase plane and show that it can generate oscillations (like the cycling during an episode) and bistability, depending on the average neuronal excitability.

As a preliminary step, reconsider Figure 3C and suppose that we sweep n slowly back and forth over a range that includes the bistable regime. This would lead to alternating phases of high and low activity, with sharp transitions between phases. If n were decreasing while the activity is high, the system's state point would track the upper branch, moving leftward, until we reached the S-curve's left knee. Then the state point would fall to the lower branch. The activity would drop precipitously. Now at low activity, if n were increasing, the system would track the lower branch, with activity rising very gradually, until the right knee was reached. Then the state point would jump back to the upper branch of high activity, completing a cycle, and the process would repeat.

We really want a description so that such cycling happens spontaneously, rather than being driven as above. With this goal we make the effective connectivity a dynamic variable by multiplying n by the variable d, the kinetics of which depends on activity, guaranteeing that d automatically decreases when activity is high and increases for low activity. In this formulation, d is defined as the fraction of nondepressed synapses from a representative cell to its targets (e.g., the fraction of releasable transmitter at those synapses). Its value ranges from 0 (synapses totally depressed) to 1 (all synapses fully available), so the effective connectivity (n·d) ranges from 0 to n (n, being an arbitrary measure of connectivity can be >1). In all of the following figures, we set n = 1 unless mentioned otherwise. Equation 1' thus becomes:
&tgr;<SUB><UP>a</UP></SUB><A><AC>a</AC><AC>˙</AC></A>=a<SUB>∞</SUB>(n·d·a)−a. ()
This idealized model does not distinguish between glutamatergic excitatory synapses or depolarizing gabaergic/glycinergic synapses.

We use a simple first-order kinetics for d, reminiscent of the gating kinetics in Hodgkin-Huxley-like models for voltage-gated ionic currents:
&tgr;<SUB><UP>d</UP></SUB><A><AC>d</AC><AC>˙</AC></A>=d<SUB>∞</SUB>(a)−d, ()
where dinfinity (a) is a decreasing sigmoid function (see Materials and Methods). When the system is at low activity, d tends to increase, whereas for high activity d decreases. d tends toward dinfinity (a) but lags behind changes in a with a time constant tau d. Hence, the system can potentially oscillate, but whether it does depends on parameter values. For example, if the negative feedback from depression is too fast, the oscillations are precluded. The condition (parameter relationship) for which oscillations emerge in the model is given in the Appendix.

We need to examine the dynamics of this a - d system over a range of theta  values, because in later sections theta  will become a slow modulatory variable. For theta  small, the network is highly excitable, and it operates at a steady high-activity level. If the network is held at a low activity level with synapses fully available and then released (Fig. 4A), it tends to its desired high level; in this case, it shows damped oscillations around the steady state because the depression is delayed (in the same way as the "delayed rectifier" potassium current does not turn on immediately after a voltage increase).



View larger version (18K):
[in this window]
[in a new window]
 
Figure 4.   Dynamical behavior of the a - d system for low firing threshold, theta  (0.15). A, Time course of activity a (solid curve) and depression d (dashed curve) showing that the system reaches a steady state after some damped oscillations. B, Phase-plane representation of the dynamics, showing the nullclines and trajectory of network state. The a- and d-nullclines are defined by a = 0 and d = 0. There is only one intersection, therefore the system has only one steady state. It is a stable steady state; as shown by the time course (A) or phase-plane trajectory. The damped oscillations evident in the time course (A) are seen here as spiralling toward steady state. n = 1.

The system's time course from Figure 4A can be viewed in another way, as a trajectory by plotting in the (ad) plane simultaneously the points defined by the values of a and d at consecutive and equally spaced times (Fig. 4B, dotted line). Although time is not represented explicitly on such a diagram, velocity along the trajectory can be deduced from the spacing between consecutive points. Here, on release from low activity the network quickly (large spacing between points) moves toward a state of high activity with the synapses depressed, slowing down (small spacing between points) as it approaches the steady state, in a spiraling fashion (the damped oscillation).

Using the a - d plane, we can analyze and predict qualitatively the system's dynamical behavior and the consequences of changing various parameters [for an introduction to phase plane analysis, see Rinzel and Ermentrout (1998) or Strogatz (1994)]. The power of this graphical/geometrical approach comes largely from understanding the two solid curves in Figure 4B. These curves are called the nullclines. The a-nullcline is defined as the relationship between a and d that is obtained by setting a = 0. Similarly, the d-nullcline comes from setting d = 0. Given these nullclines, we can map out the flow directions in the phase plane, according to the different areas delimited by the nullclines. For example, in the area on the right of the a-nullcline and above the d-nullcline, we know that a must be increasing and d decreasing (hence the arrow pointing up and left). Also, we know that a must be at a local minimum or maximum when the trajectory crosses the a-nullcline and so on. Where the nullclines intersect, neither a nor d is changing, and this corresponds to a steady state of the system.

The a-nullcline is defined implicitly according to a = ainfinity (n·d·a). For a physiological interpretation, consider d as a parameter so that points on the a-nullcline are the steady-state values of a at which decay of a balances the regenerative input from recurrent excitation in Equation 1. Of course, the result is similar to the S-curve in Figure 3C. The d-nullcline is simply the graph of dinfinity (a). The nullclines are drawn in Figure 4B, showing that there is only one possible steady state (the unique intersection of the two curves). In general there can be one or three (exceptionally two) intersections, as we will show below for some values of theta . The a-nullcline has the characteristic switchback form (S-shape) of autocatalysis in nonlinear excitable and oscillatory systems. Although the steady state is stable in this case, it could be made unstable by increasing tau d to slow the negative feedback process, in which case the activity would become oscillatory.

For sufficiently large values of theta , there is again only one steady state, but now at a low activity level: in Figure 5B, the unique intersection of the nullclines near the abscissa at d = 0.923 (vertical arrow). Note that the a-nullcline's low-activity branch is very close to the abscissa. This network model is excitable. If the system is stimulated briefly but adequately while at this quiescent, resting state, it exhibits a transient response of high activity and then returns to rest. Figure 5A shows the threshold effect by plotting a(t) and d(t) for two different initial conditions. In one case (gray curves), the activity starts just below threshold and directly returns to the steady state. In the other case (black curves), the activity is initially just above threshold and then increases in a regenerative way before depression returns it to the steady state. Note that this threshold for network excitability is approximately given by the value of the middle branch of the a-nullcline and should not be confused with theta , which is the average cellular threshold. This phase plane portrait gives us a graphical way to see the effect of increased threshold, by comparison with Figure 4B. The a-nullcline's right shift is understandable intuitively. For a given value of depression variable d, one expects with increasing theta  to lose the possibility of a high-activity state in the a-subsystem. With enough increase of theta  the high-activity state in the full a - d system is precluded.



View larger version (17K):
[in this window]
[in a new window]
 
Figure 5.   Dynamical behavior of the a - d system for high theta  (0.28). A, Time course for two different initial conditions. In one case, a(0) is above a threshold, and the network exhibits a transient high-activity phase and then falls to quiescence (black curves). In the other case, the network goes directly to "rest" (gray curves). a is indicated by the solid curves; dashed curves represent d. B, Phase-plane representation of the dynamics, showing trajectories for the two cases illustrated in A. There is only one steady state, defined by the intersection of the d-nullcline and the lower branch of the a-nullcline (for a sime  0)---the quiescent state. If the system is perturbed away from steady state by a sudden increase in activity, the threshold for generation of a transient cycle of high activity is approximately given by the middle branch of the a-nullcline (it is actually slightly above). n = 1.

Finally, for some intermediate values of theta , there can be three steady states, as shown on Figure 6, A and B. In both A and B, the lower steady state is stable, and the intermediate one (very close to the low state) is unstable. However, the stability of the high-level steady state depends on the particular values of the other parameters. An important parameter, which does not change the location of the steady states, is tau d. The high-level steady state is stable in Figure 6A, for a value of tau d = 1. Hence, the system is bistable. If the value of tau d is increased, the negative feedback caused by depression is more delayed, which can lead to oscillations around the steady state (now unstable), as shown in Figure 6B. The system is bistable again, but now we have a low-activity steady state coexisting with a high-activity oscillatory state. Depending on the initial conditions (or external inputs), the system can either reach the low steady state or the high activity state, either steady (Fig. 6A) or oscillatory (Fig. 6B,C).



View larger version (21K):
[in this window]
[in a new window]
 
Figure 6.   Dynamical behavior of a - d system for intermediate theta  (0.2). A, Phase-plane representation of the dynamics for tau d = 1. The high and low steady states are stable. Any trajectory eventually reaches one of them, depending on initial conditions. A perturbation away from the low steady state can switch the system to the high state, if it is above the middle branch (black trajectory). Otherwise, the system falls back to the lower state (gray trajectory). B, Phase-plane representation of the dynamics for tau d = 2. This higher value destabilizes the high steady state and leads to the appearance of a stable oscillation (closed orbit) around it. Again, a sufficient perturbation can switch the system from the low-activity steady state to oscillations. C, Time courses for the trajectories shown in B. Superthreshold perturbation leads to an oscillation (black curves), whereas after a subthreshold perturbation, activity returns directly to the low steady state (gray curves). The dashed curves represent d. n = 1.

We summarize this collection of possible behaviors in a bifurcation diagram (Fig. 7A) by plotting the a-values for the steady and oscillatory states when the parameter theta  is varied. The Z-shaped curve (black solid line and dashes) for steady states is analogous to the S-curve of Figure 3C; the orientation is left-right-reversed because increasing theta  reduces network excitability, like decreasing n. Starting at theta  low with the unique high-activity steady state and moving rightward along this curve, we see that this state loses stability at theta  > 0.181 (solid turns to dashed). At this point, a stable oscillation emerges, around the steady state. We plot the maximum and minimum values of a during a cycle (thick gray curves). At emergence, the oscillation's amplitude is very low, but then grows as theta  increases. This type of oscillation birth is called a Hopf bifurcation (Rinzel and Ermentrout, 1998). At a critical value of theta  (theta c = 0.207), the oscillation encounters the unstable state on the Z-curve's middle branch and disappears. That is, the minimum activity during a cycle is just equal to the excitation threshold, and the oscillation cannot be maintained if theta  is increased any farther. For larger theta  only the low-activity steady state is stable.



View larger version (13K):
[in this window]
[in a new window]
 
Figure 7.   A, Bifurcation diagram showing the possible behaviors of the a - d fast subsystem for different values of theta . For low values of theta , there is only one steady state, at high level, which becomes unstable for increasing theta . For intermediate values of theta  (around 0.2), the system is bistable: there is a stable high-activity oscillation and a stable, low steady state. With larger values of theta , the system can only be in a low steady state. Circles represent the cases shown in Figures 4B, 5B, and 6B. Thick black curves: stable steady states; dashed, thin curve: unstable steady states; thick gray curve: maximal and minimal values of a for the periodic oscillations. B, Period of the oscillations.

The system is bistable for theta  between the Z-curve's left knee and theta c. Over this range, it can be either at a low, steady activity level, or at a high, oscillatory level. These high-activity oscillations mimic the cycling observed during an episode of activity in the chick spinal cord. In the next sections we will see how slow autonomous dynamics for this modulatory variable theta  can sweep the a - d subsystem back and forth over this bistable regime to account for the cord's episodic rhythms. Thus, although Figure 7A does not show motion per se, it forms the skeleton over which episodic dynamic patterns can be understood. For example, as shown from Figure 6A,B, the value of theta  where the cycles emerge depends on parameters such as tau d. If tau d is decreased, the oscillations emerge for higher values of theta . Consequently, when an episode is started, the system first jumps toward a high stable steady state, then with increased theta  this high state becomes unstable and oscillations emerge. Thus, our model when the slow variable is included may exhibit, depending on parameter values, episodes that have phases of near steady high-activity before oscillations start, as sometimes seen experimentally (Fig. 1, bottom panel).

Figure 7B shows the variations of oscillation period with theta . For the lowest possible value of theta  allowing oscillations, the period is finite. The period lengthens as theta  increases, first slowly and then abruptly just before the oscillations disappear (it becomes infinite at the exact point where the low oscillatory level meets the intermediate, unstable state). Except for the very end of the range, the cycle period ranges between 5 and 10 units of time (approximately 7 for the activity shown on Fig. 6C). Therefore, because the typical cycle period measured experimentally is ~1 sec for embryonic day (E) 10 embryos, tau d = 2 predicts that the time constant of the depression variable d is of the order of 200-400 msec. This is comparable to time scales reported in other systems (Chance et al., 1998).

Slow activity-dependent depression leads to episodic behavior

In the chick embryo, after an episode of activity, evoked monosynaptic and polysynaptic potentials are depressed and recover with a time constant of ~1 min (Fedirchuk et al., 1999). In addition, after an episode, the membrane potential of ventrally located spinal neurons hyperpolarizes and subsequently recovers as a depolarizing ramp until the next episode. Collectively, these observations indicate the presence of a slow (time scale minutes) activity-dependent form of depression that regulates network excitability and has been proposed to control the occurrence of episodes (O'Donovan and Chub, 1997). Although the mechanisms for these changes are unknown, they could occur through a modulation of cellular properties or a prolonged "slow" form of synaptic depression. In the following, we formulate and analyze models for each of these two possibilities.

The two-variable system presented in the previous section models a bistable, recurrent excitatory network that can exhibit persistent states of low and high activity, the latter being oscillatory. In contrast, the chick's spinal cord shows spontaneous alternation between two such states. Here, we add a slow, third variable that enables the model to generate sustained, repetitive episodes of cycling at high activity separated by low-activity phases, thereby mimicking the spontaneous, episodic rhythm of the cord. We show that the mathematical structure that underlies the episodic behavior is the same for both types of slow activity-dependent depression. Therefore, the two models behave in similar ways, say, when perturbed by a brief stimulation or noisy input, and they are affected similarly by some parameter variations. In a later section, however, we show that the models lead to different predictions for variations in connectivity and therefore the models can be distinguished experimentally.

Slow modulation of cellular excitability: theta -model

Intracellular recordings have shown that there is a slow ramp depolarization in motoneurons of chick embryonic spinal cord (N. Chub and M. J. O'Donovan, unpublished observation), suggesting that after an episode the functional threshold (theta ) of neurons is raised and slowly decreases until the next episode. We can therefore imagine that theta  increases during an episode and decreases during the silent phase. This means that the growth and decay of theta  is activity-dependent, analogous to a slow adaptation mechanism. During high activity the threshold rises, eventually terminating the activity, and then it recovers during the interepisode quiet phase. When the threshold decreases adequately, the network's small amount of spontaneous activity can trigger another episode. We formulate below a kinetic model for theta  that yields this behavior.

From our dynamical systems viewpoint and the previous sections, we can predict the model system's behavior as follows (we will confirm this below in Figure 8 with a simulation). First, because theta  will move very slowly, we expect that the a - d subsystem will evolve relatively rapidly to one of its two attractors (Fig. 7A) and then track it as theta  moves. Suppose the network is in the low-activity state, as after an episode. Under this condition theta  will slowly decrease, and the a - d subsystem will slowly track the Z-curve's lower branch until the left knee is reached. The system then rapidly moves to the oscillatory high-activity state, starting an episode. During the episode theta  now increases, and as indicated by the thick gray curves on Figure 7A, the cycle amplitude also gradually increases. Moreover, the cycle period increases too (Fig. 7B), but only slightly at the beginning of the episode, then significantly near the end of the episode. Note, this significant drop in period just before episode termination is a generic feature of this type of model and is typical of many experiments in chick cord (Fig. 1, top panel). Eventually, theta  reaches the critical level theta c where activity can only jump back to the lower steady state, terminating the episode. As theta  then decreases again (recovery), the process will be repeated. Thus, the presence of bistability of the a - d subsystem for a range of theta  values suggests a mechanism for the slow occurrence of episodes, in the same way that the bistability of the simple recurrent network over a range of connectivity values suggested a mechanism for the generation of cycles during an episode. Our prediction of the slow rhythm's a - theta trajectory rests on the assumptions (i.e., on our construction) that (1) the fast a - d subsystem is bistable, (2) theta  evolves very slowly, and (3) theta  evolves in the proper directions during the active and quiescent phases. These features form the essence of the fast/slow dissection technique for analyzing multiple time scale oscillations and excitability of this sort (Rinzel and Ermentrout, 1998).



View larger version (23K):
[in this window]
[in a new window]
 
Figure 8.   A, Time courses of a and theta  for the theta -model. During an episode theta  increases, until the episode stops. B, Corresponding trajectory in the (atheta ) plane, superimposed on the bifurcation diagram. Also plotted is the theta -nullcline (dot-dashed curve). Below this nullcline, the trajectory flows leftward (<A><AC>&thgr;</AC><AC>˙</AC></A> < 0, recovery during the silent phase); above it the trajectory flows rightward (<A><AC>&thgr;</AC><AC>˙</AC></A> > 0, depression during an episode). The points numbered 1, 2, 3, and 4 correspond to the same states of the network on both A and B. During an episode, the trajectory passes sequentially 1-2-3-4; during a silent phase it goes from 4 to 1.

As for a and d, we model the variations of theta  with a simple, first-order kinetics:
&tgr;<SUB>&thgr;</SUB><A><AC>&thgr;</AC><AC>˙</AC></A>=&thgr;<SUB>∞</SUB>(a)−&thgr;, ()
with theta infinity an increasing (sigmoidal) function of a and tau theta being two to three orders of magnitude greater than tau a. Thus theta  will increase from low values if activity is high and decrease from high values if activity is low.

The time courses of a and theta  for the full three-variable model resulting from the interaction of the a - d fast subsystem and the slow theta  subsystem are shown in Figure 8A. As expected, the behavior is episodic, with rhythmic cycling during episodes. The system's trajectory actually lives in the three-dimensional phase space, a - d - theta . In Figure 8B we project the trajectory onto the (a - theta ) plane, enabling us to see the dynamic behavior along with the static bifurcation diagram from Figure 7A. We see qualitatively what we anticipated in the preceding paragraphs, with theta  covering the bistable range back and forth. Although the trajectory does not exactly follow the bifurcation diagram, it would if we had made tau theta larger.

It is important to emphasize that both onset and termination of the episodes are controlled by the value of theta , as shown from the trajectory in the (a - theta ) plane. Also shown in Figure 8B is the theta -nullcline, the set of points for which <A><AC>&thgr;</AC><AC>˙</AC></A> = 0, that is, theta  = theta infinity (a). Below this nullcline, <A><AC>&thgr;</AC><AC>˙</AC></A> < 0, therefore theta  is decreasing, and conversely, for any point above the nullcline theta  is increasing. Thus, theta  decreases during the silent phase and increases during an episode. Obviously, the position of the theta -nullcline relative to the a - d bifurcation diagram is critical for the existence and timing of episodic rhythmicity (see below).

This graphical representation allows one to easily predict/understand the effect of parameter variations on the episodic rhythm. Suppose the cellular adaptation is recruited more easily (say, theta theta is decreased). This lowers the theta -nullcline. If it drops too much it will intersect the Z-curve's lower branch. At this intersection, all variables, including theta , are at a stable steady state of low activity; the network will remain silent. On the other hand, if the theta -nullcline is too high, theta  can stop increasing before the end of the episode, and the network might enter continuous oscillation. Therefore, the position of the nullcline, or in other words, the value of a for which the average cellular excitability in the network starts to get depressed, has to be within a certain range for episodic oscillations to occur.

The episode durations and silent phase lengths depend on two factors. First, they depend on the range of values that theta  covers between the beginning and end of the episodes, that is, the horizontal distance between points 4 and 1 on Figure 8B. Increasing this range dilates both the episode and silent phase in the same way, because the paths 1-2-3-4 and 4-1 are increased accordingly. As illustrated below, this range can be changed by varying any single parameter of the fast subsystem. Second, the episode durations and silent phase lengths depend on the velocity with which theta  varies. Because this velocity is inversely proportional to tau theta , an increase in this time constant increases similarly the time it takes to depress the network and the recovery period, so both duration and silent phase change by the same factor. However, this velocity is also proportional to the horizontal distance between the current point in the a - theta plane and the theta -nullcline, because tau theta <A><AC>&thgr;</AC><AC>˙</AC></A> = theta infinity (a- theta . Therefore, if the nullcline is changed [that is, the function theta infinity (a) is changed], velocities are modified. However, this change in the velocity depends on the activity state (high or low) of the system. For example, if the nullcline is moved down, the velocity decreases during the silent phase because theta infinity (a- theta is reduced, but increases during the episode because theta infinity (a- theta is increased. Therefore, episode and silent phase lengths may be affected in opposite ways when a parameter of the slow subsystem is changed. The distinction between the effects of parameters of the fast versus slow subsystems will be important subsequently to distinguish the two models.

Long-term synaptic depression: s-model

As was considered previously, it is also possible to interpret the slow form of activity-dependent depression as arising synaptically rather than from a slow modification of cellular thresholds. Indeed, this mechanism has been proposed to account for the slow recovery of synaptic potentials after an episode (Fedirchuk et al., 1999). Calling s the new, slow variable, the system becomes:
&tgr;<SUB><UP>a</UP></SUB><A><AC>a</AC><AC>˙</AC></A>=a<SUB>∞</SUB>(n·s·d·a)−a

&tgr;<SUB><UP>d</UP></SUB><A><AC>d</AC><AC>˙</AC></A>=d<SUB>∞</SUB>(a)−d

&tgr;<SUB><UP>s</UP></SUB><A><AC>s</AC><AC>˙</AC></A>=s<SUB>∞</SUB>(a)−s.
As for d, a high value of s (near 1) means "not depressed," so sinfinity is a decreasing function of a. Note that the fraction of nondepressed synapses is now s·d, and the effective connectivity is n·s·d. The activity obtained with this model is qualitatively similar to the one obtained with the theta -model, as shown in Figure 9. Here again, the onset and termination of episodes are controlled by the slow depression variable (s), therefore silent phase and episode durations are determined in the same way. The reason the theta -model and the s-model are so similar is that variations of s induce the same qualitative changes in the dynamics of the fast subsystem as do variations of theta . Note that because depression is associated with smaller values of s, in contrast to the theta -model where large theta  means larger depression, the bifurcation diagram of Figure 9B has opposite (left/right) orientation to that of Figure 8B.



View larger version (23K):
[in this window]
[in a new window]
 
Figure 9.   A, Time courses of a and s for the s-model. During an episode s decreases, until the episode stops. As for the theta -model (Fig. 7B), the cycle period increases near the end of an episode. A little after t = 400, a stimulation was applied, which triggered an episode. This episode was shorter because recovery of s was not complete. B, Corresponding trajectory in the (as) plane, superimposed with the bifurcation diagram and the s-nullcline (dot-dashed curve). C, Recordings from an E10 cord. Stimulating quickly after a spontaneous episode triggers a short episode. Black triangles indicate time of stimulation.

Simulation experiments with the models

To further characterize the models' behaviors and compare them with experimental data, a few manipulations were conducted. The model network can be perturbed by a brief stimulation. It is also important to understand how the model is affected by varying its parameters, because the parameter variations may be identified with experimental manipulations or with evolving system conditions that occur during development. In addition to the comparison with experimental data, the following simulations offer some predictions that are testable experimentally.

Activity resetting. In the chick spinal cord, it is easy to trigger an episode during a silent phase with a brief electrical stimulus. In addition, if the stimulation is applied shortly after a spontaneous episode, the evoked episode is much shorter than the average spontaneous episode, as illustrated on Figure 9C. This effect was observed on all preparations tested (>20). It is to be expected if the silent phase is a period of recovery from depression, and it is indeed a characteristic feature of the models presented here. We have illustrated this point for the s-model on Figure 9A,B. Again exploiting the graphical representation of Figure 9B, we see that to evoke an episode the stimulus must be substantial enough to bring the activity above the S-curve's middle branch. Otherwise, activity returns immediately to the low steady state. Because the vertical distance between the unstable steady state and low-activity state is highest just after an episode and progressively declines, the smaller the amount of time after an episode, the harder it is to successfully evoke a new episode, in agreement with experiments. These results are also observed with the theta -model because it relies on the same mechanism of episode generation.

Changing a parameter of the fast subsystem. The main effect of changing one parameter of the fast a - d subsystem is to change the values of the slow variable at the beginning and/or end of episodes (or, equivalently, at the end and beginning of silent phases) by deforming the shape of the bifurcation diagram and, in many cases, without moving it much relative to the slow variable's nullcline. Therefore, as mentioned above (section describing the theta -model), silent phase and episode duration are affected in the same way, because the slow variable runs through the same range during the episode and the silent phase.

Figure 10 illustrates the effects of changing parameters in the fast subsystem: the steepness of a<