Research ArticleNEUROSCIENCE

Synaptic mechanisms for motor variability in a feedforward network

See allHide authors and affiliations

Science Advances  19 Jun 2020:
Vol. 6, no. 25, eaba4856
DOI: 10.1126/sciadv.aba4856


Behavioral variability often arises from variable activity in the behavior-generating neural network. The synaptic mechanisms underlying this variability are poorly understood. We show that synaptic noise, in conjunction with weak feedforward excitation, generates variable motor output in the Aplysia feeding system. A command-like neuron (CBI-10) triggers rhythmic motor programs more variable than programs triggered by CBI-2. CBI-10 weakly excites a pivotal pattern-generating interneuron (B34) strongly activated by CBI-2. The activation properties of B34 substantially account for the degree of program variability. CBI-10– and CBI-2–induced EPSPs in B34 vary in amplitude across trials, suggesting that there is synaptic noise. Computational studies show that synaptic noise is required for program variability. Further, at network state transition points when synaptic conductance is low, maximum program variability is promoted by moderate noise levels. Thus, synaptic strength and noise act together in a nonlinear manner to determine the degree of variability within a feedforward network.


When one repeats a motor act, the behavior is often variable, even with extensive practice (16). Behavioral variability can be manifested in various forms, one of which is the cycle-to-cycle variability sometimes observed during rhythmic motor behaviors (4, 7, 8). It has been suggested that this and other forms of variability have adaptive benefits or computational advantages (3, 5, 912). Although previous work has determined how variability in spike trains can arise at the single-cell level in the cortex (2, 5, 9, 13, 14) and spinal cord (15) of vertebrates and in the nervous systems of invertebrates (16), it is virtually unknown how variability at the circuit level is generated. Moreover, previous studies have focused, to a large extent, on the variability that is observed when neurons receive concurrent excitation and inhibition. In contrast, our studies explore how variability can be generated in a feedforward motor network where circuit elements are driven primarily by excitatory inputs (1722).

We study the highly tractable Aplysia feeding network, which has a small number of well-defined cell types/synaptic connections (22). Aplysia feeding behavior is repetitive, with each cycle consisting of a radula protraction-retraction sequence. Cycle-to-cycle variability has been demonstrated both in vivo and in vitro (4), but the mechanisms responsible for this variability have not been identified. The feeding network is a feedforward circuit consisting of at least three “levels” (2224): command-like neurons [cerebral-buccal interneurons (CBIs)], pattern-generating interneurons (e.g., B34), and motoneurons (e.g., B8) (Fig. 1A). This type of circuit is common and is present in both invertebrates and vertebrates (1722, 2527) . In the Aplysia feeding and other networks, the initiation of motor programs depends on feedforward excitation of all three types of elements. Although inhibition is also present in the network, it serves to terminate, rather than activate, neurons. Command-like neurons are cerebral ganglion neurons that are driven by appropriate excitatory sensory inputs under physiological conditions. In this study, we induced analogs of feeding behavior by direct stimulation of command-like neurons in the cerebral ganglion. This allowed us to focus on mechanisms underlying variability in the motor circuit itself rather than taking afferent input as a source of noise, as has been done in other studies (5, 11, 28).

Fig. 1 Feedforward circuit, CBI-10 activity during feeding, and motor programs elicited by CBI-10 and CBI-2.

(A) Schematic diagram of the feedforward circuit in the Aplysia feeding system: from command-like neurons, to pattern-generating neurons [e.g., protraction interneurons (PIs), B34 and B63 and retraction interneuron (RI), B64], and protraction (PM) and retraction (RM) motoneurons. Each cycle of the protraction-retraction sequence is mediated primarily by the alternating activity of B63 and B64. Connection symbols: open triangle, excitation; filled circle, inhibition; s, slow connections. (B) Morphology of a left CBI-10 revealed by carboxyfluorescein injection. Its main axon projects anteromedially, travels ventrally and turns back toward the soma, and then projects posteromedially before making a U-turn to the cerebral-buccal connective (CBC) (see fig. S1). Note the extensive processes in the medial region near the soma. The front image originated from three fluorescent images at different depths stacked together. The background image was a bright-field image stacked with the fluorescent image. AT, anterior tentacular; LLAB, lower labial nerve; ULAB, upper labial nerve; CPe, cerebral-pedal connective; CPl, cerebral-pleural connective. (C and D) CBI-10 in a semi-intact preparation. Feeding behavior was monitored by buccal mass “pressure.” (C) Upon delivery of dried seaweed near the mouth, CBI-10 was activated during feeding episodes. Completion of the feeding sequence was verified by exiting of seaweed from the cut end of the esophagus. (D) Stimulation of a single CBI-10 through DC current (bar) elicited four feeding-like responses. (E and F) Multiple cycles of motor programs elicited by CBI-10 (E) or CBI-2 (F) from a single in vitro preparation. Open bar, protraction; filled bar, retraction. (G to I) Paired t test comparing the coefficient of variation (CV) of protraction duration (G) (t7 = 3.861, **P = 0.0062), duty cycle (H) (t7 = 5.147, **P = 0.0013), or interspike interval (ISI) of B34 (I) (t4 = 5.964, **P = 0.004) when CBI-10 or CBI-2 was stimulated at 10 Hz.

We demonstrate that two command-like neurons (CBI-2 and CBI-10) in Aplysia initiate the same type of behavior/program, but there is a prominent difference in the variability of evoked motor programs. This provides an intrinsic advantage. Variability can be studied by making a within-preparation comparison (i.e., within a single animal, programs induced by CBI-2 can be compared to programs induced by CBI-10). We focused on these two command-like neurons because other neurons either have been extensively studied (19, 23, 24, 29) or are involved in other types of behavior/program (30). Our data strongly suggest that the difference in variability results from differential activation of a pivotal pattern-generating interneuron (B34). B34 has a relatively high spiking threshold, so summation of excitatory postsynaptic potentials (EPSPs) must occur. CBI-2 and CBI-10 both induce EPSPs in B34 that vary in amplitude (i.e., synaptic noise is observed at both the CBI-2 to B34 and CBI-10 to B34 synapses). Although this synaptic noise in the biological network could be due to presynaptic factors such as fluctuations in neurotransmitter release and/or randomness in postsynaptic responsiveness, we successfully modeled it as a type of presynaptic noise. The CBI-2 to B34 and CBI-10 to B34 synapses differ in that CBI-10–induced EPSPs are significantly smaller. In this range, variability in EPSP amplitude effectively translates into variability in synaptic activation (i.e., the B34 firing frequency is more variable when activity is triggered by CBI-10). Our computational analysis demonstrates the essential importance of both synaptic noise and low synaptic strengths in generating variable motor programs. Further, we show that synaptic noise can expand the working range of synaptic strength by which a command-like neuron uses to drive programs. Because feedforward circuits that include command-like neurons and pattern-generating neurons are common (1722, 2527, 31), our findings are likely to be broadly applicable.


Characterization of the command-like neuron CBI-10

