help button home button Biophys. J.
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
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 Google Scholar
Google Scholar
Right arrow Articles by Oprisan, S. A.
Right arrow Articles by Canavier, C. C.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Oprisan, S. A.
Right arrow Articles by Canavier, C. C.
Biophysical Journal 84:2919-2928 (2003)
© 2003 The Biophysical Society

Dynamics from a Time Series: Can We Extract the Phase Resetting Curve from a Time Series?

S. A. Oprisan *, V. Thirumalai {dagger} and C. C. Canavier *

* Department of Psychology, University of New Orleans, New Orleans, Louisiana; and {dagger} Volen Center for Complex Systems, Brandeis University, Waltham, Massachusetts

Correspondence: Address reprint requests to Carmen C. Canavier, University of New Orleans, Dept. of Psychology, 2000 Lakeshore Dr., GP 2001, New Orleans, LA 70148. Tel.: 504-280-6775; Fax: 504-280-6049; E-mail: ccanavie{at}uno.edu.


   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION AND CONCLUSIONS
 APPENDIX
 ACKNOWLEDGEMENTS
 REFERENCES
 
Recordings of the membrane potential from a bursting neuron were used to reconstruct the phase curve for that neuron for a limited set of perturbations. These perturbations were inhibitory synaptic conductance pulses able to shift the membrane potential below the most hyperpolarized level attained in the free running mode. The extraction of the phase resetting curve from such a one-dimensional time series requires reconstruction of the periodic activity in the form of a limit cycle attractor. Resetting was found to have two components. In the first component, if the pulse was applied during a burst, the burst was truncated, and the time until the next burst was shortened in a manner predicted by movement normal to the limit cycle. By movement normal to the limit cycle, we mean a switch between two well-defined solution branches of a relaxation-like oscillator in a hysteretic manner enabled by the existence of a singular dominant slow process (variable). In the second component, the onset of the burst was delayed until the end of the hyperpolarizing pulse. Thus, for the pulse amplitudes we studied, resetting was independent of amplitude but increased linearly with pulse duration. The predicted and the experimental phase resetting curves for a pyloric dilator neuron show satisfactory agreement. The method was applied to only one pulse per cycle, but our results suggest it could easily be generalized to accommodate multiple inputs.


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
METHODS
 RESULTS
 DISCUSSION AND CONCLUSIONS
 APPENDIX
 ACKNOWLEDGEMENTS
 REFERENCES
 
Central pattern generators (CPGs) are networks of neurons that are capable of producing rhythmic output (Beer et al., 1999Go; Chiel et al., 1999Go; Golubitsky et al., 1998Go; Kopell and Ermentrout, 1988Go; Marder and Calabrese, 1996Go) in the absence of both input from higher centers and sensory feedback. CPGs are believed to mediate certain rhythmic, stereotyped behaviors, such as walking, flying (Robertson and Pearson, 1985Go), swimming (Arshavsky et al., 1985Go; Satterlie, 1989Go), chewing, feeding, and scratching (Pearson, 1993Go). The pyloric network of the stomatogastric ganglion is an example of a well-studied CPG that produces a three-phase rhythm driving muscles in the stomach of lobsters and crabs (Selverston and Moulins, 1987Go). This CPG develops a stable rhythmic pattern over a frequency range from 0.3 to 3 Hz. The pacemaker unit is formed by the anterior burster (AB) neuron electrically coupled with two pyloric dilator (PD) neurons. The isolated AB neuron oscillates at a frequency in the normal range of AB/PD network (Bal et al., 1988Go; Hooper and Marder, 1987Go; Marder and Meyrand, 1989Go; Miller and Selverston, 1982Go) whereas the isolated PD neuron has irregular oscillations with a much longer period than that of the normal network (Abbott et al., 1991Go).

Bursting activity is characterized by slow oscillations in membrane potential that alternate between a silent hyperpolarized interburst and a burst characterized by the firing of a succession of action potentials (Baxter et al., 2000Go; Butera et al., 1995Go; Canavier et al., 1991Go; Chay and Keizer, 1983Go; Ermentrout, 1986Go, 1996Go; Liu et al., 1998Go). The dynamics are comprised of at least two distinct time scales, often including a slow one associated with the underlying subthreshold oscillations in the membrane potential, as well as a fast one associated with the generation of action potentials. Analytical techniques take advantage of the separation in the time scales, although the separation is not complete because an action potential may perturb the slow state variables that in turn control the firing of action potentials. Bursting neurons can be modeled as nonlinear oscillators, and some of the behavior of coupled nonlinear oscillators is independent of at least some details of the properties of the individual neurons and of coupling, but rather is characteristic of the general architecture of the network (Canavier et al., 1997Go, 1999Go; Collins and Stewart, 1993Go; Collins and Richmond, 1994Go; Kopell and Ermentrout, 1988Go). As a component of a network, nonlinear oscillators can be characterized by their phase resetting curve (PRC; see Abramovich-Sivan and Akselrod, 1998aGo,bGo; Dror at al., 1999Go; Ermentrout, 1996Go; Perkel et al., 1964Go; Pinsker, 1977Go). A PRC tabulates the effects on the period of the current cycle as a result of brief perturbations applied to an uncoupled oscillator at various phases in the burst (Canavier et al., 1997Go, 1999Go; Ermentrout, 1985Go; Ermentrout and Kopell, 1991Go; Kopell and Ermentrout, 1988Go; Murray, 1993Go; Pavlides, 1973Go). To predict the behavior of the circuit from the response of the uncoupled oscillator, the perturbation used to generate the PRC must resemble as closely as possible the input that the oscillator would receive in the circuit. A powerful computational tool based on the PRC method was developed by Canavier and co-workers (Canavier et al., 1997Go, 1999Go) to investigate stable pattern emergence in ring circuit dynamics. According to Winfree (1980)Go, there are two types of PRCs—type 1 and type 0. Ermentrout (1996)Go subdivided Winfree's type 1 PRCs into Type I and Type II, which he then correlated with different excitability characteristics.

The characteristics of a type 1 PRC are relatively small phase shifts and a continuous transition between delay and advance. Type 0 PRC displays larger phase shifts and discontinuities between delays and advances, and there exist phases that will never result from a perturbation. Type 0 PRC results from a hard perturbation, which requires a stronger resetting stimulus (intensity and/or duration). Our previous work (Oprisan and Canavier, 2001Go, 2002Go) and that of others focused on type 1 resetting, whereas in the experimental data presented in this article we encountered type 0 resetting.

