Phase Resetting and Phase Locking in Hybrid Circuits of One Model and One Biological Neuron

Electrophysiology

H. americanus were purchased from Yankee Lobster (Boston, MA) and maintained in artificial seawater at 10–12°C until used. The stomatogastric nervous system (STNS) was dissected out and pinned out in a dish coated with Sylgard (Dow Corning, Midland, MI), and the STG was desheathed with fine forceps. Throughout the experiments, the stomatogastric nervous system was superfused with chilled (9–14°C) saline containing (in mM) NaCl, 479.12; KCl, 12.74; CaCl₂, 13.67; MgSO₄, 20; Na₂SO₄, 3.91; HEPES, 5; pH 7.45. Extracellular recordings were made with stainless steel pin electrodes in Vaseline wells on the motor nerves and amplified with a differential AC amplifier (Model 1700, A-M Systems, Carlsborg, WA). Intracellular recordings from cells in the STG were obtained with an Axoclamp 2B amplifier (Axon Instruments, Foster City, CA) in discontinuous current clamp mode using microelectrodes filled with 0.6 M K₂SO₄ and 20 mM KCl; electrode resistances were in the range of 20–40 MΩ. Extracellular and intracellular potential traces were digitized with a Digidata 1200A board (Axon Instruments), recorded using Clampex 8.0 software (Axon Instruments), and analyzed with in-house software. The PD and lateral pyloric (LP) motor neurons were identified based on their membrane potential waveforms, the timing of their activity in the pyloric rhythm, and their axonal projections to the appropriate motor nerves. The only synaptic feedback to the pyloric pacemaker group through the LP to PD inhibitory synapse was removed by applying 10⁻⁵ M picrotoxin in the bath. The pharmacologically isolated pyloric pacemaker was monitored by impaling one of the PD neurons; it served as the biological oscillator in the experiments reported here.

Dynamic clamp

We recorded the membrane potential of the AB/PD complex and used the dynamic clamp (Sharp et al., 1993a,b) to replace the synaptic input from LP onto PD with artificial synaptic inputs: the membrane potential V at the PD cell body was amplified, fed into a Digidata 1200A board (Axon Instruments), and digitized at a rate of 1.7 kHz with in-house software modified from a C++ program kindly provided by R. Pinto (Pinto et al., 2001). The dynamic clamp program detected bursts in the ongoing PD rhythm and monitored the instantaneous period. Artificial dynamic clamp synaptic inputs were generated at different phases of the PD rhythm by setting the synaptic activation to the desired value for the desired duration. During the synaptic input, the program computed the momentary synaptic current as described below. To inject this synaptic current into the PD neuron, the program computed the corresponding command voltage, which was turned into an analog voltage by the Digidata board and sent to the electrode amplifier.

Model neurons

The model neurons had a single compartment with eight Hodgkin-Huxley type membrane currents and an intracellular calcium buffer. The membrane currents were based on voltage-clamp experiments on lobster stomatogastric neurons (Turrigiano et al., 1995) and included a fast sodium current (I_Na), a fast and a slow transient calcium current (I_CaT and I_CaS), a fast transient potassium current (I_A), a calcium-dependent potassium current (I_KCa), a delayed rectifier potassium current (I_Kd), and a leak current (I_leak). The model neurons used here had identical current dynamics and differed only in the maximal conductances of their eight membrane currents; these conductances were chosen to produce different burst periods, durations, and duty cycles in the different model neurons. The maximal conductances of the eight currents and the corresponding intrinsic periods and burst durations for all seven model neurons considered in this study are listed in Table 1 in the Supplementary Material. The model was described in detail in Prinz et al. (2003), and similar models have been used before (Goldman et al., 2001; Liu et al., 1998). The model was implemented in a C++ program, and all differential equations were integrated with an exponential Euler method at a time resolution of ∼50 μs, corresponding to 10 updates of the model neuron for every voltage value read and current command written by the dynamic clamp.

Artificial synapses

For all synapses simulated here, the synaptic current was I_syn = g_synm (V_postsyn − E_rev) where g_syn is the conductance of the synapse, m is the synaptic activation, V_postsyn is the membrane potential of the postsynaptic neuron, and E_rev is the synaptic reversal potential, which was set to −90 mV for all synapses based on voltage clamp experiments on the LP to PD synapse (Thirumalai, 2002). All synapses were instantaneously activating and deactivating. The synaptic activation function, , was chosen such that the synapse is fully activated during the burst and deactivated between bursts.

Because the potential between spikes reaches very hyperpolarized levels, and the activation of the synapses is instantaneous, we used a filtered version, V_filt, of the presynaptic membrane potential to compute the synapse activation to smooth it. For every time step Δt, V_filt was updated according to V_filt(t + Δt) = V_filt(t) + (V_filt(t + Δt) − V_filt(t))Δt/τ_filt, where τ_filt = 1 ms for the biological neuron. The value of V_1/2 = −45 mV for synapses from the biological neurons was set below the minimal filtered voltage between spikes within a burst. We chose s = 0.1 mV for all synapses, which ensured rapid activation and deactivation of the synapse during the up- and downstrokes of each presynaptic burst. For synapses from model neurons, the values of V_1/2, and the corresponding filtering constants, τ_filt, were chosen such that the synapse remained fully activated throughout the burst, even during the hyperpolarized intervals between spikes. The chosen values ranged from −47 mV to −40 mV for V_1/2 and from 10 mS to 30 ms for τ_filt (Table 2 in the Supplementary Material).