In Aplysia, there is more than one feeding command-like neuron (fig. S1) (19, 23, 24, 29). CBI-2, identified previously, reliably drives feeding-like behavior and motor programs (18, 19, 23, 24, 29). Here, we describe a previously uncharacterized CBI, CBI-10. Similar to other CBIs including CBI-2, CBI-10 projects its axon to the buccal ganglion through the ipsilateral cerebral-buccal connective (CBC) (Fig. 1B and fig. S1). To record activity of CBI-10 during feeding behavior, we conducted experiments in a semi-intact preparation where normal feeding episodes can be elicited with seaweed delivered to the mouth (19). A single CBI-10 was activated during seaweed-induced feeding and fired at an average frequency of 6.77 ± 0.34 Hz (indicated by open bars in Fig. 1C) (n = 18 cycles from four preparations). This was somewhat lower than the frequency reported for CBI-2 (10.6 Hz) (19). In the same semi-intact preparations, feeding-like movements could also be triggered in the absence of seaweed by directly activating a single CBI-10 (Fig. 1D).

CBI-2 and CBI-10 elicited programs with different degree of variability

In the isolated central nervous system (CNS), feeding motor programs can be triggered by stimulating CBI-2 at 10 Hz (Fig. 1F) and are monitored by recording activity in the B8 motoneuron (as a monitor of radula closing/retraction) and by recording from the I2 nerve (as a monitor of radula protraction) (18, 19, 29, 32). For comparative purposes, we stimulated CBI-10 at the same frequency, i.e., 10 Hz, and found that continuous stimulation of CBI-10 produced multiple cycles of programs that were similar in that each cycle had two phases: Radula protraction was followed by retraction (Fig. 1E). Programs differed, however, in that activity induced by CBI-10 appeared to have a more variable protraction duration and duty cycle. To quantify this, we computed the coefficient of variation (CV) of protraction duration and the duty cycle (Fig. 1, E and F). Both CVs were larger when programs were triggered by CBI-10 (Fig. 1, G and H). Further, we computed the CV of the interspike interval (ISI) of B34. B34 is a protraction-interneuron that is activated during motor programs triggered by both CBIs (Fig. 1, A, E, and F). The CV of the ISI is a representation of the variability of neuronal firing (2, 5, 9, 13, 14). The CV of the B34 ISI was higher when activity was triggered by CBI-10 (Fig. 1I).

To determine whether programs might be more variable if we stimulated CBI-2 at a lower frequency, we triggered programs by stimulating CBI-2 at 5 to 7 Hz. This is the minimal firing frequency (i.e., threshold) required for eliciting a motor program. There was no change in the CV of protraction duration, and there was a decrease (rather than an increase) in the CV of the duty cycle (fig. S2). These data suggest that the regularity of CBI-2–elicited programs was not dependent on the CBI-2 firing frequency. Thus, in a single neural network, the same motor program elicited by different command-like neurons can have different variability. Motor variability has been challenging to study because variability is present between individuals/preparations (4, 7, 8), but the different degree of variability triggered by the two CBIs in the same preparations allows us to overcome this difficulty.

Circuit mechanisms underlying motor variability: Roles of pattern-generating interneurons (B34)

To examine why protraction duration and duty cycle are less variable when motor programs are triggered by CBI-2, we sought to determine whether CBI-2 more strongly and less variably activates protraction interneurons. We focused on B34 since previous studies have shown that changes in its activity affect protraction duration (19, 32), and it did show less variable firing in CBI-2–elicited programs (Fig. 1I). We quantified B34 firing frequency during CBI-2– and CBI-10–evoked programs. B34 fired at a higher frequency when programs were triggered by CBI-2 (Fig. 2A). To determine whether there was a correlation between B34 activity and variability, we plotted the B34 firing frequency versus the CV of protraction duration (Fig. 2B) and the CV of duty cycle (Fig. 2C) in motor programs elicited by CBI-2 or CBI-10. In both cases, data were inversely correlated.

Fig. 2 Correlational and causal roles of pattern-generating interneurons (B34) in the generation of motor programs with different degree of variability.

(A) Activity of B34 during CBI-10 (n = 19) versus CBI-2 (n = 22) evoked programs (unpaired t test, t39 = 5.136, ****P < 0.0001). (B and C) B34 activity was inversely correlated with the CV of protraction duration (B) (r = −0.6615, r2 = 0.4376, P < 0.0001, n = 58) and duty cycle (C) (r = −0.6245, r2 = 0.39, P < 0.0001, n = 58). Lines, linear regression lines; r, correlation coefficient. (D to G) Subthreshold depolarization of B34s (E) (bars under B34 recordings, n = 5) in the programs elicited by CBI-10 (8 Hz). Paired t test comparing the CVs under control conditions with those obtained with depolarization (dep) [(F) protraction duration: t4 = 3.991, *P = 0.0163; (G) duty cycle: t4 = 3.078, *P = 0.037]. (H to K) Hyperpolarization of B34s (I) (bars under B34 recordings, n = 8) in programs elicited by CBI-2 (5 Hz). Paired t test comparing the CVs under control conditions with those obtained with hyperpolarization (hyp) [(J) protraction duration: t7 = 5.198, **P = 0.0013; (K) duty cycle: t7 = 6.797, ***P = 0.0003]; error bars, SEM. c-B34, contralateral B34. Both ipsilateral and contralateral B34s were depolarized (E) or hyperpolarized (I), and the two B34s are coupled electrically (44). Thus, subthreshold depolarization of B34s (E) made the programs elicited by CBI-10 less variable, whereas hyperpolarization of B34s (I) made programs elicited by CBI-2 more variable.

Last, we determined whether manipulating B34 activity would alter program variability. During CBI-10–induced programs, depolarizations that increased B34 firing frequency reduced the CVs of both protraction duration and duty cycle, recapitulating the properties of CBI-2–elicited programs (Fig. 2, D to G). During programs evoked by stimulating CBI-2 at threshold frequency, hyperpolarizations that prevented B34 from firing had the opposite effects, i.e., variability increased (Fig. 2, H to K). These data indicate that B34 plays an essential role in determining whether motor programs are variable.

Difference in synaptic strength potentially important for different degree of variability

What are the synaptic mechanisms and intrinsic properties that could account for the weaker B34 activity during CBI-10–elicited programs? A previous study (24) demonstrated that CBI-2 makes a monosynaptic excitatory connection with B34 and that the EPSPs are chemical, cholinergic, and facilitated. We found that CBI-10 also made a monosynaptic excitatory connection with B34 and that the EPSPs were chemical, cholinergic, and facilitated as well (Fig. 3, A to C, and fig. S3, A and B).

Fig. 3 Strength and noise in synapses to B34 from CBI-10 versus CBI-2.