As a first step toward understanding the oscillatory patterns produced by such circuits, we study how inputs from other neurons in the circuit affect the frequency of bursting neurons. As an example, we have chosen to examine the effects of simulated synaptic perturbations on the PD neuron, which is electrically coupled to the AB neuron. Although the AB neuron sets the rhythm for the pyloric circuit in the stomatogastric ganglion of the lobster Homarus americanus, the PD neuron is experimentally more accessible. We have assumed that the electrically coupled pair can be considered as a functional unit and a single nonlinear oscillator for the purpose of our analysis. The theoretically predicted phase resetting for a perturbation that takes the form of a pulse in synaptic conductance was compared with the experimental results obtained using conductance pulses of various amplitudes and durations. The present study addresses the PRC prediction in the case of hard inhibitory perturbations. The preliminary results regarding the applicability of the method to excitatory perturbations are encouraging but we need further experimental data for a quantitative comparison with the theoretical predictions. We considered only inhibitory perturbations in this study because they predominate in the pyloric circuit. In practice, only perturbations strong enough to hyperpolarize the membrane below its most hyperpolarized level in a free running experiment were considered. Thus we were guaranteed hard perturbations during the depolarized phase, and the hard nature of the perturbation was confirmed because the old phase mapped only to a subset of the possible new phases. Such hard perturbations are physiologically realistic in the context of this circuit. The resulting theoretical predictions were within the limits of the variability in the experimental data.

Previously, PRCs have been utilized to analyze certain simple circuits, specifically unidirectional rings comprised entirely of bursting neurons (Canavier et al., 1997Go, 1999Go). The unidirectional ring architecture arises from the assumption that each neuron receives input from exactly one other neuron during each cycle. The methods developed by Canavier et al. (1997Go, 1999)Go enable the design of circuits capable of generating multiple patterns of periodic behavior and offer efficient mechanisms for the control of multistability. Pattern prediction based on the PRC assumes that the complete structure of the solution space of the unidirectional ring can be predicted from the PRC. The two basic hypotheses on which this method is built are: 1), the input received in the circuit by each component oscillator from its presynaptic oscillator has the same effect in the closed loop circuit condition as an identical input delivered in an open loop condition such as that used for generating PRC curves, and 2), the effect of such an input effectively dies out by the time the next input is received. Although the ring geometry has been utilized extensively in biophysical modeling studies (Collins and Stewart, 1993Go; Collins and Richmond, 1994Go), we would like to extend these methods to include bidirectional coupling and arbitrary geometries that allow for more than one input per neuron per cycle, and for the inclusion of neurons that do not burst endogenously when isolated from the circuit. A first step in generalizing circuit analysis is to understand the phenomenology of phase resetting so that we can predict the effect of the multiple and/or overlapping inputs within a single cycle.

We model the PD/AB complex as a limit cycle oscillator, a simple periodic oscillator in which all variables (such as membrane potentials) repeat themselves exactly during each periodic cycle. A convenient way to study periodic behavior is by using a phase representation. The term phase has been utilized in the literature in several different ways. Here we use the definition of phase given by Winfree (1980)Go: the elapsed time measured from an arbitrary reference divided by the intrinsic period. A limit cycle oscillation requires a minimum of two variables (Guckenheimer and Holmes, 1983Go; Hoppensteadt and Izhikevich, 1997Go; Murray, 1993Go; Winfree, 1980Go). The limit cycle oscillator can exhibit a range of dynamic activity that depends upon the characteristic time scales of the two variables (Rinzel and Lee, 1986Go; Rinzel and Ermentrout, 1998Go). If they are comparable in magnitude, the qualitative dynamics of a phase oscillator arise, and the two variables tend to vary smoothly in unison. If, however, one variable varies much more slowly than the other, a relaxation oscillator can be produced, in which alternating rapid increases and decreases in the fast variable (jumps) are followed by slow relaxation in the slow variable. We found that the PD/AB complex is best modeled as a relaxation oscillator with two solution branches, one depolarized (the burst) and one hyperpolarized, with hysteresis with respect to the slow variable. Previous studies (Oprisan and Canavier, 2001Go, 2002Go) conjectured that the PRC could be retrieved from a time series record of membrane potential. We have developed a mapping method that allows us to determine the PRC from a membrane potential time series, based on the assumed phase model described above.

In the present study, we use a consistent mathematical definition of phase resetting (Canavier et al., 1997Go, 1999Go). We choose the instant at which the neuron reaches the spiking threshold as the reference point. A synaptic conductance perturbation applied at time t changes the intrinsic period Pi of the current cycle to P1. In other words, a stimulus applied at (temporal) phase changes the duration of the current cycle to P1({varphi}). The normalized first-order PRC is so that a positive value represents a phase delay and a negative value a phase advance (Fig. 1).



View larger version (21K):
[in this window]
[in a new window]
  		FIGURE 1  The membrane potential record from a free running PD neuron with an intrinsic period of oscillation, P0 (continuous line), is perturbed at time t by applying an inhibitory synaptic conductance pulse. As a result, the next burst is advanced and occurs after P1 < P0 (dashed line, only first spike shown). The perturbation lasts for 250 ms (horizontal bar) and its strength is 100 nS.

 
The terms phase advance and delay are based on a conceptual framework developed for soft perturbations that reset the phase tangentially along the limit cycle: advances move the trajectory in the direction of its natural motion, speeding up the limit cycle, whereas delays move the trajectory in the opposite direction, slowing down the limit cycle. Previous theoretical work by ourselves (Oprisan and Canavier, 2002Go) and others (Ermentrout, 1996Go; Izhikevich, 2000Go) focused on such tangential perturbations. The hard perturbations in this study caused the trajectory to move in a direction that was normal to the limit cycle. Thus the concept of advance or delay does not apply in the original sense of tangential motion, but the definition that we use preserves the idea that an advance causes the next burst to occur sooner, whereas a delay causes the next burst to occur later.


   METHODS
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
RESULTS
 DISCUSSION AND CONCLUSIONS
 APPENDIX
 ACKNOWLEDGEMENTS
 REFERENCES
 
Experimental methods
Adult Homarus americanus were obtained from Commercial Lobsters, Boston, MA and maintained in artificial sea water tanks at 11°C. Lobsters were anesthetized on ice for 15 min and then the complete stomatogastric nervous system consisting of the paired commissural ganglia, esophageal ganglion, and the stomatogastric ganglion and their connecting and motor nerves was dissected and pinned out in a transparent Sylgard coated dish (Dow Corning, Midland, MI) containing chilled (9–13°C) saline. Saline composition (in mM): 479.12 NaCl; 12.74 KCl; 13.67 CaCl2; 10 MgSO4; 3.91 Na2SO4; 5 HEPES; pH 7.45 (Richards et al., 1999Go). The stomatogastric ganglion was desheathed and Vaseline wells were made on the motor nerves for extracellular nerve recordings. The preparation was continuously superfused with chilled saline. The lateral pyloric and the pyloric dilator (PD) neurons were identified using standard procedures (Selverston and Moulins, 1987Go).