PRCs for biological and model neurons

To generate the PRC of a biological neuron in response to synaptic input from a given model neuron and for a given synapse strength g_model>bio, we first determined the intrinsic period P₀ and burst duration b₀ of the biological neuron from the filtered membrane potential recording of the unperturbed neuron of at least 20 s. We define the intrinsic period P₀ of the biological neuron as the time between two successive crossings of the −45 mV voltage threshold with positive slope and the intrinsic burst duration b₀ as the time between the voltage threshold crossing with positive slope and the following crossing with negative slope (Fig. 1 C). The intrinsic periods of the four biological neurons immediately before different couplings to model neurons ranged from 1.12 s to 2.5 s and their burst durations were between 0.44 s and 0.88 s (Table 3 in the Supplementary Material). Once the intrinsic period was determined, the membrane potential of the biological neuron in response to conductance pulses of amplitude g_model>bio was recorded, using a stored activation profile of the artificial synapse, which corresponds to the activation produced by a single presynaptic burst in a model neuron in the open loop condition (no feedback). Individual stimuli were delivered ∼10 s apart to ensure that the biological neuron returned to its unperturbed activity between stimuli. The stimulus interval was computed by multiplying the desired stimulus phase by the period of the preceding cycle. For a full PRC, we delivered stimuli at 20 (for biological neuron 1) or 40 (for the other biological neurons) equally spaced phases between 0 and 1. Each PRC was recorded a total of four times. The experimental protocol described above was repeated for the same model neuron for all synaptic conductances g_model>bio listed in Table 1 (third column). Once all synaptic conductances g_model>bio were completed, the next model neuron was selected.

A similar experimental protocol was used to extract the PRCs for the model neurons at 40 equally spaced phases between 0 and 1. To obtain the PRCs of the bursting model neurons, we applied a synaptic input with a stored activation profile of the artificial synapse, which in this case corresponds to the activation produced by a single presynaptic burst in the biological neuron in the open loop condition (no feedback), at different times during the ongoing simulated rhythm. We assigned phase zero to the burst onset, which was specified by a voltage threshold of −42 mV (Fig. 1 D). Each PRC was recorded only once.

Dynamic-clamp generated hybrid circuits

We used the dynamic clamp to generate 164 different mutually inhibitory hybrid circuits between one of four biological neurons and one of seven model neurons (Fig. 1 B). The range of synaptic conductances used in our experiments encompasses the estimated physiological range of synaptic conductances onto the AB/PD complex in the STG (between 10 nS and 60 nS according to Thirumalai, 2002) and extends beyond the 500-nS limit considered in the existing implementations of analog two-neuron hybrid networks (Pinto et al., 2000). For each circuit, we recorded the membrane potential of the biological and the model neuron for at least 20 s before and 40 s after the synapses were switched on. The recording in the closed loop configuration was used to extract the experimental values of the relative burst timing of the biological and the model neuron as well as the phase-locked period.

Theoretical method

The theoretical method for predicting phase-locked modes has two components, existence and stability. Only modes that satisfy the periodicity constraints can exist, hence the existence criterion, and only those that are stable can be observed in the presence of the small perturbations that result from the unavoidable presence of noise, hence the stability criterion. The periodicity constraints (Canavier et al., 1997, 1999; Luo et al., 2004) include the assumptions that the steady-state phase-locked period of each neuron is determined by the closed loop phase resetting produced by the input from its partner (right-hand side of Eq. 1) and that this resetting is equal to that observed in the open loop condition. In addition, there is the obvious requirement that in 1:1 phase locking the period of each neuron must be the same and equal to the sum of the stimulus interval (t_s) and the recovery interval (t_r) for each neuron. The equivalence of the sum of the stimulus and recovery intervals with the period in the closed circuit can be stated as

(1)

,where the index j stands for bio or model neuron, and ϕ_bio and ϕ_model, represent the phase of each neuron, corresponding to the position on the appropriate unperturbed limit cycle at which an input is received. The stimulus interval (t_s,j) is defined as the interval between the initiation of a burst in neuron j and the beginning of an input onto neuron j, and the recovery interval is defined as the interval between the beginning of an input onto neuron j and the initiation of the next burst in neuron j (Fig. 2 A). The first and second order PRCs can be used to predict the steady-state values of t_s,j and t_r,j as follows:

(2)