(A to C and E to G) Synaptic strengths and facilitation in physiology (A to C) and in the model (E to G). Both CBI-10 (A and E) and CBI-2 (B and F) stimulations (10 Hz for 2 s) elicited monosynaptic EPSPs in B34, but the EPSP amplitude was larger for CBI-2 than for CBI-10. (C and G) Group data from physiological experiments (n = 8) (C) and the model (n = 6) (G). Paired t test comparing average amplitudes of the EPSPs from CBI-10 versus CBI-2. In the physiology (C): first EPSPs, t7 = 2.512, *P = 0.040; last EPSPs, t7 = 6.462, ***P = 0.0003; in the model (G): first EPSPs, t5 = 11.55, ****P < 0.0001; last EPSPs, t5 = 10.39, ***P = 0.0001. Paired t test comparing average amplitudes of the first and last EPSPs. In the physiology (C): CBI-10, t7 = 3.494, *P = 0.010; CBI-2, t7 = 7.318, ***P = 0.0002; in the model (G): CBI-10, t5 = 2.914, *P = 0.0332; CBI-2, t5 = 7.579, ***P = 0.0006; error bars, SEM. (A, B, D to F, and H) Synaptic noise in physiological experiments (A, B, and D) versus the model (E, F, and H). CBI-10 or CBI-2 was stimulated at 10 Hz for 2 s for 10 times (only five alternate examples from the 10 trials are shown). In the physiology (D), the average CVs for the last EPSPs from CBI-10 to B34 were 0.221 (n = 6) and from CBI-2 to B34 were 0.194 (n = 6); paired t test, t5 = 0.739, P = 0.493; n.s., not significant. Error bars, SEM. In the model (H), synaptic noise is modeled as a type of presynaptic noise, and the average CV for the last EPSPs from CBI-10 to B34 was 0.206 with SD of noise σ = 0.18 (n = 6) and from CBI-2 to B34 was 0.184 with σ = 0.15 (n = 6) (paired t test, t5 = 1.135, P = 0.308), matching physiological data. For group data in (G) and (H), six simulations of 10 time series of EPSPs for 2 s were performed. The first simulation result was obtained with the default values of maximum conductance and plasticity parameter as defined in table S3. Five more sets of these two parameters were generated by adjusting them up or down randomly (up to 5% of their default values). The schematic diagram illustrates the circuit elements and their connections in the simulation. Connection symbols: open triangle, excitation.

CBI-10–induced EPSPs were, however, smaller (Fig. 3, A to C). This was apparent when the two CBIs were stimulated at the same frequency for the same duration in the same preparation (Fig. 3, A to C). We also determined intrinsic properties of B34. B34 had a relatively high spike threshold and low input resistance compared with another protraction interneuron B63 (fig. S3, C and D), suggesting that it is relatively difficult to activate. Thus, to activate B34 at high frequency, a strong excitatory synapse from a CBI, e.g., CBI-2, is required. Overall, these data partially explain why B34 fires at a higher frequency when motor activity is triggered by CBI-2.

Synaptic noise in the biological network and in model networks

The above data demonstrate that there is a correlation between program variability and B34 activity (i.e., when programs are triggered by CBI-10, there is less B34 activity and programs are more variable, and the reverse is true for programs triggered by CBI-2). A further question is why this correlation was observed. Previous studies demonstrated that B34 activity affects protraction duration (19, 32). However, these studies did not analyze variability.

Our data show that there is variability in the B34 activity, i.e., the B34 ISI is variable (Fig. 1I). This suggests that there is some type of synaptic noise in the network. To characterize it, we scrutinized synaptic potentials induced in B34. We found that EPSPs induced by CBI-10 and CBI-2 did vary in amplitude (Fig. 3, A to C). That is, when EPSPs facilitate, a progressive increase in amplitude is expected. This was not always the case, e.g., an occasional decrease or unexpectedly large increase in EPSP amplitude was observed. To formally characterize this synaptic “noise,” we stimulated both CBIs at 10 Hz for 2 s 10 times. For the 10 trials, we measured the amplitude of the last EPSP and computed the CVs for CBI-10 (0.221 ± 0.018) and CBI-2 (0.194 ± 0.025) datasets (Fig. 3D). There was no statistically significant difference. Thus, EPSPs from both CBIs are variable.

To demonstrate a causal role for this synaptic noise, it is necessary to manipulate it, e.g., remove it. This is difficult to accomplish in a biological network. We therefore constructed a simplified computational model of the Aplysia feeding network (fig. S4A). This model includes compartmental models of CBI-2 and CBI-10, B34, and two identified plateau-generating interneurons that are essential for generating protraction and retraction (see Materials and Methods and tables S1 to S3). To build the model, we also identified the synaptic connection from CBI-10 to B63 as excitatory (fig. S4B). The modeled CBIs-B34 synapses were adjusted so that synaptic conductances and facilitation mimicked physiological data (Fig. 3, E to G). Synaptic noise was modeled as presynaptic by multiplying the transmitter variable with a Gaussian random variable (see Materials and Methods). To compare the resulting synaptic noise, we elicited EPSPs in B34 by stimulating both CBIs at 10 Hz for 2 s 10 times. Similarly, we computed the CVs for the last EPSPs in each dataset. As expected, these CVs (Fig. 3H) were similar to the CVs calculated from physiological data. These results indicate that it is possible to model variability in synaptic input from the CBIs to B34 using presynaptic noise. However, we cannot exclude the possibility that there is a postsynaptic contribution in the biological network.

We also considered the possibility that synaptic noise was generated via a different type of mechanism, i.e., by some type of “background input” to B34. To explore this possibility, we modeled background input as a synaptic current(s) originating from a spike train of a neuron or neurons outside the network, and the simulated spike train obeyed a Poisson distribution (see Materials and Methods and fig. S5). Specifically, we modeled three types of background input in the absence (fig. S5, A to F) or presence (fig. S5, G to L) of presynaptic noise [as in Fig. 3 (E to H)]: (i) as concurrent excitatory and inhibitory synaptic currents, each originating from a Poisson spike train with a rate of 2 Hz (fig. S5, A, B, G, and H); (ii) as an excitatory synaptic current originating from a Poisson spike train with a rate of 10 Hz (fig. S5, C, D, I, and J); or (iii) as an inhibitory synaptic current originating from a Poisson spike train with a rate of 10 Hz (fig. S5, E, F, K, and L). Each type of simulation was performed at two conductance levels (higher in the left two panels and lower in the right two panels). We also determined the CV of the last EPSPs in B34 when CBI-10 or CBI-2 was stimulated at 10 Hz for 2 s 10 times (bar graphs in fig. S5). In all cases, background noise modeled in this fashion did make the EPSPs variable, as shown by the CVs. The CV tended to be higher for the CBI-10–B34 synapse than the CBI-2–B34 synapse, but the CVs were similar in the biological network (Fig. 3D). Moreover, background input also tended to cause large baseline changes in B34 membrane potential, which were especially apparent when CBI-10 was stimulated. These baseline shifts in membrane potential were not apparent in physiological recordings. Thus, these data suggest that this type of background noise is likely not a main factor in the Aplysia network. Consequently, presynaptic rather than background noise is used in all subsequent simulations described below.

Synaptic mechanisms underlying motor variability: Roles of synaptic noise and strength