The lateral pyloric neuron was filled with Lucifer Yellow (Molecular Probes, Eugene, OR) with -20 nA for 45 min and then illuminated with blue light from a 100 watt mercury lamp through the MZFLIII Leica binocular microscope to photoinactivate the neuron (Miller and Selverston, 1979Go). A PD neuron was impaled with two microelectrodes filled with 0.6 M K2SO4 plus 20 mM KCl (20–30 M{Omega} resistance). The dynamic clamp (Pinto et al., 2001Go; Sharp et al., 1993Go, 1996Go) was used to create an artificial inhibitory synaptic conductance in the PD neuron during different phases of its burst using software written in LabWindows (National Instruments, Austin, TX). Square pulses in gsyn of various amplitudes and durations were applied to generate phase response curves. The maximal synaptic conductance gsyn used varied from 100 nS to 1000 nS. The current injected was calculated from the conductance using the relationship

where Erev = -90 mV. Data were acquired using a Digidata 1200 data acquisition board (Axon Instruments, Foster City, CA) and analyzed using programs written by Dr. Bill Miller (Vegasci, http://www.vegasci.com) in MATLAB (The MathWorks, Natick, MA). Phase resets were calculated offline and plotted using Microsoft Excel. A phase of zero was assigned to the spiking threshold (~-47 mV).

Theoretical methods
Previous theoretical work on phase resetting (Oprisan and Canavier, 2001Go, 2002Go) focused on mathematical models of oscillators, for which the system equations are known. This type of approach is not particularly useful in the analysis of neural oscillators, for which the system equations are not known precisely. On the other hand, it is well known that the geometry of the solution structure of physical oscillators can be recovered from a time series via delay embedding reconstruction of the attractor (Casdagli et al., 1991Go; Eckmann and Ruelle, 1985Go; Fraser and Swinney, 1986Go; Guckenheimer and Holmes, 1983Go; Hegger at al., 1999Go; Packard et al., 1980Go; Schouten et al., 1994Go; Takens, 1981Go; Wolf et al., 1985Go). We found that such a reconstruction provides insights regarding the phenomenology of phase resetting in neural oscillators.

We plan to characterize the AB/PD complex as a limit cycle oscillator, and a limit cycle oscillator requires at least two variables. Because we have access to only a single variable, the time series recording of the membrane potential, we need a second variable, preferably a slow one to capture the relaxation oscillator dynamics that characterize the AB/PD complex. We used both a low-pass filtered version of the membrane potential, and a spline fit to the minima during the burst to produce an envelope. Both methods gave similar results, so we used the envelope of the membrane potential to give us a second, slower, variable-to-plot against the membrane potential (Fig. 2). We justify the use of only two variables in the Appendix. Briefly, we used the package TISEAN (http://www.mpipks-dresden.mpg.de) that utilized time-delayed versions of the voltage system to reconstruct a limit cycle attractor corresponding to bursting. The embedding dimension was revealed to be three, and the spikes caused a displacement in each one of the time-shifted versions of voltage (see Appendix). The underlying envelope was two-dimensional, however, so we reduced the system to two variables, enabling us to assume that the perturbations only affected one variable (membrane voltage), which greatly simplifies the analysis. To better capture the relaxation oscillator characteristics of bursting in the AB/PD complex, we used a version of the membrane potential envelope shifted by the time lag suggested by the TISEAN software as the slow variable, and the membrane voltage as the fast variable (Fig. 2).



View larger version (18K):
[in this window]
[in a new window]
  		FIGURE 2  The time delayed envelope of the membrane potential Ve(t) provides the slow (envelope) variable (dashed line) for the reduced two-dimensional attractor reconstruction (A). The two-dimensional attractor retains the spikes of the fast variable (membrane potential, V(t)) and is smooth in respect with the slow variable (B). The optimum time lag for the slow variable was found using the package TISEAN and is set, in this case, to {Delta}t = 340 ms.

 

   RESULTS
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
DISCUSSION AND CONCLUSIONS
 APPENDIX
 ACKNOWLEDGEMENTS
 REFERENCES
 
Predictions for short pulses
We were able to predict the PRC for short pulses using the reconstructed limit cycle attractor by making the following assumption, that a strong inhibitory perturbation instantaneously switches the trajectory from the depolarized (spiking) branch on the right to the hyperpolarized (silent) branch on the left. Thus in Fig. 3 A we have plotted the reconstructed limit cycle, and the dotted line indicates the effect of applying a sufficient perturbation at point 1 (phase {approx} 0.21), which is to shift the phase to 0.67 (point 2) on the hyperpolarized branch of the limit cycle. The neuron will be able to spike again after 0.33 instead of 0.79 of its intrinsic period. This gives the -0.46 phase advance plotted in Fig. 3 C for a phase of 0.21. Fig. 3 B shows a mapping between the phase on the depolarized ({varphi}D) and hyperpolarized ({varphi}H) branches of the limit cycle at a constant value of the slow envelope Ve(t). Thus only V(t) and not Ve(t), the assumed slow variable, is assumed to be perturbed by a hyperpolarization. The normal component of phase resetting ({Delta}{varphi}N) is given by {Delta}{varphi}N = {varphi}H - {varphi}D on the depolarized branch, and by {Delta}{varphi}N = 0 on the hyperpolarized branch. This assumes the trajectory quickly relaxes to the nearest branch (in this case the hyperpolarized branch) after a perturbation.



View larger version (15K):
[in this window]
[in a new window]
  		FIGURE 3  The instantaneous PRC. The superthreshold perturbation of the limit cycle dynamics (the horizontal heavy dotted line that connects the points marked 1 and 2 on A) shifts the figurative point 1 from the depolarized branch to the hyperpolarized branch of the unperturbed limit cycle, point 2 on A. The vertical dotted line on A separates the depolarized and hyperpolarized branches of the limit cycle. The mapping between the phases on the depolarized ({varphi}D) and hyperpolarized ({varphi}H) branches of the limit cycle (B) contains two distinct regions: a), between the onset of the burst and the end of maximum value of slow variable, which corresponds to a phase advance; and b), between the minimum of the slow variable and the onset of the burst, which corresponds to a phase delay. The dotted line on B is plotted along zero phase resetting region. The phase resetting due to normal displacement ({varphi}N) is a large phase advance for points (phases) on the depolarized branch and nearly zero phase resetting for most of the hyperpolarized branch (C).

 