.Since both of the above quantities depend solely upon the phase at which an input is received, one can plot the pair of points (t_r,model, t_s,model) for each value of ϕ_model as well the pair of points (t_s,bio, t_r,bio) for each value of ϕ_bio, discarding any resultant negative values for the intervals because they violate causality. Due to the choice of axes for these plots, the intersection(s) of these two curves at the periodicity criterion t_r,bio = t_s,model and t_s,bio= t_r,model will give the values of the stimulus and recovery times for which the periodicity constraints in Eq. 1 are satisfied for both values of j: model and bio. If there is a fixed time difference between the initiation of a presynaptic burst and the activation of the synaptic coupling to the postsynaptic neuron, then the receipt of the input by the postsynaptic neuron does not coincide with the initiation of the presynaptic burst. The stimulus and recovery intervals relative to burst onset are shown in Fig. 2 B for a particular case in which the coupling is activated before the first spike in the burst occurs, which could occur in the case of graded synaptic coupling. The periodicity criterion is now t_s,model +δ = t_r,bio and t_s,bio +δ = t_r,model, so one must plot the pairs (t_r,model, t_s,model + δ) and (t_s,bio + δ, t_r,bio) to obtain the intersection that satisfies this periodicity criterion.

The stability criteria were initially formulated by ignoring the contribution of second order resetting and using the first order PRC to define a mapping between the phases at which inputs are received in successive cycles in terms of the dependence of φ_j[n] on φ_j[n − 1] where n indicates the current cycle and n − 1 indicates the previous cycle. A mapping in terms of the stimulus and recovery intervals t_s,j[n] and t_r,j[n] is shown in Fig. 2 C, and due to the implicit dependence of these quantities on φ_j[n], the mapping can be written in terms of the φ_j[n]. If a steady state is assumed for both φ_model[∞] and φ_bio[∞], the mapping also gives the deviation of each phase, Δφ_j[n], from its steady-state value. This mapping can be linearized by assuming that the change in cycle period that occurs in cycle n due to a small perturbation in the previous cycle, Δφ_j[n − 1], can be given by m_1,jΔφ_j[n − 1] where m_1,j is the slope of the first order PRC F′_1,j(φ_j[∞]) at the phase at which neuron j receives an input in the assumed steady 1:1 phase-locked mode. This leads immediately to the result of Dror et al. (1999) that Δφ_j[n] = (1 − m_1,bio)(1 − m_1,model)Δφ_j[n − 1]. The quantity (1 − m_1,bio)(1 − m_1,model) is a multiplier λ that determines whether a perturbation from the assumed steady 1:1 phase-locked mode will increase or decrease. If −1 < λ < 1, the perturbation will decrease to zero, guaranteeing the stability of the assumed 1:1 phase-locked mode. Thus the stability criterion if F_2,j(φ_j) = 0 is −1 < (1 − m_1,bio)(1 − m_1,model) < 1. However, if F_2,j(φ_j) ≠ 0, then two previous cycles (n − 1 and n − 2) must be taken into account. Using the same methodology, and defining m_2,j as the slope of the second order PRC F′_2,j(φ_j[∞]) at the phase that neuron j receives an input in the assumed steady 1:1 phase-locked mode, the future of a perturbation is given by a higher order discrete map because an additional cycle must be taken into account (Oprisan and Canavier, 2001):

(3)

.To guarantee stability, meaning that perturbations die out and go to zero, the solutions to the characteristic equation of Eq. 3, i.e.,

(4)

,which are λ₁ and λ₂, both must be in the range from −1 to 1 for the presumed 1:1 phase-locked mode to be stable. Note that by neglecting the contribution of the second order PRC in Eq. 4, i.e., setting both m_2,model and m_2,bio equal to zero, we recover the result of Dror et al. (1999) that λ = (1 − m_1,bio)(1 − m_1,model).

Phase resetting curves

We extracted all PRCs from the membrane potential recordings offline, using a burst threshold of −45 mV for all biological neurons. For each synaptic coupling g_model>bio, the average of the four PRCs was computed and approximated with a polynomial fit. The upper and lower envelopes reflect the trial-to-trial variability of biological neuron PRCs (dashed lines in Fig. 3, A1 and A2). The upper (lower) envelope is the polynomial fit through the maximal (minimal) values of the phase resetting measured in the four experiments versus the stimulus phase. Typical first order PRCs obtained for the biological neuron in the open loop setup are shown in Fig. 3 A1, and typical second order resetting curves are given in Fig. 3 A2. The results from the four different trials described in the Methods are indicated by different symbols, and the average resetting and both envelopes are shown. In general, the first order resetting consisted of advances early in the cycle associated with burst truncation followed by a shortened cycle, a flat region in the middle of the cycle, and delays near the end of the cycle as the subsequent burst was delayed by the prolonged inhibition. The observed first order resetting is consistent with previous work (Oprisan et al., 2003; Prinz et al., 2003), and because the applied perturbation is essentially a square pulse, the dependency of the PRC in the biological neuron on the particular model neuron used to drive it reduces to the different burst durations of the various models. Again consistent with previous work, the longer pulse durations reduce the extent of the nearly flat region in the middle in favor of a more uniformly linear PRC. The second order PRC for the biological neuron was in general essentially flat and equal to zero. Additional examples of PRCs for biological neurons are given in the Supplementary Material.