When synaptic noise as modeled in Fig. 3 (E and F) was included in the model network (fig. S4A and table S1), both CBIs could trigger programs. As in the biological network, CBI-10–elicited motor programs were more variable than those elicited by CBI-2 [Fig. 4, A and C to F; see also Fig. 5 (A and D)]. Further, the variability of programs could be manipulated by changing overall synaptic drive to B34, i.e., variability could be converted to CBI-2 program variability by depolarizing B34 during CBI-10–evoked programs (Fig. 4, A, B, G, and H) or vice versa, i.e., converted to CBI-10 program variability by hyperpolarizing B34 during CBI-2–induced programs (Fig. 4, I to L).

Fig. 4 The feeding network model behaves like the biological network.

(A and C) Comparison of programs elicited by CBI-10 (A) and CBI-2 (C) with synaptic noise (CBI-10 and CBI-2 were stimulated at 10 Hz). (B) Depolarization of B34 (bar) during CBI-10 programs with noise made the programs less variable. (D to F) Comparison of CVs of protraction duration (D) (paired t test, t9 = 7.963, ****P < 0.0001), duty cycle (E) (paired t test, t9 = 6.073, ***P = 0.0002), and B34 ISI (F) (paired t test, t9 = 14.97, ****P < 0.0001) from programs elicited by CBI-10 versus CBI-2. (G and H) Comparison of the CVs of protraction duration (G) (paired t test, t9 = 4.855, ***P = 0.0009) and duty cycle (H) (paired t test, t9 = 5.157, ***P = 0.0006) from programs elicited by CBI-10 when B34 was not depolarized (Normal) versus when B34 was depolarized (B34 dep). (I to L) Effects of B34 hyperpolarization on CBI-2–elicited programs. The programs were elicited by stimulation of CBI-2 at threshold frequency (6.5 Hz) (I and J). Hyperpolarization of B34 (B34 hyp) (bar) made the programs more variable (J). The current used to hyperpolarize B34 was set at 80% of the current that was just enough to completely prevent firing of B34. This was because 100% of the current prevents normal program generation in the model. B63 activity represents protraction, whereas B64 activity represents retraction. (K and L) Group data. Protraction duration: paired t test, t9 = 5.579, ***P = 0.0003; duty cycle: paired t test, t9 = 7.462, ****P < 0.0001. Error bars, SEM. All CVs were derived from averages of 10 simulations with a duration of 200 s each (there were five or more cycles in each simulation).

Fig. 5 Program variability in the computational model: Necessity of synaptic noise and roles of synaptic strength.

(A to H) Representative examples of motor programs elicited by CBI-10 or CBI-2 at different conductances [parameter analysis data are shown in (I) to (P)]. (A to D) With noise: Normal conductance from CBI-10 to B34, 0.68 and to B63, 1.59 (A) or high conductance from CBI-10 to B34, 1.1 and to B63, 1.59 (B); low conductance from CBI-2 to B34, 2.4 and to B63, 1.69 (C) or normal conductance from CBI-2 to B34, 3.4 and to B63, 1.69 (D). (E to H) Without noise: The same parameters as in (A) to (D). B63 activity represents protraction, whereas B64 activity represents retraction. (I to P) Parameter analyses of program variability for synaptic conductance from CBI-10 (I to L) or CBI-2 (M to P) to B34. The parameter (synaptic conductance from CBIs to B34) range was divided into orange (I to L) or purple (M to P) and white areas, with the data plotted as blue lines (with noise) or red lines (without noise). Orange area: No programs elicited by CBI-10; purple area: CBI-2–elicited programs were generated only by B63. White ones: Normal programs generated by both B34 and B63. Note the shift of the white areas with noise to the left relative to the data without noise for CBI-10 programs. One hundred data points were used for each plot. Each data point was derived from one simulation with a duration of 200 s (each simulation contains five or more programs). Labels in plots depict conductance parameters for examples shown in (A) to (H).

An important goal of the modeling was to determine the impact of removing synaptic noise when programs were triggered by CBI-10. Specifically, synaptic noise as modeled in Fig. 3E was removed. When we performed this manipulation, motor programs were no longer generated (Fig. 5E). These data indicate that synaptic noise can be essential for motor program generation. We hypothesized that removing synaptic noise from the CBI-10 circuit would be possible if the synaptic connection with B34 were stronger. We therefore performed a parameter analysis in which we systematically varied the strength of the CBI-10–B34 synaptic conductance and measured program variability in the presence and absence of synaptic noise. Even in the presence of synaptic noise, output patterns transitioned between two different states: a state in which motor programs were not generated [orange in Fig. 5 (I to L)] and a state in which they were [white in Fig. 5 (I to L)]. With synaptic noise, programs could, however, be triggered at a lower conductance. When motor programs were triggered, variability was apparent with synaptic noise but essentially absent without synaptic noise (see also fig. S6). The addition of synaptic noise had its greatest impact close to the point at which state transitions occurred, i.e., when the conductance was near the lowest value necessary for program generation.

To determine whether the CBI-2 circuit would behave similarly, we performed a similar conductance analysis. With CBI-2 stimulation, it was possible to trigger motor programs without synaptic noise at all modeled conductances. Again, however, output patterns transitioned between two states. At relatively low conductances, motor programs were abnormal (i.e., not all pattern-generating neurons were consistently active; B63 was active, but B34 was inactive; see Fig. 5G) [purple in Fig. 5 (M to P)]. At high conductances, programs were normal [white in Fig. 5 (M to P)]. As with CBI-10 stimulation, variability was apparent with synaptic noise but essentially absent without synaptic noise. Similarly, the addition of synaptic noise had its greatest impact close to the point at which there was a state transition. Together, these analyses indicate that synaptic noise is essential for generating variable motor programs. Further, the addition of synaptic noise is most effective when synaptic connections are relatively weak. At lower than normal conductances (Fig. 5C), CBI-2 can elicit programs that are as variable as those evoked by CBI-10 at normal conductances (Fig. 5A).

Impact of different levels of synaptic noise on motor variability

Last, we sought to determine how varying the level of synaptic noise would affect program variability. To vary the noise level, we altered σ, SD of the Gaussian variable. When programs were triggered by CBI-10 and the strength of the CBI-10–B34 synapse was normal, variability was highest at moderate noise levels (Fig. 6, A and B). When noise was increased further, variability did not increase; instead, it decreased. Note that when we extended the range of noise level (fig. S7, A and B), program variability reached a plateau with increasing noise. Other modeling experiments determined whether similar results would be obtained in the CBI-2 circuit. When the strength of the CBI-2–B34 synapse was normal, variability slightly increased as did σ (Fig. 6, G and H, and fig. S7, C and D).

Fig. 6 Program variability in the computational model: Roles of different levels of synaptic noise.