There is good agreement between the predicted and experimental phase resetting curves for short pulses (Fig. 4). Fig. 4, A and B are for 100 and 250 ms pulses in the same neuron. Fig. 4 C is for a 200 ms pulse in a second neuron. Note that the PRC is flat (zero phase resetting) in the region that corresponds to the hyperpolarized branch {varphi} (0.6, 0.9). Presumably, additional hyperpolarization that does not delay the jump to the depolarized branch (see next section) does not appreciably affect the time course of the assumed slow variable. One candidate for a slow process during hyperpolarization is the calcium dynamics. We can infer that if calcium is indeed the slow variable, then the additional hyperpolarization does not significantly change the rate of the removal of calcium accumulated during a burst. The limit cycles shown in Figs. 2 and 4 are for the first neuron, and analogous figures apply to the second neuron. The PRCs resulting from the application of a synaptic pulse of 100 ms duration with a pulse strength in the range from 100 nS to 1000 nS displays a good agreement with our predicted PRC (Fig. 4 A). It is striking that the amplitude of the pulse does not seem to affect the amount of phase resetting produced. This threshold effect, in that a pulse is either sufficient to cause a jump between branches or it is not, suggests a rule for predicting the effect of simultaneous pulses. If they are both subthreshold, then neither will have an appreciable effect, but if their sum is above threshold, a large phase resetting will be produced. If they are both above threshold then the combined pulse will have the same effect as a single pulse. The stimulus strength must be sufficient to switch branches, but otherwise does not affect the PRC. If the membrane potential perturbation is strong enough to escape from the depolarized branch to the hyperpolarized branch, then the PRC can be predicted using a mapping between points on the limit cycle that have the same value of the slow variable (Fig. 3 B).



View larger version (19K):
[in this window]
[in a new window]
  		FIGURE 4  Experimental and predicted PRCs for the instantaneous normal component only. For short synaptic conductance pulses and a wide range of synaptic strengths the agreement between experimental and predicted phase resetting is satisfactory. A and B refer to the data from the first neuron. C is the common plot of the experimental and predicted PRCs for the second neuron. The pulse duration {tau} ranged from 100 ms to 250 ms with the synaptic strength gsyn between 100 and 1000 nS. As the pulse duration increases, the instantaneous PRC fails to give an accurate prediction of the phase resetting. Data from five experiments were analyzed but here only the results from two of them are presented. The other experiments gave similar results.

 
Predictions for long pulses
The physiologically realistic synaptic perturbations are far from infinitesimal. For pulses longer than 250 ms the fit was poor using only {Delta}{varphi}N (not shown). As a result, we corrected the above relationship of phase resetting {Delta}{varphi}T to take into account the effect of the duration of the perturbation {tau}. The finite duration of the perturbation induces additional delay to the total phase resetting. Since the hyperpolarizing pulse is strong, it will not allow the neuron to jump from the hyperpolarized to depolarized branch until the perturbation terminates. Thus if the recovery time (time after the normal resetting until the next burst is expected 1 - ({varphi}s + {Delta}{varphi}N), where {varphi}s is the phase at which the stimulus was applied and Pi is the intrinsic period) is shorter than the normalized pulse duration {tau}/Pi, then an additional component of phase resetting {Delta}{varphi}w must be added for the time spent waiting for the perturbation to end. Therefore, {Delta}{varphi}w = {tau}/Pi - (1 - ({varphi}s + {Delta}{varphi}N)) if the recovery time is shorter than the pulse duration, and {Delta}{varphi}w = 0 if it does not. The total phase resetting is then given by {Delta}{varphi}T = {Delta}{varphi}N + {Delta}{varphi}w. The above equations produced a very good agreement of the predicted PRCs with the experimental results even for very long pulses (Fig. 5). Fig. 5, A–D are from the first neuron, Fig. 5 E is from the second, and a similar procedure was applied with similar results to data from three other neurons (not shown) for a total of five neurons (n = 5). Note that Fig. 5 B uses the same data as Fig. 4 B but produces a better fit due to the {Delta}{varphi}w correction, and the same applies to Fig. 5 C and Fig. 4 C. The main effect of an inhibitory pulse longer than the recovery time consists in shrinking the region of zero phase resetting on the hyperpolarized branch. This happens because such pulses will always induce a burst as the neuron is released from inhibition upon pulse termination. For an inhibitory perturbation with a pulse duration equal to the silent regime (interburst) duration, no matter what the phase when the perturbation is applied, the phase after perturbation is zero (burst initiation), resulting in a constant cophase (new versus old phase) plot and a linear PRC (Demir at al., 1997Go). For pulse durations equal or longer than the interburst interval the PRC is a continuous, linear curve (the horizontal region of zero phase resetting progressively shrinks until it completely disappears). Moreover, if the pulse duration exceeds the duration of the interburst interval, the PRC is linear and vertically shifted with a phase resetting equal to the difference between the normalized pulse duration and the normalized interburst interval. We can now formulate a rule for adding pulses whose durations overlap or are sequential, regardless of amplitude, providing that the combined amplitude is always suprathreshold. The total duration of the summed pulses is the effective duration, and the total resetting can be calculated by using {Delta}{varphi}T = {Delta}{varphi}N + {Delta}{varphi}w, with {Delta}{varphi}N due only to the first pulse but {Delta}{varphi}w due to the combined duration.



View larger version (30K):
[in this window]
[in a new window]
  		FIGURE 5  Experimental and predicted PRCs using both components of resetting. The total phase resetting {varphi}T is given by the sum of the normal resetting {varphi}N (as for instantaneous PRCs) and the additional delay induced by finite duration of the perturbations {varphi}w. The agreement between experimental and predicted phase resetting is good. Resetting is insensitive to the value of gsyn but quite sensitive to duration of the perturbation {tau}. A–D are the PRCs for the first neuron. E is the PRC for the second neuron. The pulse duration {tau} ranged from 100 ms to 750 ms with the synaptic strength gsyn between 100 and 1000 nS. Data from five neurons were analyzed but here only the results from two neurons are presented. Analysis of the data from the other neurons gave similar results.

 

   DISCUSSION AND CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION AND CONCLUSIONS
APPENDIX
 ACKNOWLEDGEMENTS
 REFERENCES
 
The present study is a computational approach to PRC reconstruction based entirely on the membrane potential record. We showed that, for the case of hard inhibitory perturbations, we can successfully reconstruct the PRC using the experimentally recorded membrane potential from a PD neuron. These PRCs suggest that the effective control parameter that determines the shape of the PRC is not pulse amplitude but pulse duration. Demir and co-workers (Demir et al., 1997Go), based on a model study, also observed that "stimulus amplitude has a minimal effect on the shape of the PRC, whereas stimulus duration evokes significant changes in the PRC." Demir and co-workers (Demir et al., 1997Go) also used a phase plane representation of the limit cycle, but their analysis was more complicated because their model had not just one slow variable, as we assumed, but two. If the membrane potential perturbation {Delta}V is weak (subthreshold) then the PRC might be obtained using tangential displacement along the limit cycle (Ermentrout, 1996Go; Kopell and Ermentrout, 1988Go; Oprisan and Canavier, 2001Go, 2002Go), but here we have not examined this possibility. There is an additional situation to consider. A subthreshold pulse may move the trajectory near threshold for a switch to the hyperpolarized branch without actually crossing the threshold. In this case, the trajectory may travel slowly back to the depolarized branch causing a significant delay. This type of slow relaxation is difficult to deal with analytically (Oprisan and Canavier, 2001Go) and may limit the applicability of the method presented in this article.