First order resetting curves are given for model neuron 1 in Fig. 3 B1 and second order in Fig. 3 B2. Here, only a single trial is shown, but trials at different values of the conductance g_bio>model are indicated by different symbols. The first order resetting is typical of the model neuron PRCs in that at 1 nS for the g_bio>model conductance the PRC has some curvature, but as the conductance is increased to 10 nS and above, the PRC saturates and becomes linear. This is likely due to the inhibition driving the model neuron to a fixed point (Demir et al., 1997) such that the burst follows the end of the perturbation by a fixed interval. As with the biological neuron, advances are observed early in the cycle and delays late in the cycle. The second order resetting in Fig. 3 B1 is typical of model neurons in that it is essentially flat and thus independent of phase with a slope near zero. It is also typical in that the value of the second order resetting is near zero at g_bio>model = 1 nS, but at a higher value (10 nS), a delay is observed largely due to a lengthened burst in the second cycle. The second order resetting shown in Fig. 3 B2 is atypical in that the delay that constitutes the second order resetting usually saturates at higher conductances, but in the case of model 1 at large conductance values (above 100 nS) an advance is observed. This advance is largely due to burst shortening in the subsequent cycle, as observed for the same model neuron driven by a different biological neuron in Fig. 1 D. Additional examples of resetting curves for the model neurons are shown in the Supplementary Material.

Summary of predictions

Our intent was to predict which of the 164 distinct circuits created by utilizing different biological neurons reciprocally coupled to different models at varying values of the synaptic conductances in each direction would exhibit stable 1:1 phase locking, and further to predict the actual values of the stimulus and recovery intervals for each neuron as well as the phase-locked period for the circuit. Fig. 4 illustrates the type of firing pattern that one would expect if no stable 1:1 phase locking existed: the period, stimulus intervals, and even firing order are variable since the model neuron often, but not always, fires twice before the biological neuron fires. Table 1 indicates for each model circuit that 1:1 locking was observed by an “a” before the slash and that it was predicted by an “a” after the slash, such that “a/a” indicates a successful prediction. In this table, “c/−” also indicates a success because a complex mode such as that shown in Fig. 4 was observed, and the dash indicates that no 1:1 mode was predicted to exist. The three failures are indicated in bold by “c/a,” which signifies that a complex mode was observed despite a prediction of 1:1 phase locking. In each of these cases, g_bio>model was less than or equal to 2 nS, so the difficulty seems to be related to the low value of this parameter. In addition, these were the only cases in which the burst in AB/PD would have been truncated in the predicted mode. We will illustrate the prediction method using two successful specific examples and then summarize the accuracy of our predictions quantitatively.

Example of a prediction that did not require consideration of second order resetting

Fig. 5 gives an example of the prediction method using a hybrid circuit with g_model1>bio4 = 100 nS and g_bio4>model1 = 100 nS. The polynomial fit to the average of the first order PRCs (Fig. 5 A1) and to the average total resetting (Fig. 5 A2) were used to compute the dependence of recovery interval on the stimulus interval for both the biological (long dashes) and the model (short dashes) neuron as explained in the Methods, with the contribution of second order resetting ignored in Fig. 5 B1 but considered in Fig. 5 B2. The contribution of second order resetting in the biological neuron is negligible (Fig. 3 A2), but the effect of including the essentially constant second order resetting (triangles pointing up in Fig. 3 B2) of about −0.25 is to shift the linear PRC for the model in a downward direction. It also produces a downward shift of the recovery intervals for the model neuron (short dashes) between Fig. 5, B1 and B2. This shift causes the short dashed graph to appear shortened, because negative (acausal) stimulus intervals were not plotted (see Methods).

The intrinsic period of the biological neuron 4 was 1.29 ± 0.05 s (Fig. 5 C1), and the intrinsic period of model neuron 1 was 1.09 s (Fig. 5 C2), and Fig. 5 D shows the 1:1 phase locking that was actually observed when the two were coupled in the hybrid circuit. The predictions using only first order resetting (the intersection point in Fig. 5 B1) were quite accurate, giving a value of 0.651 s for t_rmodel1, which is a 2% error compared to the average observed value of 0.674 s, and a value of 0.656 s for t_rbio4, which has a 3% error compared to the average observed value of 0.637 s. Adding these two quantities results in a predicted phase-locked period of P_e = 1.307 s, which has only a 0.4% error compared to the average observed value of 1.312 s. Using the second order PRC contributions to the total phase resetting moved the intersection point only slightly to the left in Fig. 5 B2, resulting in a predicted value of 0.659 s for t_rmodel1, which gives a 3% error rather than a 2% error, and a predicted value of 0.658 s for t_rbio4 = 0.658 s (a 2% rather than a 3% error) and an identical prediction for the phase-locked period of 1.317 s. Therefore, for this particular example, the prediction based on both F₁ and F₂ is not significantly different from the prediction based only on the first order resetting, F₁.

The reason for the lack of impact of the contribution of F₂ results from the resemblance of the graph of the stimulus versus recovery interval for the model neuron to vertical lines in Fig. 5, B1 and B2. This is a consequence solely of the linear form (Demir et al., 1997; Oprisan et al., 2003; Prinz et al., 2003) of the first order PRC for the model neuron (Fig. 3 B1). The PRC has essentially unit slope (see m_1,model in Table 2). Therefore, the first order resetting is equal to the phase minus a constant we can denote ϕ₀. Substituting ϕ − ϕ₀ for F₁(ϕ) in the expression for the recovery interval of the model neuron (see Methods) produces t_rmodel = P_0model (1 − ϕ₀), so the x coordinate of the short dashed graph in Fig. 5, B1 and B2, is very nearly constant and independent of the stimulus phase ϕ. On the other hand, the contribution of F₂ shifts the y coordinate, the model stimulus interval, downward by ∼250 ms for each recovery interval. The example in the next section shows that although this shift is quantitatively unimportant for the example in Fig. 5, it can be qualitatively very significant if it causes the intersection to appear or disappear entirely.