(A to D) Effects of different levels of noise (SD, σ) on the CVs of CBI-10–elicited programs. CVs of protraction duration and duty cycles were the highest at moderate noise levels near transition points when the CBI-10 to B34 conductance was normal at 0.68 (A and B). CVs increased with noise levels when conductance was high at 1.1 (C and D), although the highest CVs were lower than the low CVs in (A and B). The red vertical dashed lines in (A and B) indicate the SDs used to generate data in Fig. 5 (I and J). Orange areas, no programs. (E to H) Motor programs elicited by CBI-2. The CBI-2 to B34 conductance was low at 2.4 (E and F) or normal at 3.4 (G and H). CVs were the highest at moderate noise levels near transition points when the CBI-2 to B34 conductance was low (E and F). CVs increased with noise levels when the CBI-2 to B34 conductance was normal (G and H), although the overall CVs were low compared with the data in (E) and (F). The red vertical dashed lines indicate the SDs used to generate data in Fig. 5 (M and N).

We hypothesized that these different findings were a reflection of the fact that under normal conditions, the connection between CBI-2 and B34 is stronger than that between CBI-10 and B34. To test this hypothesis, we determined the relationship between σ and variability in the CBI-10 circuit at a higher than normal CBI-10–B34 conductance (Fig. 6, C and D). In this situation, variability did increase when synaptic noise increased (although the highest CV was lower than the lowest CV observed when the CBI-10–B34 conductance was normal), similar to the behavior of CBI-2 circuit at its normal conductance (Fig. 6, G and H). In addition, we determined the relationship between σ and variability in the CBI-2 circuit at a lower than normal CBI-2–B34 conductance (Fig. 6, E and F). In this situation, after its maximum, variability decreased toward a plateau as noise increased, similar to the behavior of CBI-10 circuit at its normal conductance (Fig. 6, A and B). Thus, our results show that when the networks are near transition points (e.g., with relatively low conductance), they can behave in a strongly nonlinear manner and be very sensitive to noise. Consequently, a moderate level of noise actually promotes higher program variability than a higher level of noise.


We performed physiological and computational studies to determine mechanisms for the different degree of variability in the activity of a small neural network in Aplysia. Others studying variability have concentrated on a specific situation, one in which neurons receive balanced excitatory and inhibitory synaptic inputs (2, 5, 9, 1316, 33). This occurs in the cortex and other brain regions/circuits (9) and promotes the variability in spike train of single neurons. In this situation, variability is essentially generated postsynaptically when synaptic integration occurs, e.g., the combination of excitation and inhibition that a neuron receives results in fluctuations in its membrane potential.

In contrast, we study variability in a different situation, where neurons are primarily driven by excitatory input. In this situation, we showed that presynaptic factors (i.e., synaptic noise in the command-like neurons—pattern-generating interneuron connection) play an essential role in generating variability. Specifically, at the circuit level, we established a causal role for a pivotal circuit element (B34) in determining whether motor output is variable (Fig. 2, D to K) and showed that high or low motor variability in the Aplysia feedforward circuit is determined by activity levels of B34 (Fig. 2). If B34 is weakly activated [as is the case when network activity is driven by one command-like neuron (CBI-10)], then variability in network activity is high. If B34 is strongly activated [as is the case when network activity is driven by a second command-like neuron (CBI-2)], then variability in network activity is low. At the synapse level, synapses from both command-like neurons to B34 vary across trials, and the CBI-2 connection was stronger than the CBI-10 connection (Fig. 3). Putting the two together, variability is high when synapses are both weak and noisy (Figs. 4 and 5), i.e., when activity is triggered by CBI-10. Postsynaptic factors do, however, also appear to play a role in our network. We showed that B34 has a relatively high activation threshold [it has both a high spike threshold and low input resistance (fig. S3, C and D)]. A high activation threshold is presumably important since it allows “room” for synaptic integration to occur and different degrees of variability to be expressed.

One noteworthy aspect of our research is that we demonstrated that synaptic noise in conjunction with low synaptic strength can be a major source of motor variability in the situation where neurons are primarily driven by excitatory input. These situations are clearly common. For example, excitatory input appears to be prevalent in central pattern-generating networks where a half-center oscillator (34) provides excitatory inputs to follower motoneurons and is itself often driven by excitatory inputs from higher-order neurons. Central pattern generators (CPGs) underlie many invertebrate and vertebrate behaviors (1722) and may play a role in the operation of cortical networks as well (31).

Another noteworthy aspect of our work is that we performed parameter analyses in which we altered the amount of synaptic noise and demonstrated that the impact of synaptic noise on program variability depends on the amplitude of synaptic strength (Fig. 6). Specifically, as expected, when synaptic conductance is high, higher levels of noise do tend to promote higher motor variability, but the overall variability is quite low (even lower than the lowest variability when synaptic strength is low). In contrast, when synaptic strength is low and motor variability is at a high level, maximum variability is actually promoted by moderate rather than high levels of noise. Although other studies have not explicitly manipulated noise levels, they have shown that higher energy with which higher noise is often associated promotes more stable oscillations or less variability in molecular and cellular networks (35). The promotion of maximum motor variability at moderate noise levels is also akin to the broadly defined stochastic resonance that is present in a variety of nonlinear systems (11, 36), including the nervous system.

A further question of interest is whether variability in motor behavior is desirable. It has been suggested as an optimal behavioral strategy in the face of an uncertain environment (3, 5, 10, 12). We also showed that presynaptic noise underlying motor variability expands the working range of synaptic strengths by which a command-like neuron is capable of driving functional motor programs (Fig. 5, I to L). A similar role of noise has been observed in other circuits (37) and a molecular network (38). Moreover, in other systems, variability in neural activity can change during sensory perception (39), learning (6, 40), or development (7, 41). For example, songbird brains actively exploit variability to facilitate reinforcement learning of songs (42). Different activity states in two Caenorhabditis elegans interneurons correspond to reliable or variable chemotactic behavior (43).

Specific synaptic mechanisms underlying the variability in songbirds and C. elegans are largely unknown. Because songbirds use feedforward circuits (42), they may also use the mechanisms characterized here. There are correlative data suggesting that strengthening and pruning action–specific connections during song development reduce the sensitivity of motor control circuits to variable inputs and thereby make the motor output less variable (41). Although variability in the C. elegans network involves feedback connections between interneurons, chemical synapses between the interneurons are apparently important for variability. Thus, an interaction of synaptic noise and synaptic strength could play a role in this network as well. More generally, given that command-like neurons and associated feedforward circuits are present in many other animals (1719, 21, 22, 25) including vertebrates (20, 26, 27), we expect these findings to be of general interest. In particular, we predict that the synaptic mechanisms identified here will play a central role in neural networks in which circuit elements are driven primarily by excitatory synaptic inputs.


Subjects and electrophysiology

Experiments were performed on Aplysia californica (100 to 300 g) obtained from Marinus (Newport Beach, CA, USA). Aplysia are hermaphroditic (i.e., each animal has reproductive organs normally associated with both male and female sexes). Animals were maintained in circulating artificial seawater (ASW) at 14° to 16°C.