Consistent with the description in our results of how the flat region in the PRC shrinks with increasing pulse duration, Demir and co-workers (Demir et al., 1997Go) showed, using a model neuron, that for pulses of sufficient amplitude and duration, the limit cycle trajectory actually converges to a rest state, or hyperpolarized fixed point, from which it recovers quickly when released from inhibition. Prinz and co-workers (Prinz et al., 2003Go) have also reported that in conductance-based models of bursting and spiking oscillators, PRCs saturate as synaptic input conductance is increased. This saturation occurs even before synaptic reversal is reached. This would result in a linear PRC for very long pulses, i.e., Fig. 5 C and D. Currently, our method predicts a linear PRC for pulses that arrive on the hyperpolarized branch and are longer than the recovery time on that branch: {Delta}{varphi}T = {varphi}s + {tau}/Pi - 1. The fact that pulses longer than the duration of the silent period do tend to produce linear PRCs (equivalent to a flat phase plot) could be an indication of the existence of a rest state and the only invariant set in the phase space. The waiting time concept {Delta}{varphi}w is an intuitive way to take into account the phase space dynamics toward the invariant sets (fix point, limit cycle, etc.) without assuming a priori knowledge about their topological nature. The computational results of Demir and co-workers (Demir et al., 1997Go) and experimental results of Prinz and co-workers (Prinz et al., 2003Go) offer strong evidence for such an invariant set when the pulse duration exceeds a given threshold, which is, according to our findings, the duration of the interburst interval.

Sources of errors
The major source of errors in predicting phase resetting from a time series is the trial-to-trial variability inherent in the experimental determination of the PRC. In the experimental setup described in this article, the AB/PD complex still receives many modulatory inputs, since only the lateral pyloric was killed. For example, gastric mill events cause a transient change in the period of the AB/PD complex (Mulloney, 1977Go; Bartos and Nusbaum, 1997Go; Marder et al, 1998Go; Nadim et al, 1998Go, 1999Go; Bartos et al, 1999Go; Thuma and Hooper, 2002Go). Other sources of noise besides activity in other neurons are fundamentally random, such as channel noise, shot noise, and noisy release of transmitter. Another potential source of error is in the reconstruction of the limit cycle. Although we are not guaranteed a correct mapping between the points on the limit cycle with respect to the assumed second variable, the results are robust whether the slow variable is a low-pass filtered version of the membrane potential or the envelope, and for somewhat different values of the time lag used for delay embedding. Finally, we have assumed instantaneous normal resetting to the other side of the limit cycle. In reality, activation and inactivation processes associated with voltage-gated channels all have kinetics associated with them, so the actual path may be curved and/or multidimensional. We expect that this will cause more problems with subthreshold pulses, which were not examined in this study. In this vein, it is somewhat surprising that the very long pulses in this study did not require corrections for their effect on the assumed slow variable.

Novelty of the approach
The method we propose is novel in two salient ways. We have predicted the PRC, not from model equations, but from a time series record of membrane potential. Second, we have examined perturbations in a normal rather than a tangential direction along the limit cycle. Previous methods, both analytic (Ermentrout, 1986Go; Izhikevich, 2000Go; Kopell and Ermentrout, 1988Go) and geometric (Oprisan and Canavier, 2002Go) tried to predict phase resetting from model equations. Both assumed that the coupling is in the form of a current perturbation, and that this perturbation in current can be converted to a perturbation in membrane voltage. The displacement in voltage is then projected back onto the limit cycle (which is assumed to have a constant slope for the duration of the perturbation), to determine the phase resetting. We refer to this type of resetting as tangential to the limit cycle. In addition to Winfree's classification of the phase resetting in the type 0 and type 1 (Winfree, 1980Go), there is another classification in terms of Type I and Type II that is correlated to the oscillatory mechanism that determines the excitability of the neuron. The tangential resetting approach works well (Oprisan and Canavier, 2001Go) for Type I oscillators, which have an onset of oscillation at arbitrarily low period due to an underlying saddle-node bifurcation, because they function as integrators (Izhikevich, 2000Go), so the effect of a perturbation is converted to a change in the velocity around the limit cycle. On the other hand, the tangential approach does not work well (Oprisan and Canavier, 2001Go) for Type II oscillators, which have an onset of oscillation at a finite period due to an underlying Hopf bifurcation, because they function as resonators, so the effect of a perturbation in current is to change the shape of the limit cycle with minimal effect on velocity.

Impact on circuit analysis
The method in this article suggests a simple rule for adding hard inhibitory perturbations of the type encountered in this article. The amplitude of the synaptic conductance was not a determinant of the phase resetting, because all pulses received during a burst were above the threshold required to terminate the burst and move the trajectory to the other side of the limit cycle. Thus the effect of an overlapping second perturbation with same timing and duration, of arbitrary amplitude, would be the same as a single pulse alone (Prinz et al., 2003Go). If the pulses are of different durations, then the longer duration alone would determine the phase resetting of the combined pulse. The effect of two pulses whose durations did not overlap would be a simple sum of the effects of the two single pulses. A neuromodulatory substance could affect phase resetting by changing the shape of the limit cycle, most effectively by changing the duty cycle which would change the mapping from the depolarized branch (bursting) to the hyperpolarizing branch (interburst). The understanding of phase resetting gained in this study will help us analyze more complex circuits than was previously possible.


   APPENDIX
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION AND CONCLUSIONS
 APPENDIX
ACKNOWLEDGEMENTS
 REFERENCES
 
Phase space reconstruction using delay embedding
There are two equivalent methods that allow phase space reconstruction based on a time series: delay embedding and derivative embedding method (Abarbanel et al., 1993Go; Takens, 1981Go). To examine a single-valued time series in its true multidimensional form, we use a process called delay coordinate embedding to create a trajectory in the n-dimensional space. We can generate n-dimensional points as follows:

where {tau} is called the embedding lag, and D is the embedding dimension (Packard et al., 1980Go). Takens' theorem states that such a reconstruction exists, and the theorem applies for almost all delays, as long as an infinitely long series of noiseless observations is used.