A close inspection of the actual membrane potential recordings in the phase-locked mode (Fig. 5 D) shows that although the duration of the burst in the biological neuron remained near its intrinsic value of 0.527 s, the model neuron burst duration was 0.12 s during phase-locked mode compared to an intrinsic value of 0.28. We conclude that near the particular phase 0.523 (marked by squares in Fig. 5, A1 and A2) at which the biological neuron receives an input from the model neuron in the hybrid circuit, the phase resetting is not very sensitive to burst duration, and hence the prediction algorithm is very robust in this case to striking changes in burst duration from the open to the closed loop condition. Since insensitivity to burst duration is not always the case (Oprisan et al., 2003; Prinz et al., 2003), this is a potential source of error (see Discussion).

Fig. 6 summarizes the predictions made using the hybrid circuit with g_model1>bio4 = 100 nS at seven values of g_bio4>model1. A shows the recovery interval for the model neuron, B shows the recovery interval for the biological neuron, and C shows their sum, which is the period. The experimentally measured values are given in all panels by black squares connected by black lines. The predictions using the average resetting are given by red circles connected by a red line, and in this case overlay quite well with the experimental data. The predictions based on the upper and lower (triangles on the blue and green curves) envelopes largely bracket the experimental results. Note that the predicted phase-locked period (Fig. 6 C) is always close to the intrinsic period of the biological neuron (1.29 s) at all values of g_bio4>model1, suggesting that the flat region in the PRC of the biological neuron (Fig. 5, A1 and A2) may be a globally attracting set of the network dynamics.

Thus far, we have illustrated only the portion of the analysis that uses periodicity constraints to predict modes that can exist because they satisfy these constraints. To predict that a mode will be observed, we must also predict that the mode will be stable, and as described in the Methods, the contribution of F₂ to the stability criterion can either be considered or ignored. In general, the hybrid circuit did not provide a good test case to determine whether the stability criteria correctly predict the effect of F₂, because in general, the slopes of the first order PRCs were between 0 and 1, and the slope of the second order PRC was near 0. These values guarantee stability regardless of whether the contribution of F₂ is considered. For the example given in Fig. 5, the entry in Table 2 for g_bio4>model1 = 100 nS shows that ignoring F2 gives a multiplier λ equal to zero, whereas considering F₂ produces λ₁ and λ₂ equal to 0.108 and 0.037, all of which have absolute values smaller than 1. Hence both methods predict a stable 1:1 phase-locked mode, as they do for all seven values of g_bio4>model1 = 100 nS considered in Fig. 6. The very small values of λ imply that perturbations decay quickly, producing rapid convergence to the limit cycle for the component neurons as required by our assumptions.

Example of a prediction that did require consideration of second order resetting

Fig. 7 shows an example of a hybrid circuit in which the consideration of second order resetting was crucial to a correct analysis, the hybrid network formed with the biological neuron 1 and model neuron 1 (g_model1>bio1 = 100 nS, g_bio1>model1 = 50 nS). The first order resetting alone is shown for both the biological (long dashes) and model (short dashes) neuron (Fig. 7 A1). The addition of second order resetting (Fig. 7 A2) produces only a small change in the PRC for the biological neuron, but shifts the PRC for the model neuron upward, corresponding to a constant delay of ∼200 ms, the opposite of what was observed in Fig. 5, A1 and A2. When second order resetting is ignored (Fig. 7 B1), there is no intersection point between the two recovery versus stimulus interval curves of the neurons, and hence no 1:1 mode is predicted to exist. The contribution of F₂ to the stimulus interval for the model neuron 1 causes the dotted curve in Fig. 7 B2 to be shifted upward, as explained above, such that an intersection with the curve for the biological neuron 1 does indeed occur. Fig. 7 C1 shows the average intrinsic period of the biological neuron to be 1.5 s and that of the model 1 neuron to be 1.09 s. Fig. 7 D shows the 1:1 phase-locked mode obtained when the circuit was closed. Using the intersection point in Fig. 7 B2 resulted in a predicted value of 0.773 s for t_rmodel1 compared to an actual average value of 0.776 s, corresponding to a 0.4% error, a predicted value of 0.656 s for t_rbio1 compared to the actual average value of 0.746 s, corresponding to a 12% error, and a prediction for the phase-locked period of 1.429 s compared to the actual average period of 1.522 s (6% prediction error). A close inspection of the actual membrane potential recordings in the phase-locked mode (Fig. 7 D) shows that although the duration of the burst in the biological neuron remained near its intrinsic value of 0.55 s, the model neuron burst duration was 0.5 s during phase-locked mode compared to an intrinsic value of 0.28 s. The significant increase in the burst duration of the model neuron 1 for this particular hybrid network (g_model1>bio1 = 100 nS, g_bio1>model1 = 50 nS) during the phase-locked mode causes the biological neuron to be released from inhibition at a later time than would be predicted by the intrinsic burst duration, potentially resulting in a longer recovery interval for the biological neuron 1 in closed loop than predicted. This in turn could lead to the observed systematic underestimation of the recovery of the biological neurons, which is evident in the summary in Fig. 8 B. The lack of an intersection that satisfies all periodicity constraints in the absence of consideration of the effects of F₂ as shown in Fig. 7 B1 is not unique to this example but rather occurred in 21 out of the 164 tested circuits, including all those listed in Tables 1 and 3 with g_bio1>model1 ≥ 10 nS at two values of g_model1>bio1 (100 nS and 500 nS) and model neuron 5 coupled to biological neuron 1 (g_model5>bio1 = 100 nS).