Intracellular recordings were made using single-barrel electrodes (5 to 10 megohms) filled with 0.6 M K2SO4 and 60 mM KCl. Intracellular signals were acquired using an Axoclamp 2B or 900A amplifier (Molecular Devices), a Neuroprobe amplifier model 1600 (A-M Systems), or a Getting model 5A amplifier. A Grass model S88 stimulator was used for stimulation. Extracellular signals were acquired from polyethylene suction electrodes using a differential AC amplifier (model 1700; A-M Systems). Recordings were made in ASW [460 mM NaCl, 10 mM KCl, 55 mM MgCl2, 11 mM CaCl2, and 10 mM Hepes buffer (pH = 7.6)] unless otherwise indicated. All chemicals were purchased from Sigma-Aldrich (St. Louis, MO).

Cell identification

Most neurons studied in this report (fig. S1) have been identified using criteria described previously. These include the cerebral neuron CBI-2 (18, 19, 23, 24, 30, 44); buccal pattern-generating interneurons B34, B63 (19, 32, 44), and B64 (45); and buccal motoneuron B8 (18). CBI-10 (Fig. 1B and fig. S1), like other CBIs, projects its axon to the buccal ganglion through the CBC and was named previously (46) on the basis of backfills of the CBC. The morphological, physiological, and network properties of CBI-10 were described in the present paper. To reveal the morphology of CBI-10, an electrode containing 3% 5(6)-carboxyfluorescein dye in 0.1 M potassium citrate was used to ionophoretically inject dye (currents: −5 to −8 nA for 10 to 20 min). A fluorescence microscope (Nikon or Olympus) was used to view and photograph the ganglion. Three to six pictures were taken through the Z axis to focus on the soma or the processes at different depths, and these pictures were overlaid using Photoshop (Auto-Blend Layers/Stack Images).

Semi-intact preparations and feeding behavior

We used a semi-intact preparation in which activity of the CBIs can be recorded while the animals eat (19). Animals (100 to 200 g) were anaesthetized by injections of 333 mM MgCl2 (~33% of body weight). An incision was made from the dorsal side. The gut was separated from the rest of the animal at the level of the esophagus. The head, including tentacles, rhinophores, and buccal mass, together with the CNS (cerebral ganglion, buccal ganglion, and pleural-pedal ganglion), was cut from the rest of the animal at the level of the rhinophores (dorsally) and slightly caudal to the beginning of the foot (ventrally).

The head structure remained innervated by the upper labial nerve and the anterior tentacular nerve of the cerebral ganglion, and the buccal mass remained innervated by all buccal nerves (fig. S1). The preparation was set in a two-chamber dish lined with Sylgard (Dow Corning). The intact head (including the buccal mass) was situated in the larger deeper chamber, while the CNS was pinned in the smaller, shallower chamber. The sheath over the E cluster, which contains CBI-10, was removed.

The foot artery was ligated. The buccal artery was cannulated to allow continuous perfusion of the head and buccal mass with cooled fresh ASW (plus 10 mM glucose) at ~2 ml/min using a bubble separator in the perfusion line (inlet). The buccal mass pressure was monitored with a pressure transducer whose probe was placed in the bubble separator. Buccal mass pressure, a noninvasive means of measurement, has been shown to correspond well to feeding movements (19). The preparation was maintained at 14° to 16°C.

Feeding behavior was elicited by touching the mouth with cut pieces of dried seaweed held by a pair of forceps. Pieces of seaweed were approximately 4 mm by 7 mm. A typical feeding sequence consisted of biting-swallowing responses that were initiated when the seaweed was successfully grasped. Completion of feeding was verified by the visual observation of the seaweed moving out of the cut end of the esophagus.

Isolated CNS preparations and motor programs

Electrophysiological recordings from CNS preparations (including the cerebral and buccal ganglia) were performed as described previously (19, 29). Animals were anesthetized by injection of 333 mM isotonic MgCl2 (~50% of body volume), and the cerebral and buccal ganglia were dissected out. Ganglia were desheathed, transferred to a recording chamber (lined with Sylgard) containing ~1.5 ml of ASW, continuously perfused at 0.3 ml/min, and maintained at 14° to 17°C. To suppress polysynaptic pathways, a high divalent cation saline (HiDi) was used containing the following: 368 mM NaCl, 10 mM KCl, 13.8 mM CaCl2, 101 mM MgCl2, and 10 mM Hepes (pH 7.6). This saline does not alter PSP amplitude (30). To block chemical synaptic connections, we used a 0 Ca2+ saline containing the following: 368 mM NaCl, 10 mM KCl, 101 mM MgCl2, and 10 mM Hepes (pH 7.6). This solution had no obvious effects on electrical connections.

Hexamethonium chloride was dissolved in HiDi immediately before each application. The HiDi with the drug was perfused into the recording chamber. During all pharmacological and 0 Ca2+ experiments, we made sure that the membrane potential throughout experiments remained the same by applying appropriate hyperpolarizing or depolarizing currents. The cell hyperpolarization/depolarization was performed using single-electrode current-clamp technique, and care was taken to correctly balance the electrode resistance.

Motor programs were elicited via intracellular stimulation of CBI-2 or CBI-10 using square wave current pulses at various frequencies as specified. Each current pulse was set to trigger a single action potential. CBI stimulation was manually terminated after the protraction phase. Each cycle of CBI-2– and CBI-10–elicited programs has primarily two phases, i.e., the protraction phase precedes the retraction phase. These two phases mediate the protraction-retraction movements of the radula. For the most part, the protraction phase was monitored using the I2 nerve, which contains the axons of the radula protraction motoneurons B61/62 (18, 23, 44). The retraction phase was monitored either by the appearance of activity in the retraction interneuron B64 (45) or by the hyperpolarization of protraction interneurons, such as B34, which were inhibited by B64. In some cases, retraction was also monitored by B8 depolarization after protraction termination.

To determine the effects of depolarization and hyperpolarization of B34 on motor programs elicited by CBI-2 or CBI-10, we continuously stimulated a CBI to induce 5 to 6 cycles of programs for each trial and waited for 9 min between each trial. To characterize synaptic noise, we stimulated CBI-2 or CBI-10 at 10 Hz for 2 s for 10 trials and waited 3 min between each trial. The stimulation duration was set at 2 s because this allows the expression of synaptic facilitation without evoking too many spikes. The amplitude of the last EPSP at each trial was quantified, and the CVs were computed for the 10 EPSPs.

Mathematical framework for the Aplysia feeding network–like model: Neuronal elements and synaptic connections