To estimate the time delay, {tau}, we used the following methods (Abarbanel et al., 1993Go; Broomhead and King, 1986Go; Eckmann and Ruelle, 1985Go; Grassberger and Procaccia, 1983Go; Schouten et al., 1994Go; Wolf et al., 1985Go):

  1. Autocorrelation function, which measures the linear dependence of the time series values at one time on the values at another time. The autocorrelation function may be used to detect deterministic components masked in a random background because autocorrelation functions of deterministic data persist over all time displacements, while autocorrelation functions of stochastic processes tend to zero for large time displacement.
  2. Average mutual information function, which is defined as the average information about x[i] gained when observing x[i + {tau}]. The mutual information can be interpreted as the difference between the uncertainty of x[i] and the remaining uncertainty of x[i] after observing x[i + {tau}]. In other words, it is the reduction in uncertainty of x[i] gained by observing x[i + {tau}]. The time delay must be a large enough that independent information about the system is contained in each component of the reconstructed vector. However, it must not be so large that the components of the reconstructed vector are independent with respect to each other. Conversely, if the time delay is too short, the vector components will not be independent enough and will not contain any new information (Fraser and Swinney, 1986Go).

To estimate the embedded dimension, D, we used the following methods:

  1. False nearest neighbors. Suppose that in an m-dimensional delay space the reconstructed attractor is a one-to-one image of the attractor in the original phase space. In this reconstructed space, the topological properties are preserved. Thus, the neighbors of a given point are mapped onto neighbors in the delay space. Due to the assumed smoothness of the dynamics, neighborhoods of the points are mapped onto neighborhoods again. Suppose an embedding in a k-dimensional space with k < m. Due to this projection, the topological structure is no longer preserved. Points are projected into neighborhoods of other points to which they would not belong in higher dimensions. These points are called false neighbors. The dimension should be increased until no further reduction in false nearest neighbors can be achieved.
  2. The spectrum of the fractal dimensions, which describes the number of variables necessary to define the system, as does the nonfractal dimension of a simple object. From the spectrum of fractal dimensions we evaluated only: a) The capacity dimension, which lacks information on frequency of orbit visits to each square or box and is only a geometric measure of complexity. b) The information dimension, which quantifies the non-uniformity of points distribution on an attractor by weighing the boxes proportional to the number of data points contained in each. If all data points were distributed uniformly among boxes, the information and capacity dimensions would be equal. c) The correlation dimension, which is based on pairwise distances. One implementation calculates distances between every pair of data points to determine the number of pairs less than a distance r. A second implementation constructs spheres of radius r at each point and counts the number of points in each sphere. The global embedding dimension, D, is the minimum number of time-delay coordinates needed so that the trajectories x(t) do not intersect in a D-dimensional space. In dimensions less than D, trajectories can intersect because their projected down into too few dimensions. Subsequent calculations, such as predictions, may then be corrupted. If it is too large, noise and other contamination may corrupt other calculations because noise fills any dimension (Casdagli et al., 1991Go; Sauer et al., 1991Go; Wolf at al., 1985Go).

We used the data from the experiments on the first and second neurons discussed in the Results section. In both cases, the sampling time for the experimental data was 0.5 ms and average period is ~1500 ms. Using the TISEAN package (Hegger et al., 1999Go) we obtained, for the first set of experimental data, the time lag {Delta}t1 = 680 ({approx}340 ms) and for the second set, {Delta}t2 = 600 time steps. For both cases the estimated embedding dimension was D = 3. A reconstructed phase space (Fig. 6 A) displays unrealistically fast oscillations of all state variables (shown only for the first data set). This is an artifact of the delay embedding method that the rapid changes corresponding to an action potential are reflected in every state variable, and every state variable has the same time scale. This leads to a major drawback of the delay embedding method which is its inability to capture the dynamics of phase space. For example, a relaxation oscillator, that frequently characterizes bursts, will always map into a phase oscillator after reconstruction. This fact is particularly important when we are interested not only in capturing the geometry but also the dynamics of the phase space. One possible way to recover a slow variable is to use the envelope of the membrane potential record (Fig. 2 A) or a low-pass filter (Fig. 6 B). The average topology of the reconstructed attractor from the PD/AB oscillator does not change (Fig. 6 B) and reveals the fact that the attractor is relatively flat with only the excursions caused by action potentials resulting in movement of the plane of the envelope of the membrane potential. This means that, despite the fact that the dynamics is described by three independent variables, we can use only two normal coordinates in the plane of the limit cycle and the third variable normal to the limit cycle.



View larger version (31K):
[in this window]
[in a new window]
  		FIGURE 6  Three-dimensional limit cycle unfolding. The membrane potential record from a PD neuron with a sampling time of 0.5 ms and intrinsic period of ~1500 ms is used to unfold the limit cycle. The TISEAN package predicted an optimum time lag of 680 time units (340 ms) and an embedding dimension D = 3. The unfiltered phase space attractor (A) displays spurious spiking phenomenon due to reconstruction artifacts. A low-pass filter indicates that the three-dimensional attractor is "flat," suggesting that only two variables could be used to describe the system (B).

 
For this purpose, our fast variable is the membrane potential record (Fig. 6 A) and the slow variable is its envelope (Fig. 2 A). The phase space reconstruction based on a delay embedding maps the limit cycle attractor into a phase oscillator.


   ACKNOWLEDGEMENTS
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION AND CONCLUSIONS
 APPENDIX
 ACKNOWLEDGEMENTS
REFERENCES
 
The authors thank Dr. Eve Marder for her support and for her comments on an earlier draft of the manuscript.

This research was supported by National Institutes of Health grant NS-17813 (to E.M.) and National Science Foundation grant IBN-0118117 (to C.C.C.).

Submitted on September 20, 2002; accepted for publication January 21, 2003.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION AND CONCLUSIONS
 APPENDIX
 ACKNOWLEDGEMENTS
 REFERENCES
 
Abarbanel, H. D. I., R. Brown, J. J. Sidorowich, and L. S. Tsimring. 1993. The analysis of observed chaotic data in physical systems. Rev. Mod. Phys. 65:1331–1392.

Abbott, L. F., E. Marder, and S. L. Hooper. 1991. Oscillating network: control of burst duration by electrically coupled neurons. Neural Comput. 3:487–497.

Abramovich-Sivan, S., and S. Akselrod. 1998a. A PRC-based model of a pacemaker cell: effect of vagal activity and investigation of the respiratory sinus arrhythmia. J. Theor. Biol. 92:219–243.

Abramovich-Sivan, S., and S. Akselrod. 1998b. A phase response curve based model: effect of vagal and sympathetic stimulation and interaction on a pacemaker cell. J. Theor. Biol. 192:567–579.[Medline]

Arshavsky, Y. I., I. N. Beloozerova, G. N. Orlovsky, Y. V. Panchin, and G. A. Pavlova. 1985. Control of locomotion in marine mollusc Clione limacina. I. Efferent activity during actual and fictitious swimming. Exp. Brain Res. 58:255–262.[Medline]

Bal, T., F. Nagy, and M. Moulins. 1988. The pyloric central pattern generator in crustacean: a set of conditional neuronal oscillators. J. Comp. Physiol. 163:715–727.