Fig. 8 summarizes the predictions made using the hybrid circuit with g_model1>bio4 = 100 nS at eight values of g_bio1>model1. A shows the recovery interval for the model neuron, B shows the recovery interval for the biological neuron, and C shows their sum, which is the period. The experimentally measured values are given in all panels by black squares connected by black lines. There is no experimental value at g_bio1>model1 = 2 nS, because according to the entry in Table 1, a complex mode was observed although 1:1 phase locking was predicted at this value. The predictions using the average resetting are given by red circles connected by red lines, and in this case do not overlay quite as well with the experimental data as the example in Fig. 6, but the discrepancy is still within the variability of the data. The predictions based on the upper and lower (triangles on the blue and green curves) envelopes still largely bracket the experimental results. On average, the recovery interval of the model neuron 1 was overestimated by 8%, whereas the recovery interval of the biological neuron was underestimated by ∼10% (Table 3). The prediction error for the phase-locked period was only ∼5% (see next section). Note that the observed phase-locked period (Fig. 8 C) is again always close to the intrinsic period (1.50 s) of the biological neuron at all values of g_bio1>model1.

Quantitative summary of results

Table 3 gives a quantitative summary of the accuracy of the predictions. For each of the 20 combinations of a given model neuron and a given biological neuron at a fixed value of g_model>bio, we summed the error over all values of g_bio>model and reported the statistics. The values for the examples given in Figs. 6 and 8 are shown in bold. The first example (Fig. 6) shows that over a range of seven g_bio4>model1, the average error was 1.3 ± 0.8 and −0.3 ± 1.1, respectively, for the biological and model neuron recovery intervals. These numbers are at the lower end of those we observed and reflect the excellent fit shown in Fig. 6, A and B. The error when F₂ was not considered was −0.7 ± 1.3 and −2.2 ±1.9, which indicates that the contribution of F₂ did not improve (or worsen) our prediction in this case. The other example, given in Fig. 8 over a range of eight values of g_bio1>model1 for which phase-locked modes were observed, gives average errors that were at the high end of those observed with the correction for F₂ included, −10.2 ± 3.1 and 8.0 ± 3.3, respectively, for the biological and model neuron recovery intervals. In addition, this example illustrates an extreme case in which the existence of 1:1 phase-locked modes could not be predicted without the contribution of F₂. The consideration of second order effects was required for the prediction of the observed 1:1 phase locking in a total of three of the combinations (specifically in 21 individual circuits as described above). In several cases, such as for model neuron 7 combined with biological neuron 2 or 4, considering the contribution of F₂ drastically reduced the average error in the recovery interval of the biological neuron, whereas in others, such as model neuron 2 or 6 coupled with biological neuron 3, the average error in the recovery interval of the model neuron was drastically reduced. In all but two combinations, the mean prediction error of both recovery intervals was smaller when the contribution of the second order PRC to the total phase resetting was considered. In all but five combinations, the SE of both recovery intervals was smaller when the contribution of the second order PRC was considered.

Although in general this data set did not lend itself to a test of the predictive efficacy of including F₂ in the stability criterion, the case of the hybrid network g_bio1>model5 = 1 nS and g_model5>bio1 = 100 nS did provide such a test. The steady phase of the biological neuron in the predicted 1:1 mode was ϕ_bio1 = 0.754, and the slopes of the PRCs at that phase were m_1bio1 = 0.4, m_2bio1 = 0.0, m_1model5 = −0.733, and m_2model5 = 0.263. Based on the first order stability criterion, the absolute value of the characteristic root λ is 1.33, which is greater than 1, so that the steady mode is predicted to be unstable. However, considering the contribution of the second order PRC decreases the maximal value of λ to 0.78, which correctly predicts that the observed 1:1 phase-locked mode will be stable. In a total of 9 out of 164 circuits, stable phase locking would be incorrectly predicted to be unstable if the effects of the second order phase resetting on the stability criterion were not taken into account. This provides further support for the importance of considering the contribution of second order resetting.