In constructing a compartmental model of the Aplysia feeding circuit (Fig. 1A), we considered two major issues: (i) the cyclic or oscillatory activity generated by pattern generating elements and (ii) the feedforward inputs originating from command-like neurons, CBIs. The cyclic, biphasic protraction-retraction sequence of the Aplysia feeding behavior is primarily generated by the protraction neurons B63/B31-32 (44, 47) and the retraction interneuron B64 (45). Although there is modeling work that has sought to determine how the electrically coupled B63/B31-32 complex may mediate protraction (47), no modeling on the generation of the protraction-retraction sequence in Aplysia is available for others to easily replicate. In contrast, a computational study (21) on a feeding circuit of a closely related gastropod, Lymnaea stagnalis, simulated the generation of an analogous protraction-retraction sequence that is mediated by negative feedback between the protraction neuron N1m and retraction neuron N2v (fig. S4C). The firing patterns of N1m and N2v during motor programs are quite similar to those of the homologous neurons, B63 and B64, in Aplysia. Thus, instead of building an oscillatory network composed of the B63/B31-32 complex and B64 by identifying and tuning parameters, we adapted the basic modeling framework of N1m and N2v and their neuronal and synaptic parameters to model B63 and B64, respectively. Some parameters were modified to recapitulate Aplysia recordings, as described later in the “Parameter setting in the model” section.

Specifically, B63 and B64 were modeled with two compartments (fig. S4A). One compartment represents the soma, which is nonspiking but generates a plateau potential, and the other compartment represents the axon, which is spiking but does not generate a plateau potential. In some ways, the two-compartment model of B63, with less complexity in the model, could be considered a simplified version of the B63/B31-32 complex in Aplysia because B63 is spiking, whereas B31-32 somata are nonspiking but contain a plateau potential (47). B64 has also been shown to be a plateau-generating interneuron (45).

This is a reasonable approach because the circuit above is only used to generate the basic cyclic activity in a pattern-generating network. The primary focus of the modeling was to examine whether differences in CBI input to the pattern-generating network determine variability. It is notable that whereas the Lymnaea network (21) was modeled with graded chemical synaptic transmission, we only used spike-mediated synaptic transmission (see the “Chemical synapses” section below), as there is no evidence for graded transmission in Aplysia. In addition, we added a relatively weak inhibitory synapse from B63 to B64 (fig. S4A and table S3) because of its presence in Aplysia (44). No synaptic inhibition from N1m to N2v exists in Lymnaea (fig. S4C).

On top of the above basic framework, we added three neurons, CBI-10 or CBI-2, and B34, to complete the feedforward circuit, and used this network as a minimal model (fig. S4A) for the Aplysia feeding circuit (Fig. 1A). In the model, the CBIs act as command-like neurons. Protraction interneuron B34 is a feedforward node. These cells were all modeled with a single compartment, representing the soma. Note that CBI-10 and CBI-2 have the same intrinsic parameters but exert different synaptic strengths on B34 and B63. The parameters for CBI input were all obtained from Aplysia data (see also the next section).

All compartments were modeled as followsCmdVdt=IinjIleakIxIecIcs(1)where V is the membrane potential of the soma or axon, and Cm represents membrane capacitance. Iinj is the external current and only applied to CBIs to generate network activity. Ileak = gL × (VEL) is the leakage current, where EL is the leakage reversal potential.