Bartos, M., and M. P. Nusbaum. 1997. Intercircuit control of motor pattern modulation by presynaptic inhibition. J. Neurosci. 17:2247–2256.[Abstract/Free Full Text]

Bartos, M., Y. Manor, F. Nadim, E. Marder, and M. P. Nusbaum. 1999. Coordination of fast and slow rhythmic neuronal circuits. J. Neurosci. 19:6650–6660.[Abstract/Free Full Text]

Baxter, R. J., C. C. Canavier, H. A. Lechner, R. J. Butera, A. A. DeFranceschi, J. W. Clark, and J. H. Byrne. 2000. Coexisting stable oscillatory states in single cell and multicellular neuronal oscillators. In Oscillations in Neural Systems. D. Levine, V. Brown, and T. Shirley, editors. Erlbaum Associates, New Jersey. pp. 51–77.

Beer, R. D., H. J. Chiel, and J. C. Gallagher. 1999. Evolution and analysis of model CPGs for walking. II. General principles and individual variability. J. Comp. Neurosci. 7:119–147.[Medline]

Broomhead, D. S., and G. P. King. 1986. Extracting qualitative dynamics from experimental data. Physica D. 20:217–236.

Butera, R. J., J. W. Clark, D. A. Baxter, C. C. Canavier, and J. H. Byrne. 1995. Analysis of the effects of modulatory agents on a modeled bursting neuron: dynamic interactions between voltage and calcium-dependent mechanisms. J. Comput. Neurosci. 2:19–44.[Medline]

Canavier, C. C., R. J. Butera, R. O. Dror, D. A. Baxter, J. W. Clark, and J. H. Byrne. 1997. Phase response characteristics of model neurons determine which patterns are expressed in a ring circuit model of gait generator. Biol. Cybern. 77:367–380.[Medline]

Canavier, C. C., J. W. Clark, and J. H. Byrne. 1991. Simulation of the bursting activity of neuron R15 in aplysia: role of ionic currents, calcium balance, and modulatory transmitters. J. Neurophysiol. 66:2107–2124.[Abstract/Free Full Text]

Canavier, C. C., D. A. Baxter, J. W. Clark, and J. H. Byrne. 1999. Control of multistability in ring circuits of oscillators. Biol. Cybern. 80:87–102.

Casdagli M., E. Stephen, J. D. Farmer, and J. Gibson. 1991. State space reconstruction in the presence of noise. Physica D. 51:52–98.

Chay, T. R., and J. E. Keizer. 1983. Minimal model for membrane oscillations in the pancreatic cell. Biophys. J. 42:181–190.[Abstract/Free Full Text]

Chiel, H. J., R. D. Beer, and J. C. Gallager. 1999. Evaluation and analysis of model CPGs for walking. I. Dynamical models. J. Comp. Neurosci. 7:1–20.

Collins, J. J., and I. N. Stewart. 1993. Coupled nonlinear oscillators and the symmetries of animal gaits. J. Nonlin. Sci. 3:349–392.

Collins, J. J., and S. A. Richmond. 1994. Hard-wired central pattern generators for quadrupedal locomotion. Biol. Cybern. 71:375–385.

Demir, S. S., R. J. Butera, A. A. DeFranceshi, J. W. Clark, and J. H. Byrne. 1997. Phase sensitivity and entrainment in a modeled bursting neuron. Biophys. J. 72:579–594.[Medline]

Dror, R. O., C. C. Canavier, R. J. Butera, J. W. Clark, and J. H. Byrne. 1999. A mathematical criterion based on phase response curves for the stability of a ring network of oscillators. Biol. Cybern. 80:11–23.

Eckmann, J. P., and D. Ruelle. 1985. Ergodic theory of chaos and strange attractors. Rev. Mod. Phys. 57:617–620.

Ermentrout, G. B. 1985. The behavior of rings of coupled oscillators. J. Math. Biol. 23:55–74.[Medline]

Ermentrout, G. B. 1986. Parabolic bursting in an excitable system coupled with a slow oscillation. SIAM J. Appl. Math. 46:233–253.

Ermentrout, G. B., and N. Kopell. 1991. Oscillator death in systems of coupled neural oscillators. SIAM J. Appl. Math. 29:195–217.

Ermentrout, G. B. 1996. Type I membranes, phase resetting curves, and synchrony. Neural Comput. 8:979–1001.[Abstract]

Fraser, A. M., and H. L. Swinney. 1986. Independent coordinates for strange attractors. Phys. Rev. 33A:1134–1140.

Golubitsky, M., I. Stewart, P.-L. Buono, and J. J. Collins. 1998. A modular network for legged locomotion. Physica D. 115:56–72.

Grassberger, P., and I. Procaccia. 1983. Estimation of Kolmogorov entropy from a chaotic signal. Phys. Rev. A. 28:2591–2599.

Guckenheimer, J., and P. Holmes. 1983. Nonlinear Oscillations, Dynamical Systems, and Bifurcation of the Vector Fields. Springer-Verlag, Berlin, Heidelberg, New York.

Hegger, R., H. Kantz, and T. Schreiber. 1999. Practical implementation of nonlinear time series methods: the TISEAN package. Chaos. 9:413–435.[Medline]

Hooper, S., and E. Marder. 1987. Modulation of the lobster pyloric rhythm by the peptide proctolin. J. Neurosci. 7:2097–2112.[Abstract]

Hoppensteadt, F. C., and E. M. Izhikevich. 1997. Weakly Connected Neural Networks. Springer-Verlag, New York.

Izhikevich, E. M. 2000. Neural excitability, spiking, and bursting. Int. J. Bif. Chaos. 10:1171–1266.

Kopell, N., and G. B. Ermentrout. 1988. Coupled oscillators and the design of central pattern generators. Math. Biol. 90:87–109.

Liu, Z., J. Golowasch, E. Marder, and L. F. Abbott. 1998. A model neuron with activity-dependent conductances regulated by multiple calcium sensors. J. Neurosci. 18:2309–2320.[Abstract/Free Full Text]

Marder, E., and P. Meyrand. 1989. Chemical modulation of an oscillatory neural circuit. In Neural and Cellular Oscillators. J.W. Jacklet, editor. Marcel Dekker, New York. 317–338.

Marder, E., and R. L. Calabrese. 1996. Principles of rhythmic motor pattern generation. Physiol. Rev. 76:687–717.[Abstract/Free Full Text]

Marder, E., Y. Manor, F. Nadim, M. Bartos, and M. P. Nusbaum. 1998. Frequency control of a slow oscillatory network by a fast rhythmic input: pyloric to gastric mill interactions in the crab stomatogastric nervous system. Ann. N. Y. Acad. Sci. 860:226–238.[Abstract/Free Full Text]

Miller, J. P., and A. I. Selverston. 1979. Rapid killing of single neurons by irradiation of intracellularly injected dye. Science. 206:702–704.[Abstract/Free Full Text]