If there is an offset between the activation of the coupling with the postsynaptic neuron and the initiation of a burst in the presynaptic neuron, the Methods explain how to implement the periodicity constraints, and an example is included in the Supplementary Material. Such an offset could in principle result from axonal delays or from graded synaptic coupling in which activation of the coupling precedes burst onset. All of the entries in Table 3 corresponding to circuits composed of biological neuron 3 and model neuron 2 were calculated in the presence of such an offset, and the errors were neither excessively high nor low compared to the remainder of the circuits that had no such offset.

In Table 3, the signs of the errors for the recovery interval in the two neurons are frequently opposite. The peculiarities of the specific PRCs utilized, including the flat region in the PRC that is often observed in a range of phase bracketing 0.5 in the biological neuron (Fig. 5 A1, long dashes) and the linearity with unit slope observed in the PRC for the model neuron, frequently combined to produce compensatory errors in the two recovery intervals as follows. In the flat region of the first order PRC for the biological neuron, F_1bio(ϕ) = φ₀ = constant, thus t_rbio = P_0bio (1 − ϕ + φ₀), which produces a region of negative unit slope in the curve indicated by dashes in Figs. 5 and 7, B1 and B2. Since the slope of the biological neuron recovery interval curve near the intersection point is −1, and the model recovery interval resembles a straight vertical line, a shift in the model recovery interval produces an equal and opposite shift in the recovery interval for the biological neuron and leaves the prediction of the period essentially constant in the region of negative unit slope in the recovery interval curve produced by the flat region in the first order PRC for the biological neuron. Additional examples of this phenomenon are given in the Supplementary Material.

Validity of the mathematical analysis

Mathematical biology requires that assumptions be made to map the mathematics onto aspects of the biological system under study. These assumptions always involve an approximation, and the validity of the mathematical model in a given situation depends upon the quality of the approximation. Most previous theoretical work on the phase-locked modes in biological oscillators has focused on weakly coupled, simple integrate and fire (IF) neurons (Hansel et al., 1995; Mirollo and Strogatz, 1990; van Vreeswijk et al., 1994), IF neurons with arbitrarily strong pulse coupling (Bressloff and Coombes, 1998), spike response methods applied to neurons with standard dynamics, i.e., to IF and type I neurons (Gerstner et al., 1996), or spiking neurons with negligible second order resetting (Goel and Ermentrout, 2002). These studies are all based on assumptions that are inappropriate for the bursting neurons, strong coupling, and relaxation oscillator dynamics observed in the hybrid circuits in this study.

The mathematical model that we present is a discrete, cycle-by-cycle mapping of the activity of two coupled neurons in the hybrid circuit that requires no knowledge of the equations governing the dynamics of each neuron, only the empirically determined PRC. The following essential assumptions were made in this study: 1), the complex composed of the soma and extended neuritic trees of both PD neurons and the AB neuron to which they are electrically coupled can be represented by a single, lumped oscillator termed the biological neuron, 2), both the biological neuron and the model neuron in the circuit are noiseless limit cycle oscillators, 3), each neuron either returns close to its original unperturbed limit cycle before the next input is received or if the cycle period changes over time, it returns to a limit cycle identical to the original one remapped to the new cycle period, and 4), the input received by each neuron has the same effect in the closed loop circuit as in the open loop circuit used to generate the phase resetting. The accuracy of the results was satisfactory for 161 of the 164 circuits (see Table 1), was quite good in many cases (see Fig. 6), and in almost every case was within the 10% variability in the period of the biological neurons. To fully understand the approximations that are introduced by the assumption, and to demonstrate where vigilance is required by others attempting to implement similar methods, we elaborate below on violations of the assumptions and their potential consequences.

The first assumption regarding the lumped AB/PD oscillator kernel was violated in several instances but fortuitously did not affect the accuracy of our predictions. In some cases, a pulse that effectively hyperpolarized the PD neuron from which we were recording did not terminate an ongoing burst, complete with spikes, presumably occurring elsewhere in the AB/PD complex due to insufficient space clamp (Fig. 9 A). In another case, a perturbation applied around a phase of 0.75 had variable results, producing an apparent discontinuity in the PRC at that phase. Some trials resulted in a missed cycle in which an attenuated subthreshold depolarization without spikes was initiated during the perturbation instead of a burst (Fig. 9 B), and on other trials no depolarization occurred until after the inhibition (Fig. 9 C). The missed cycle produced an unusually long delay compared to cases in which a rebound burst was observed just after PD was released from inhibition and hence the discontinuity in the PRC (Fig. 9 D). We theorize that no errors were introduced into our analysis by these violations of the assumption because the observed phase locking usually occurs near the middle of the AB/PD cycle, when the input conductance is low and hyperpolarization may propagate more effectively to the oscillator kernel, presumably located in the neurites of the AB neuron.