Hodgkin and Huxley ion currents. Ix in Eq. 1 represents the summation of all ion channel currents in a compartment, which all take the Hodgkin and Huxley formulationIx=gx,max·pxkx·qxlx·(VEx)(2)where x ∈ {Ion Channels}, gx, max is the maximal conductance, and Ex is the reversal potential of ion channel. px and qx are activation and deactivation variables, and their exponents kx and lx take values from the set {0,1,2,3,4}. Both px and qx follow the first-order relaxation kinetics{dpxdt=px,pxτpxdqxdt=qx,qxτqx(3)where px,∞ and qx,∞ are the steady-state values and τpx and τqx are relaxation time constants. However, it is conventional to use mNa, hNa, and nK rather than px and qx to label currents in the spiking compartment (see table S2 for details).

Connections between somatic and axonal compartments. The somatic and axonal compartments of either B63 or B64 are connected with electrical coupling (fig. S4A). Iec, post = gec × (VpostVpre) represents the current between the soma (S) and axon (A), where gec is the electrical coupling conductance (table S2). The currents flowing to soma and axon are as follows{Iec,S=gec,S×(VSVA)Iec,A=gec,A×(VAVS)(4)

Chemical synapses. Ics in Eq. 1 denotes summation of chemical synaptic currents and is given byIcs=i gcs,i·si·(VEcs,i)(5)where gcs,i is the maximum conductance and Ecs,i is the reversal potential of the ith synapse (Ecs = 0 mV for excitatory synapse, and Ecs = − 90 mV for the inhibitory one). The gating variable si follows the second-order kinetics{dridt=r,iriτcs,idsidt=risi τcs,i(6)where τcs,i is the relaxation time constant (table S3). si represents the activation of postsynaptic receptors, while ri represents the activation of presynaptic transmitter release. In this spike-mediated model that we used, the transmitter variable r∞,i equals 1 during a presynaptic spike (i.e., when the membrane potential of the presynaptic neuron crosses a spike threshold as specified in table S3) and equals 0 otherwise.

Other critical synaptic and network properties in the model

As stated earlier, some of the neuronal and synaptic parameters were adapted from the Lymnaea model (21). Here, we present more details about modeling, which are critical to the Aplysia network activity and variability generation. These aspects of modeling were largely absent in the Lymnaea model. No variability in motor programs was apparent in the Lymnaea model.

Synaptic plasticity and presynaptic inhibition. At synapses that manifest plasticity, i.e., the excitatory synapses from CBIs to B34 (Eq. 6), r∞,i = 1 should be replaced with a variable characterizing synaptic plasticity, ε, whose kinetics obeydεdt={ε2(εmaxε)τε,1, during a presynaptic spike1ετε,2,otherwise(7)(8)where τε,1 and τε,2 are two time constants.

To mimic the experimental data on the B34 EPSPs when CBI-2 and CBI-10 were stimulated at different frequencies (5 to 15 Hz), we found it necessary to modify the common synaptic-plasticity dynamics described previously (47, 48) by limiting the maximal amount of transmitter release. That is, upon arrival of a presynaptic spike, ε is set to 1, then rises up to εmax, and lastly relaxes to 1 when there is no presynaptic spiking. For simplicity, only the CBIs to B34 synapses are considered plastic, while the CBIs to B63 ones are not (table S1).

It was also found (49) that B64 presynaptically inhibits CBI-2 synaptic output to its targets in the buccal ganglia, e.g., B34. We observed a similar phenomenon for CBI-10. That is, the presynaptic inhibition of CBIs-B34 synapses is switched on whenever B64 is active. Thus, Eq. 7 should be replaced with the followingdεdt=1ετε,3,during presynaptic inhibition(9)

The specific values for the time constants in this section are presented in table S3.

Synaptic noise. There are several forms of synaptic noise (11), including presynaptic noise (33, 50), postsynaptic noise, and dynamic noise. We took the presynaptic noise into consideration (Fig. 3, E to G). Because of the ultrasensitivity in input-output curve of EPSPs, small fluctuations in presynaptic depolarization likely cause a relatively large effect on the amplitude of EPSPs. However, in our model, the membrane potentials of CBIs are deterministic and have no fluctuations when CBIs are depolarized. Thus, to simulate the presynaptic noise, we randomly varied the transmitter variable r∞,i in Eq. 6 for each CBI spike, which is assumed to be the direct consequence of fluctuations in presynaptic depolarization.

In a biological network, every synapse likely contains some noise. However, to focus our study on examining whether differences in CBI input to the pattern-generating network determine variability, only the synapses from CBIs to B34 and B63 (table S1) were considered noisy. Accordingly, the transmitter variable in Eq. 6 was modified into r∞, CBIs − B34 = ε · N(1, σ2) for CBIs-B34 synapses and r∞, CBIs − B63 = 1 · N(1, σ2) for CBIs-B63 synapses, where N(1, σ2) is a Gaussian random variable with mean of 1 and variance of σ2. The transmitter variables were only different from spike to spike but remained the same during a single presynaptic spike.

Note that as demonstrated in Results, this presynaptic noise modeled variability in EPSP amplitude very well (Fig. 3). However, we cannot exclude the possibility that there is a small postsynaptic contribution in the biological network.

Background input noise. To determine whether other forms of noise could play some role, background input noise was also considered (fig. S5). For the neurons with background noise, specifically B34, an external current is given by I = gb × (VEb) × Sb, where gb is the conductance and Eb is the reversal potential. Sb obeys the first-order kineticsdSbdt=Sbτb+k δ(ttk)(10)where {tk} is a Poisson spike train with an average rate of 10 Hz for either excitatory or inhibitory input or 2 Hz for concurrent excitatory and inhibitory inputs (fig. S5) (51). This form of synaptic noise usually originates from spiking of presynaptic neurons that are outside of the model network (fig. S5).

Parameter analyses in the model

To systematically explore the effects of synaptic strength and noise in the network, we uniformly selected 100 conductance values for both CBI-10–B34 synapse (range, 0.5 to 1.3) and CBI-2–B34 synapse (range, 2.0 to 4.0). For motor programs generated in the model, protraction is defined by bursting activity in B63, whereas retraction is defined by bursting activity in B64. On the basis of B63/B34 activity, we classified the network behavior as no programs [neither B63 or B34 is active; orange areas in Figs. 5 (I to L) and 6 (A and B)], programs by B63 only [purple areas in Fig. 5 (M to P)], or normal programs by both B63 and B34 [white areas in Figs. 5 (I to P) and 6]. The transition between different states is usually defined as bifurcation points, specifically, Hopf bifurcation in our model.

We also varied the level of synaptic noise, i.e., SD of Gaussian noise, σ. In Fig. 6, σ was uniformly distributed in the range of 0.08 to 0.3, and 100 values were taken. In fig. S7, σ was uniformly distributed in the range of 0.0 to 0.6, and 300 values were taken. We performed 10 simulations for each noise level. With given conductance and noise levels, we calculated the CVs for protraction duration, duty cycle, and B34 ISIs based on 10 simulation trials. Overall, the results from parameter analyses support the robustness of the model.

Parameter setting in the model

Summary information on synaptic properties is given in table S1. Detailed equations and parameters for each neuron are shown in tables S2. Parameters for electrical coupling between the soma and the axon of either B63 or B64 are shown in table S2 and for chemical synapses between network elements in table S3. Notably, all the currents and membrane capacitance Cm were multiplied with 1/gL for normalization.

As mentioned earlier, some parameters were adapted from the feeding CPG model of Lymnaea (21), which is similar to the feeding CPG of Aplysia. Specifically, B63 and B64 separately correspond to N1m and N2v in Lymnaea. We adjusted several parameters to mimic electrophysiological recordings in Aplysia. In particular, we added a relatively weak inhibitory synapse from B63 to B64 (fig. S4A and table S3), as described earlier. We set the synaptic conductance of both the B34-B64 (inhibitory) and B63-B64 (both excitatory and inhibitory) synapses so that the protraction duration was about 17 s. τq for the B64 soma was altered so that the retraction duration was about 7 s rather than 1 s in Lymnaea. Parameters for B34 were adapted from the axonal compartment of B63. However, B34 K-channel time constant τn was increased from 1.6 to 50.6, and the K-channel conductance was changed from 90 to 300 so that B34 afterpotential matched experimental recordings.

On the basis of our physiological experiments, what appears to be most critical to the main subject of the study, i.e., differences in variability of programs, is primarily due to the differences in synaptic connections from the two command-like neurons (CBI-10 and CBI-2) to B34. However, these three neurons are not present in the Lymnaea model (21). Thus, we chose synaptic parameters (table S3) to ensure that the EPSPs from the CBIs to B34 and B63 mimic electrophysiological recordings. Overall, the excitatory synapses from CBI-10 were weaker than those from CBI-2 (table S1). In addition, presynaptic noise from CBIs to B34 and B63, another focus of the study, was modeled (table S1).

Simulation method

We numerically integrated 35 differential equations with Python 3.6, using the Numpy and Scipy packages. We applied the scipy.integrate.ode function with parameter “dopri5,” which used an explicit fourth-order Runge-Kutta method with a time step of 0.2 ms. All simulations were run on a Linux system with 2.40-GHz Xeon E5 at the High Performance Computing Center of Nanjing University.

Data and statistical analyses

Electrophysiological recordings were digitized online with Axoscope (Molecular Devices) and plotted with CorelDRAW (Corel). Bar graphs and scatter diagrams were plotted using Prism software (GraphPad Prism 7). Data are expressed as means ± SEM. All experimental data were taken from individual animals or preparations, and n refers to the number of preparations.

Statistical tests were performed as appropriate using GraphPad Prism software. They included Student’s t test, one-way analysis of variance (ANOVA), as appropriate. Data that showed significant effects in ANOVA were further analyzed in individual comparisons with Bonferroni’s correction. Correlations between two sets of data were performed to obtain Pearson’s correlation coefficient (r). The correlation coefficient ranges from −1 to 1, with −1 indicating a perfect inverse correlation and 1 indicating a perfect positive correlation. In all statistical tests, effects were considered statistically significant when P < 0.05.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank J. Byrne for comments on an earlier version of the manuscript and J. Welsh, C.-Y. Zhang, and K. Zen for discussions. The numerical calculations in this paper have been performed on the computing facilities in the High Performance Computing Center (HPCC) of Nanjing University. Funding: This work was supported by the National Natural Science Foundation of China (grants 31671097, 31861143036, 31371104, J1103512, and J1210026) and the NIH (grants NS066587 and NS070583). Author contributions: J.J., E.C.C., and K.R.W. conceived the physiological aspects of the project. J.J., G.Z., T.-T.C., and S.-A.C. performed the experiments. W.-D.Y., Z.-W.L., K.Y., and Z.Y. helped perform the experiments. J.J., F.L., and T.W. conceived the modeling aspects of the project. T.W. and F.Y. performed the modeling studies. J.J., G.Z., T.W., T.-T.C., Y.-Y.X., S.-Q.G., K.Y., and S.-A.C. analyzed the data. J.J. and E.C.C. wrote the paper. F.L., K.R.W., G.Z., and T.W. edited the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The custom computer codes (Python) used for the simulation of the Aplysia feeding circuit are available upon request (T.W., dz1622031{at} Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article