Miller, J. P., and A. I. Selverston. 1982. Mechanism underlying pattern generation in lobster stomatogastric ganglion as determined by selective inactivation of identified neurons. II. Oscillatory properties of pyloric neurons. J. Neurophysiol. 48:1378–1391.[Abstract/Free Full Text]

Mulloney, B. 1977. Organization of the stomatogastric ganglion of the spiny lobster. V. Coordination of the gastric and pyloric systems. J. Comp. Physiol. 122:227–240.

Murray, J. D. 1993. Mathematical Biology, 2nd ed. Springer-Verlag, New York.

Nadim, F., Y. Manor, M. P. Nusbaum, and E. Marder. 1998. Frequency regulation of a slow rhythm by a fast periodic input. J. Neurosci. 18:5053–5067.[Abstract/Free Full Text]

Nadim, F., Y. Manor, N. Kopell, and E. Marder. 1999. Synaptic depression creates a switch that controls the frequency of an oscillatory circuit. Proc. Natl. Acad. Sci. USA. 96:8206–8211.[Abstract/Free Full Text]

Oprisan, S. A., and C. C. Canavier. 2001. Stability analysis of rings of pulse-coupled oscillators: the effect of phase resetting in the second cycle after the pulse is important at synchrony and for long pulses. Diff. Eqs. Dyn. Sys. 9:243–258.

Oprisan, S. A., and C. C. Canavier. 2002. The influence of limit cycle topology on the phase resetting curve. Neural Comput. 14:1027–1057.[Abstract/Free Full Text]

Packard, N. H., J. P. Crutchfield, J. D. Farmer, and R. S. Shaw. 1980. Geometry from a time series. Phys. Rev. Lett. 45:712–716.

Pavlides, T. 1973. Biological Oscillators: Their Mathematical Analysis. Academic Press, New York.

Pearson, K. 1993. Common principles of motor control in vertebrates and invertebrates. Annu. Rev. Neurosci. 16:256–297.

Perkel, D. H., J. H. Schulman, T. H. Bullock, G. P. Moore, and J. P. Segundo. 1964. Pacemaker neurons: effects of regularly spaced synaptic input. Science. 145:61–63.[Abstract/Free Full Text]

Pinsker, H. M. 1977. Aplysia bursting neurons as endogenous oscillators. I. Phase-response curves for pulsed inhibitory synaptic input. J. Neurophysiol. 40:527–543.[Abstract/Free Full Text]

Pinto, R. D., R. C. Elson, A. Szücs, M. I. Rabinovitch, A. I. Selverston, and H. D. I. Abarbanel. 2001. Extended dynamic clamp: controlling up to four neurons using a single desktop computer and interface. J. Neurosci. Meth. 108:39–48.[Medline]

Prinz, A. A., V. Thirumalai, and E. Marder. 2003. The functional consequences of changes in the strength and duration of synaptic input to oscillatory neurons. J. Neurosci. 23:943–954.[Abstract/Free Full Text]

Richards, K. S., W. L. Miller, and E. Marder. 1999. Maturation of lobster stomatogastric ganglion rhythmic activity. J. Neurophysiol. 82:2006–2009.[Abstract/Free Full Text]

Rinzel, J., and Y. S. Lee. 1986. On different mechanisms for membrane potential bursting, In Nonlinear Oscillations in Biology and Chemistry. H. G. Othmer, editor. Springer-Verlag, New York. pp.19–33.

Rinzel, J., and G. B. Ermentrout. 1998. Analysis of neural excitability and oscillations. In Methods in Neuronal Modeling: From Ions to Networks. C. Koch, and I. Segev, editors. MIT Press, Cambridge. pp.135–169.

Robertson, R. M., and K. G. Pearson. 1985. Neural circuits in the flight system of locusts. J. Neurophysiol. 53:110–128.[Abstract/Free Full Text]

Satterlie, R. A. 1989. Reciprocal inhibition and rhythmicity: swimming in a pteropod mollusk. In Neural and Cellular Oscillators. J. W. Jacklet, editor. Marcel Dekker, New York. pp.151–171.

Sauer, T., J. A. Yorke, and M. Casdagli. 1991. Embedology. J. Stat. Phys. 65:579–616.

Schouten, J. C., F. Takens, and C. M. van den Bleek. 1994. Estimating the dimension of a noisy attractor. Phys. Rev. E. 50:1851–1861.

Selverston, A. I. and M. Moulins, editors. 1987. The crustacean stomatogastric system. Springer-Verlag, Berlin.

Sharp, A. A., M. B. O'Neil, L. F. Abbott, and E. Marder. 1993. Dynamic clamp: computer-generated conductances in real neurons. J. Neurophysiol. 69:992–995.[Abstract/Free Full Text]

Sharp, A. A., F. K. Skinner, and E. Marder. 1996. Mechanisms of oscillation in dynamic-clamp-constructed two-cell half-center circuits. J. Neurophysiol. 76:867–883.[Abstract/Free Full Text]

Takens, F. 1981. Detecting strange attractors in turbulence, dynamical systems and turbulence. In Lecture Notes in Mathematics No. 898. D. Rand, and L.-S. Young, editors. Springer, Berlin. pp.366–381.

Thuma, J. B., and S. L. Hooper. 2002. Quantification of gastric mill network effects on a movement related parameter of pyloric network output in the lobster. J. Neurophysiol. 87:2372–2384.[Abstract/Free Full Text]

Winfree, A. T. 1980. The geometry of biological time. Springer-Verlag, New York.

Wolf, A., J. B. Swift, H. L. Swinney, and J. A. Vastano. 1985. Determining Lyapunov exponents from a time series. Physica D. 16:285–317.




This article has been cited by other articles:


Home page
 J. Neurophysiol.Home page
V. Thirumalai, A. A. Prinz, C. D. Johnson, and E. Marder
Red Pigment Concentrating Hormone Strongly Enhances the Strength of the Feedback to the Pyloric Rhythm Oscillator But Has Little Effect on Pyloric Rhythm Period
J Neurophysiol, March 1, 2006; 95(3): 1762 - 1770.
[Abstract] [Full Text] [PDF]

Home page
 Proc. Natl. Acad. Sci. USAHome page
R. A. Gray and N. Chattipakorn
Termination of spiral waves during cardiac fibrillation via shock-induced phase resetting
PNAS, March 29, 2005; 102(13): 4672 - 4677.
[Abstract] [Full Text] [PDF]

Home page
 J. Neurophysiol.Home page
T. I. Netoff, M. I. Banks, A. D. Dorval, C. D. Acker, J. S. Haas, N. Kopell, and J. A. White
Synchronization in Hybrid Neuronal Networks of the Hippocampal Formation
J Neurophysiol, March 1, 2005; 93(3): 1197 - 1208.
[Abstract] [Full Text] [PDF]