The second assumption was more problematic, and the approximation introduced herein is likely responsible for the bulk of the numeric errors. One source of error is hardware dependent: the dynamic clamp can introduce numerical error since the setup used has an average integration time step of ∼0.06 ms, but individual time steps have variable lengths around that value. Some variability was observed in the duty cycle of the model neurons in some simulations as a result of numerical errors; yet interestingly, much less variability was observed in the period. The biological neuron was a more significant source of error, and a conclusion of this study is that the method can, at least in some cases, tolerate levels of noise present in real biological systems. The major sources of noise in the biological portion of the experimental setup were the trial-to-trial variability inherent in the experimental determination of the PRC due to noise and fluctuations in the modulatory state of the biological circuit and the cycle-to-cycle variability of the period. Since a constant intrinsic period is assumed to calculate both the phases at which perturbations are applied and the resultant normalized delay or advance that results from the perturbation, the cycle-to-cycle variability in the period is likely to be responsible for much of the trial-to-trial variability in the PRC. In the experimental setup described in this article, the pyloric pacemaker still receives many modulatory inputs from anterior ganglia in the stomatogastric nervous system. In addition, gastric mill events can cause a transient change in the period of the pyloric pacemaker (Bartos and Nusbaum, 1997; Bartos et al., 1999; Marder et al., 1998; Mulloney, 1997; Nadim et al., 1998, 1999; Thuma and Hooper, 2002). In many cases, due to long-term trends, the period of the biological neuron immediately before the coupling in the hybrid circuit was turned on was different than when the PRC was measured. To utilize the PRC determined using a different intrinsic period to make predictions for that hybrid circuit, we were forced to make the additional assumption that the limit cycle with the new period is identical to the original one remapped to the new cycle period. This assumption guarantees a constant duty cycle, but there is insufficient evidence to support any particular mapping of the phase as cycle period is varied. The impact of the trial-to-trial variability of the PRC was addressed by using the average PRC for the predictions but also by comparing the spread of the predictions resulting from using the upper and lower envelopes of the experimental PRCs. We concluded that the trial-to-trial variability could account for much of the observed errors. The cycle-to-cycle variability of the period was 10%, and the predictions of the average period fell within that range.

The third assumption regarding the return to the limit cycle before receipt of the next input is required for the PRCs obtained at specific phases along the limit cycle in the open loop configuration to apply in the closed loop. Potential violations of this assumption could be introduced in at least two ways: higher order PRCs in response to a single pulse and temporal summation of slow processes in response to a train of pulses. In the first scenario, the deviation from the limit cycle produced by a single pulse does not die out within one trip around the limit cycle. F₁ measures the portion of the resetting that occurs during one cycle before the next burst, and F₂ presumably measures the resetting that occurs after burst initiation but still before the returning to the original point on the limit cycle at which the perturbation was received. If higher order PRCs are not negligible, then errors are introduced by the application of the open loop PRC because the phase cannot be precisely determined unless the trajectory is very near the original limit cycle. The manner in which F₂ is tabulated does not guarantee that all of the resetting occurs before a return to the original phase. Furthermore, in some hybrid networks, such as biological neuron 3 coupled with model neuron 2 (data not shown), the third order PRCs are not negligible, and could be responsible for the observed systematic errors. In this hybrid circuit the duration of the perturbation is very long compared to the cycle period in the model, so the assumption of pulsatile coupling could be violated because the trajectory cannot return to the limit cycle during the interburst. The burst duration in model neuron 3 is 1.06 s, which is 84% of the 1.26-s intrinsic period of biological neuron 2 before coupling with model neuron 3 (see Supplementary Material for an example). Although the mean prediction error of the recovery interval for this circuit was small (0.6%), the SE (9.2%, Table 3) is a better indicator of prediction accuracy. A second, distinct source of error could result if a single perturbation causes a small deviation from the limit cycle, but repeated, periodic applications of the perturbation cause the trajectory to move farther away from the limit cycle such that the distance summates temporally. In this manner, the presence of slow conductances could result in systematic errors.

The final assumption is that the input received by each neuron has the same effect in the closed loop circuit as in the open loop circuit used to generate the phase resetting. This does not require the burst to be identical in both cases, because in some instances the phase resetting can be insensitive to duration at the phase locked point (Oprisan et al., 2003), but a change in burst duration can also be an important source of error. For some model neurons, the duration of the first burst (see Fig. 1 D) after the perturbation is applied can differ from its intrinsic value measured in the open loop setup, and Figs. 5 D and 7 D show that the duration of the burst in the model neuron differs in the open and closed loop conditions. Thus, violation of this final assumption is a potentially important source of error, although in this article the errors introduced were not excessively large in magnitude (see Supplementary Material).

Applications to CPGs and other phase-locked circuits

An interesting feature of certain CPGs is that the characteristic pattern is preserved over a large range of frequencies (phase constancy; see Hooper, 1997). The promise of our methodology is that it may provide insight both into how the pattern is maintained and the circumstances that cause it to break down, such as a change in the slope of the first or second order PRC that renders the mode unstable, or a change in the magnitude of the PRC that makes it impossible for the periodicity constraints to be satisfied. Fig. 7 provides an excellent illustration of how a change in the second order resetting, for example, can cause a phase-locked mode to appear or disappear. In addition, to our knowledge, no other analysis of the existence or stability of phase locking in neural circuits has considered the effects on the second cycle after a perturbation is received. As we have shown here, such effects are particularly important for bursting neurons. However, pyramidal neocortical neurons have been shown to exhibit second order resetting as a result of a depolarizing pulse (Reyes and Fetz, 1993); therefore, theories of phase locking in cortical networks may also need to include the effects of F₂. The robustness of our methods to the presence of delays is also noteworthy.