Computational model of the distributed representation of operant reward memory: combinatoric engagement of intrinsic and synaptic plasticity mechanisms.

Operant reward learning of feeding behavior in Aplysia increases the frequency and regularity of biting, as well as biases buccal motor patterns (BMPs) toward ingestion-like BMPs (iBMPs). The engram underlying this memory comprises cells that are part of a central pattern generating (CPG) circuit and includes increases in the intrinsic excitability of identified cells B30, B51, B63, and B65, and increases in B63-B30 and B63-B65 electrical synaptic coupling. To examine the ways in which sites of plasticity (individually and in combination) contribute to memory expression, a model of the CPG was developed. The model included conductance-based descriptions of cells CBI-2, B4, B8, B20, B30, B31, B34, B40, B51, B52, B63, B64, and B65, and their synaptic connections. The model generated patterned activity that resembled physiological BMPs, and implementation of the engram reproduced increases in frequency, regularity, and bias. Combined enhancement of B30, B63, and B65 excitabilities increased BMP frequency and regularity, but not bias toward iBMPs. Individually, B30 increased regularity and bias, B51 increased bias, B63 increased frequency, and B65 decreased all three BMP features. Combined synaptic plasticity contributed primarily to regularity, but also to frequency and bias. B63-B30 coupling contributed to regularity and bias, and B63-B65 coupling contributed to all BMP features. Each site of plasticity altered multiple BMP features simultaneously. Moreover, plasticity loci exhibited mutual dependence and synergism. These results indicate that the memory for operant reward learning emerged from the combinatoric engagement of multiple sites of plasticity.

Memories are believed to be encoded as sets of changes at many sites within neuronal networks, often referred to as engrams (Semon 1904; Josselyn et al. 2015; Josselyn and Frankland 2018). Indeed, multiple loci of plasticity appear to be a common feature of memory in nervous systems (e.g., Cleary et al. 1998; Hammer and Menzel 1998; Hansel et al. 2001; Crow and Tian 2006; Strube-Bloss et al. 2011; Gao et al. 2012; Kalmbach and Mauk 2012; Mayford et al. 2012; Tsien et al. 2013; Jörntell 2016; Richards and Lillicrap 2019). However, how individual sites of plasticity contribute to memory expression remains unknown. For example, which biophysical properties of which neurons contribute to a given behavioral feature? Does each site contribute to one or multiple behavioral features? How do multiple loci interact to allow expression of a memory? Are the contributions of multiple loci of plasticity purely additive, or are they synergistic? Despite significant advances in the tools available to identify as well as activate or suppress engrams (Tonegawa et al. 2015), many of the manipulations required to address these questions are currently not experimentally feasible. Thus, we used a biologically realistic computational model to examine the contributions of previously identified sites of learning-induced plasticity, individually and in combination, to specific behavioral features altered by operant reward learning.

The present study developed a conductance-based model of a central pattern generator (CPG) that mediates feeding in Aplysia. The feeding CPG generates buccal motor patterns (BMPs), which are patterns of neural activity that mediate rhythmic movements of a food grasper during feeding. These movements consist of a forward motion (protraction) followed by a backward motion (retraction). The timing of a third motion (closure) relative to this sequence is one of the main factors determining the type of behavior emitted. BMPs have phases of activity corresponding to protraction, retraction, and closure. At least two types of BMPs are expressed and, in turn, these BMPs mediate different behaviors, such as ingestion (iBMP) and rejection (rBMP; Weiss et al. 1986; Morton and Chiel 1993a,b; Hurwitz et al. 1996; Evans and Cropper 1998; Kabotyanski et al. 2000; Jing et al. 2004; Nargeot et al. 2007; McManus et al. 2012). The connectome of the CPG is at least partly characterized as are the biophysical properties of various cells, synapses, and sites of learning-induced plasticity (Jing and Weiss 2001, 2002, 2005; Elliott and Susswein 2002; Cropper et al. 2004, 2019; Baxter and Byrne 2006; Wentzell et al. 2009; Mozzachiodi and Byrne 2010; Nargeot and Simmers 2011). For example, operant reward learning increases the frequency and regularity of ingestions in vivo and in vitro, and biases activity in the CPG toward iBMPs (Brembs et al. 2002; Nargeot et al. 1997, 2007; Mozzachiodi et al. 2008). Correlates of operant reward learning include an increase in the excitability of neurons B30, B51, B63, and B65, and increases in electrical coupling among B30, B63, and B65 (Brembs et al. 2002; Nargeot et al. 1999a,b, 2009; Sieling et al. 2014). To examine the relative contributions that these sites of plasticity make to altering motor output, a computational model of the circuit was developed that included identified cells CBI-2, B4, B8, B20, B30, B31, B34, B40, B51, B52, B63, B64, and B65. Modeling a well-described system allowed us to bridge the gap between changes in biophysical properties and their ultimate effects on the features of motor output.

Results

Simulation of biologically relevant patterns of neural activity

As a first step in examining the ways in which distributed sites of learning-induced plasticity modify rhythmic activity, a conductance-based, reduced model of the Aplysia feeding CPG was developed (Fig. 1). The model consists of a subset of the neurons in the feeding CPG that was sufficient to generate activity that resembles physiological BMPs elicited in isolated ganglia. The current model is an expansion of our previous efforts in modeling this system (Kabotyanski et al. 1994; Ziv et al. 1994; Hayes et al. 2005; Cataldo et al. 2006, 2009; Baxter et al. 2010).

Figure 1.

The connectome. The connectome of the model represents the topology of monosynaptic connections among cells. Cells highlighted in orange are generally active during the protraction phase of a BMP, whereas cells highlighted in blue are generally active during the retraction phase of a BMP (e.g., Kabotyanski et al. 1998; Jing and Weiss 2001; Jing et al. 2003, 2004; Cropper et al. 2004; Nargeot and Simmers 2012). Activity in closure motor neuron B8 (highlighted in green) shifts between the protraction phase for rBMPs and the retraction phase for iBMPs (see Fig. 3). Note, the network included one hypothetical excitatory connection from B63 to B64 (not shown) because our previous modeling studies indicate that some additional excitatory drive onto B64 is necessary to elicit the retraction phase (Cataldo et al. 2006; see Materials and Methods for details).

Figure 3.

Methods for analyzing BMPs. Four measures of network activity were used. First, the overall number of BMPs was counted. Three simulated BMPs are illustrated here. Second, BMPs were classified as being iBMPs or other pattern types, such as the depicted rBMPs (see Materials and Methods). Third, the mean, SD, and CV of the inter-burst intervals (IBIs) were calculated. The beginning of a burst was defined as the first spike in B31a, which matches methods used to identify bursts in empirical studies (e.g., Nargeot et al. 1997). Finally, activity maps were generated by counting the number of spikes in each cell during 20 nonoverlapping 1 sec time bins that ranged ±10 sec from the terminus of activity in B31a.

Figure 2A illustrates the simulated activity observed in all 13 cells upon stimulation of CBI-2. Because both axonal and somatic compartments of B31, B51, and B64 are depicted, a total of 16 traces are shown. CBI-2 is a command-like neuron that receives sensory input associated with the presence of food and that excites many CPG elements. Tonic presentation of food elicits rhythmic feeding in vivo and, similarly, tonic stimulation of CBI-2 evokes rhythmic fictive feeding in vitro (Rosen et al. 1991; Hurwitz et al. 2005). Thus, tonic stimulation of CBI-2 was used to elicit rhythmic activity in the model. As indicated in Figure 2, the model simulated rhythmic patterns of activity. The relative phases of cellular activities were similar to empirical observations (e.g., Church and Lloyd 1994; Kabotyanski et al. 1998; Jing and Weiss 2001; Nargeot et al. 2002; Jing et al. 2003, 2004; Cropper et al. 2004; Nargeot and Simmers 2011). Thus, the model reproduced the biphasic properties of BMPs, in which activity in protraction-phase neurons is always followed by activity in retraction-phase neurons. Moreover, the overall frequency of rhythmic activity was similar to empirical observations (Church and Lloyd 1994; Jing and Weiss 2001; Hurwitz et al. 2003; Sieling et al. 2014).

Figure 2.

Simulation of BMPs. Activity was elicited by injecting a sustained 1.9 nA depolarizing current into command-like neuron CBI-2. The color code of the voltage traces matches that in Figure 1. The variability among the simulated BMPs resulted from the noise that was included in all simulations. Here, both the somatic and axonal compartments of B31, B51, and B64 (e.g., B31s and B31a, respectively) are illustrated. (A) Control simulation. (B) Simulated activity after incorporating neuronal correlates of operant conditioning. See Figure 4 and text for details related to the implementation of these correlates.

Figure 4.

Implementing neural correlates of operant conditioning. Empirical studies indicate that the memory engram of operant reward learning is encoded as decreases in the input conductance of B30, B63, and B65 (A); increases in the electrical coupling between B63 and B30, and between B63 and B65 (B); increases in the excitability of B30, B63, and B65 (C); and a decrease in the input conductance and an increase in the excitability of B51 (D). In each panel, the black trace represents the control simulations, whereas the red trace represents the simulated neuronal correlate of the memory engram following operant conditioning. (A) The input conductances of B63 (A1), B30 (A2), and B65 (A3) were measured by injecting a −0.5 nA, 2-sec duration current pulse (indicated by bar below trace) into each individual cell. (B) Coupling coefficients were measured by injecting a −1 nA, 2-sec duration current pulse into B63 (left trace) and measuring the voltage deflections in B30 and B65 (traces to right; note change in scale among panels). (C) The excitabilities of B63 (C1), B30 (C2), and B65 (C3) were measured by injecting a 2-sec duration depolarizing current pulse into each cell. The magnitudes of currents were adjusted to be subthreshold in the control simulation (0.4 nA in B63, 0.73 nA in B30, and 0.74 nA in B65). Identical pulses were injected after incorporating neuronal correlates of the memory engram. (D) The input conductance (D1) and excitability (D2) of B51 were measured by injecting 1-sec duration current pulses; either a −0.5 nA pulse to measure input conductance, or a 0.25 nA pulse to measure excitability. In all examples, the cells are embedded within the connectome (see Fig. 1).

The model simulated intrinsic features of cellular activity, such as the plateau-like potential in B31 and the role of autaptic transmission in its maintenance (Hurwitz et al. 1994, 2008; Saada et al. 2009; Baxter et al. 2010). The model also reproduced properties of cellular activity that emerge due to the specific pattern of connectivity in which the neurons are embedded. For example, the model reproduced observations that activity in CBI-2 becomes rhythmic despite its tonic stimulation, and that the burst of activity in CBI-2 occurs during the protraction phase (Hurwitz et al. 2005). In addition, B52 was active at the end of each BMP and functioned to terminate activity in cells active during the retraction phase (Nargeot et al. 2002; Shetreat-Klein and Cropper 2004).

The model generated a variety of BMP types similar to that observed empirically. The simulated patterns were classified as either ingestion-like (iBMPs) or other types of pattern, which included rejection-like (rBMPs) and intermediate patterns. As illustrated in Figure 3, patterns were classified based on the relative overlap of activity in B8 with the protraction versus the retraction phases (see Materials and Methods). The examples in Figures 2, 3 indicate that the model simulated various patterns of activity similar to fictive behaviors (e.g., Jing and Weiss 2001; Jing et al. 2003, 2004), although not all differences between pattern types were captured (e.g., B4 exhibits higher activity during rBMPs, which was not observed in the model). Thus, the model appeared to be a reasonable abstraction of the feeding CPG and was used to examine the ways in which the distributed representation of an engram alters network activity.

Simulation of learning-induced changes in neural activity

At least three sets of cellular correlates of memory have been characterized following operant reward learning: (i) a decrease in the input conductance and increase in the excitability of pattern-initiating neurons B30, B63, and B65 (Nargeot et al. 2009, Sieling et al. 2014); (ii) a decrease in the input conductance and increase in the excitability of B51, a pattern-selecting neuron (Nargeot et al. 1999b; Brembs et al. 2002; Lorenzetti et al. 2008; Mozzachiodi et al. 2008); (iii) an increase in the electrical coupling between B63 and B30, and between B63 and B65 (Nargeot et al. 2009; Sieling et al. 2014). Figure 4 illustrates the implementation of these sites of plasticity.

Learning-induced changes in the input conductance and excitability of B30, B63, and B65 were simulated by modifying leakage conductances in these cells (Fig. 4A,C). Changes in the B63–B30 and B63–B65 electrical synapses were stimulated by increasing their coupling conductances (Fig. 4B). Finally, the learning-induced change in B51 was simulated by decreasing the leakage conductance of the cell (Fig. 4D; see Materials and Methods for details on specific changes made to each cell and connection). We first tested whether there is, in principle, at least one set of parameters such that changes to known sites of plasticity are sufficient to recapitulate the main features of operant reward learning. Indeed, adjusting the magnitude of the modifications described above by trial and error was sufficient to yield increases in motor pattern frequency, regularity, and bias toward iBMPs (see below), which are changes induced by operant reward learning in vivo and in vitro (Nargeot et al. 1997, 1999a, 2007, 2009; Brembs et al. 2002; Mozzachiodi et al. 2008; Sieling et al. 2014).

Computational model of operant reward memory recapitulates empirical observations

The implementation of simulated learning-induced changes led to a substantial reconfiguration of rhythmic activity generated by the simulated network (Fig. 2B). A more detailed analysis shows that this reconfiguration reproduced features of learning that have been observed empirically (Fig. 5). Empirical studies indicate that operant reward learning increases the number of BMPs and biases activity in the CPG toward expressing iBMPs (for reviews, see Baxter and Byrne 2006; Mozzachiodi and Byrne 2010; Nargeot and Simmers 2011). The control and operant memory conditions were simulated 20 times each during a 6-min stimulus to CBI-2 and the resultant BMPs were classified as being either iBMPs or other patterns (Fig. 5A). Relative to the control simulations, the operant learning simulations had an increase in both the overall rate of BMPs (from 1.71 ± 0.08 to 3.78 ± 0.13 BMPs/min; mean ± standard error) and in the rate of iBMPs (from 0.43 ± 0.04 to 1.93 ± 0.07 BMPs/min). Whereas in the control simulation ∼26% of the patterns were ingestion-like, after implementation of the operant memory that ratio increased to 51%. The model, therefore, reproduced both the increase in motor pattern rate and the shift in pattern type bias that are characteristic of operant learning.

Figure 5.

The simulated neuronal correlates of memory altered the functional properties of the network. (A) Memory-induced changes in the total number of BMPs and numbers of iBMPs versus other patterns. The number and types of BMPs were assessed during 20 simulations of 6 min of activity. In the control simulations, overall patterns occurred at a rate of 1.71 ± 0.08 BMPs/min (mean ± standard error), and iBMPs at a rate of 0.43 ± 0.04 BMPs/min. In contrast, the rates during the operant memory simulations were 3.78 ± 0.13 BMPs/min for overall patterns and 1.93 ± 0.07 BMPs/min for iBMPs. Each error bar denotes the standard error of the plot below it. (B) Memory-induced increase in the regularity of rhythmic activity. The inter-burst intervals (IBIs) were measured in control simulations and simulations following the incorporation of neuronal correlates of memory, and the CV was calculated in each case for 15 simulations that had at least 10 BMPs. Changes in CV indicate changes in the regularity with which BMPs occur. Decreases in CV represent greater regularity, whereas increases in CV represent less regularity. Histograms and fitted Gaussian curves show the distribution of all 135 IBIs in each condition (B1), and bar plots show corresponding CVs (B2). (C) Memory-induced reconfiguration of the activity map. Activity maps were generated by averaging the activity in 100 BMPs from control simulations (C1) and from simulations with the total ensemble of operant changes (C2). No distinction was made regarding the classifications of the individual BMPs. Note that, because patterns can occur in close proximity, activity maps may capture spiking from adjacent patterns (e.g., activity in protraction neurons toward the end of retraction). To highlight changes in activity that occurred after learning, the operant memory activity map was subtracted from the control map (C3). This subtraction, which is referred to as a “dynamic memory map,” revealed the changes in average spiking for each cell in the network, regardless of whether they were direct targets of modulation.

Another consequence of operant reward learning is an increase in the regularity with which BMPs are expressed (Nargeot et al. 2007, 2009). In naïve animals (or in vitro preparations), biting (or fictive feeding) occurs sporadically. Learning regularizes this behavior (or network activity). To evaluate the ways in which combined sites of plasticity affect the regularity of network activity, inter-burst intervals (IBIs) were measured. In each condition, 15 simulations that had at least 10 BMPs were analyzed. The IBIs between 10 BMPs were measured in each simulation and a histogram of the distribution of IBIs within each group of simulations was constructed (Fig. 5B1). Finally, the coefficient of variation (CV), which is equal to the standard deviation (SD) divided by the mean IBI, was calculated for each simulation. A decrease in the CV represents an increase in the regularity of network activity. The CV in control simulations was 0.167 ± 0.018, whereas implementing the full ensemble of cellular correlates reduced the CV to 0.059 ± 0.002 (Fig. 5B). These results indicate that the combination of learning-induced changes in B30, B51, B63, and B65 was sufficient to increase the regularity of rhythmic activity generated by the simulated network. Therefore, there was at least one set of parameters for which implementation of the combined sites of plasticity was sufficient to reproduce the main features of operant reward learning in Aplysia (i.e., changes in frequency, bias, and regularity).

Activity maps reveal indirect effects of the engram

To examine the ways in which activity in all 13 cells of the CPG was altered by implementation of the operant memory, activity maps were obtained by averaging 100 BMPs from control simulations (Fig. 5C1) and the same number from operant memory simulations (Fig. 5C2). The effects of plasticity were revealed by subtracting the control activity map from the operant activity map. This subtraction produced a “dynamic memory map” (Fig. 5C3), which highlighted changes in cellular activity throughout the network.

The dynamic memory map indicates that the activity in all cells was altered, regardless of whether a given cell was a direct target of modulation. Four main changes can be observed. First, the firing rate of many neurons increased during late protraction, near the transition from protraction to retraction. Second, the firing rate of several cells decreased during the early protraction phase. Direct comparison of the activity maps (Fig. 5C1 vs. C2) reveals that this change was due to a shortening of the activity in multiple cells. Third, the retraction phase primarily exhibited increases in activity, which were likely mediated by the recruitment of B51. Among other effects, this recruitment led to increased activity of B8 during retraction, thereby biasing patterns toward iBMPs. Fourth, the duration of activities increased during retraction, potentially due to a delay in termination of retraction (Fig. 5C3, reduced activity in the pattern terminator B52 at the 6–7 sec window). These results indicate that implementation of a memory engram led to population-wide reconfiguration, which comprised both direct and indirect effects of learning.

Intrinsic and synaptic plasticity make differential yet dependent contributions to the engram

One of the advantages of modeling studies is the ability to selectively incorporate individual sites of plasticity or various combinations of plasticity and then assess their relative contributions to changes in network activity. Here, comparisons were made in two steps. First, simulations characterized the relative contributions of intrinsic and synaptic plasticity to rhythmic activity. This analysis focused on pattern-initiating neurons B30, B63, and B65. Second, the individual contributions of each site of plasticity were evaluated.

Characterization of the relationship between modeled conductances and passive properties of B30, B63, and B65

Due to the presence of electrical coupling among B30, B63, and B65, changes in modeled conductances in one cell may affect the passive properties of the other two cells. Therefore, independent manipulation of each passive property (namely, input conductance and coupling ratio) required balanced changes to modeled conductances (namely, leak conductance and coupling conductance) in all three cells. We characterized the relationship among these variables by solving the equations for the input conductances and coupling ratios among B30, B63, and B65 over a range of parameters of leak conductances and coupling conductances in a three-neuron circuit (Fig. 6A; see Materials and Methods for details).

Figure 6.

Relationships among passive properties and modeled conductances. Due to the presence of electrical coupling among B30, B63, and B65, changes in one modeled conductance may affect the input conductance and coupling ratio of multiple cells. These effects were accounted for by solving the equations for the input conductances and coupling ratios among B30, B63, and B65 over a range of parameters in a three-neuron circuit (Eqs. 2–6). (A) Diagram of the three-neuron circuit and its component conductances, including leakage conductances (gleak) and coupling conductances (ges). An additional conductance (gothers) accounted for voltage-dependent conductances and coupling to neurons omitted from the three-neuron circuit (see Materials and Methods for details). (B) Effects of changes in modeled conductances on the B63–B30 coupling ratio. Although heatmaps display each variable as a function of relative changes in its two most impactful parameters, analytical solutions included all parameters (Eq. 2–6). All changes are relative to control, with “0” indicating no change and “1” a 100% change. Black circle and triangle mark respective positions of the control and operant parameters in parameter space. Black lines denote parameter values that yield no change in the plotted variable. (C) B63–B65 coupling ratio. (D) B30 input conductance. (E) B63 input conductance. (F) B65 input conductance.

Figure 6 shows the effects on each passive property of changes to the two main modeled conductances affecting it. For reference, the relative positions of control and operant memory values in the parameter space are also indicated. The slope of a color band indicates the relative contribution of each parameter for a given range of changes in the plotted passive property. Input conductances were chiefly, but not exclusively, determined by leakage conductance (Fig. 6D–F). Conversely, coupling ratios were affected similarly by coupling conductance and leakage conductance (Fig. 6B,C). Importantly, however, relative contributions were not constant throughout the parameter space, as indicated by the changing slopes of color bands, especially in the case of coupling ratios. These results confirmed that the effects of all parameters must be accounted for when targeting specific sets of changes to passive properties. Moreover, the characterization of these relationships made it possible for the parameter space to be searched for sets of values that achieve a given target set of changes to the passive properties. This approach was used in all subsequent analyses involving B30, B63, or B65 (see Materials and Methods for details).

Input conductance and electrical coupling make unique contributions to the engram

Next, the contributions of intrinsic and synaptic plasticity were compared by examining whether any particular feature of operant memory could be attributed to either form of plasticity in pattern-initiating neurons B30, B63, and B65. Two sets of simulations were run. The first set used parameters that modified the input conductances of B30, B63, and B65 to the same extent as the operant memory, but without affecting electrical coupling ratios. The second set introduced operant memory changes to coupling ratios without affecting input conductances. Neither set included changes to B51. Relative to control, changes to input conductances alone led to an increase in the overall rate of patterns (3.96 ± 0.02 BMPs/min), but also led to a decrease in the rate of iBMPs (0.10 ± 0.03 BMPs/min; Fig. 7A). Conversely, changes to coupling ratios alone resulted in a smaller increase in overall rate (2.30 ± 0.05 BMPs/min), combined with an increase in the rate of iBMPs (0.87 ± 0.09 BMPs/min; Fig. 7A). Both conditions led to reductions in the CV (0.038 ± 0.002 for changes to input conductances alone, and 0.065 ± 0.005 for changes to coupling ratios alone; Fig. 7B). Thus, each form of plasticity made unique contributions to the ensemble of features that characterize operant memory. Decreasing the input conductance of B30, B63, and B65 led to large increases in both the rate and regularity of motor patterns, but did so at the cost of a substantial reduction in bias toward ingestion-like motor patterns. Increases in coupling among these cells, on the other hand, seemed to contribute to all three characteristics of operant memory. However, coupling appeared to contribute more to regularity than it did to rate or bias, given that regularity approached operant values in the coupling-alone condition, unlike BMP rate or bias.

Figure 7.

Relative contributions of modified input conductance and electrical coupling to changes in network activity. (A) Decreasing only the input conductances of B30, B63, and B65 led to a large increase in motor pattern rate and reduced bias toward iBMPs. Conversely, increasing only the coupling ratio among these cells led to a smaller increase in rate, and to an increased bias toward iBMPs (n = 20 for each group). (B) Both changes in input conductance (Gin) and coupling ratio (CR) led to an increase in regularity (reduction in CV; n = 15 for each group). (C) Control-subtracted activity maps for changes in input conductances (C1) or coupling ratios (C2). Both sets of changes led to distinct activity reconfiguration, neither of which fully recapitulated the reconfiguration induced by operant memory (see Fig. 5D). (D) Motor pattern rate over a range of changes to input conductances and coupling ratios among B30, B63, and B65. Both variables contributed to BMP rate, although input conductance had a larger effect. Dashed black lines denote the effects of changing one variable while keeping the other constant at control values. White line denotes the diminished effect of varying coupling ratios when input conductances are reduced by 25%. Changes in coupling ratios and input conductances are relative to control, with “0” indicating no change and “1” a 100% change. A total of 95 unique sets of parameters were simulated. (E) Contribution of changes in input conductances or coupling ratios to motor pattern regularity. Both variables affected regularity. The various peaks and valleys indicate that unique combinations of changes may have unique effects on regularity. White line denotes the diminished effect of varying input conductances when coupling ratios were increased by 75%. Surface plots on panels D and E display median values for simulations ran over a subset of the parameter space. Note that CVs were only calculated for a further reduced subset of the parameter space (60 out of 95 unique parameter sets) in which simulations generated at least 10 BMPs over 6 min.

Contributions to memory are mediated by unique reconfigurations of network activity

Can the contributions of each form of plasticity be attributed to specific changes in cellular activities? Dynamic memory maps (Fig. 7C) revealed that, similarly to the complete operant memory, both manipulations reconfigured activity throughout the network. BMPs generated under each condition recapitulated a subset of the changes induced by operant memory (Fig. 5C3). Specifically, implementation of modifications to input conductance alone (Fig. 7C1) recapitulated both the increases in firing rate of several neurons during late protraction and the decreases during early protraction (e.g., compare B8, B20, B34, B63, and B65 in Fig. 7C1 and Fig. 5C3), but did not recapitulate the increases in spike rate during retraction or the lengthening of that phase (e.g., compare B52 and B64 in Fig. 7C1 and Fig. 5C3). In contrast, the coupling-alone condition (Fig. 7C2) recapitulated only the decreases in the activities of several cells during early protraction. Given that both conditions led to substantial regularization of rhythmic activity and activity decreases during early protraction (Fig. 7B,C), it seems likely that shortening the windows during which protraction neurons may be active played a role in controlling BMP regularity. By exclusion, then, the increases observed during late protraction likely contributed to the higher overall BMP rate induced by decreasing input conductance.

Effects of intrinsic and synaptic plasticity are mutually dependent

Experimental investigations suggest that the input conductances of B30, B63, and B65 contribute to motor pattern rate but not regularity, and that coupling contributes to regularity but not pattern rate (Sieling et al. 2014), in contrast to the more interrelated effects obtained with the model. To explore this difference, we examined to what extent the effects of one form of plasticity depend on the other over a wide parameter space. The characterized relationships between passive properties and modeled conductances in these three cells (Fig. 6) were used to obtain sets of parameters that yielded a range of combinations of changes to input conductance and coupling ratio (see Materials and Methods for details). Simulations were run for 95 unique sets of parameters, and the rate and CV of patterned activity were measured. The input conductances of all three cells were modified simultaneously and to the same extent, as were the B63–B30 and B63–B65 coupling ratios. Figure 7D displays the overall rate of motor pattern generation for each combination of changes. As indicated by the black dashed lines on top of the surface plot, manipulating either passive property while keeping the other at control values altered the BMP rate. However, the effect of changes in coupling ratio depended on the input conductance—as the latter increased or decreased, the effects of the former on pattern rate were diminished (e.g., white dashed line). The effect of input conductance on rate also diminished as the coupling ratio increased, but to a lesser extent. The regularity of rhythmic activity could only be assessed for the 60 (out of 95) unique sets of parameters that yielded simulations that produced at least 10 BMPs in 6 min (Fig. 7E). As was the case with BMP rate, changes in either passive property were sufficient to modify the CV when the other property was maintained at control values, but the effects of both input conductance and coupling ratio on regularity were diminished in specific regions. Moreover, whereas the BMP rate plot was smooth, the CV plot displayed various peaks and valleys, which may indicate a more complex relationship in which unique combinations of changes have unique effects on regularity. These results demonstrated that there can be, in principle, regions in parameter space in which rate and regularity are primarily affected by one passive property but not the other.

Individual sites of plasticity uniquely contribute to features of rhythmic activity

We next examined the contributions made by individual sites of plasticity in isolation. This analysis was performed by simulating six additional conditions, each including only changes either to the input conductance of one cell or to the coupling ratio between two cells. The effects of each condition on motor pattern rate, bias, and regularity are displayed in Figure 8. Compared to control values, implementing changes to the input conductance of B51 increased the rate of iBMPs (1.40 ± 0.07 BMPs/min) while slightly decreasing the overall rate of patterns (1.49 ± 0.07 BMPs/min) and leading to effectively no change in regularity (CV = 0.163 ± 0.013). Unexpectedly, modifying the input conductances of B30, B63, or B65 individually had strikingly distinct effects. The B30 change led to increased regularity (CV = 0.079 ± 0.005), a small increase in iBMP rate (0.59 ± 0.06 BMPs/min), and no effective change in overall rate (1.77 ± 0.05 BMPs/min). Modifying B63, on the other hand, increased the overall rate of BMPs (3.35 ± 0.13 BMPs/min), decreased the rate of iBMPs (0.12 ± 0.03 BMPs/min), and caused no clear change in regularity (CV = 0.192 ± 0.053). Implementing changes to the input conductance of B65 led to decreases in overall rate (1.39 ± 0.07 BMPs/min), iBMP rate (0.10 ± 0.03 BMPs/min), and regularity (CV = 0.356 ± 0.065). Thus, changes to B65 had effects that were opposite to those of the complete engram on all three features of operant memory. Finally, altering the B63–B30 and B63–B65 coupling ratios led to similar effects on two out of three features. For B63–B30 coupling, there was an increase in regularity (CV = 0.096 ± 0.012) and a small increase in rate of iBMPs (0.63 ± 0.05 BMPs/min), but the overall rate of patterns (1.80 ± 0.08 BMPs/min) was unchanged. These changes were similar to those observed when manipulating the input conductance of B30 alone. Conversely, for B63–B65 coupling, rate of iBMPs (0.68 ± 0.06 BMPs/min), regularity (CV = 0.086 ± 0.009), and overall rate of patterns (2.26 ± 0.02 BMPs/min) were all increased. These results indicated that, like the combined changes examined in the previous subsection, each specific site of plasticity made mostly unique contributions to the total set of features of operant memory. Furthermore, although most individual effects were in the same direction as the effects of the complete engram, effects in the opposite direction on at least one of the behavioral features of memory were common. Last, in at least one case there was evidence of synergism among sites of plasticity (i.e., the contributions of combined sites of plasticity could exceed the sum of the individual contributions of each site). Although B30 had the only input conductance change that increased regularity (Fig. 8B, CV = 0.079 ± 0.005), implementing changes to the input conductances of B30, B63, and B65 in combination led to a greater increase in regularity than B30 alone (Fig. 7B, CV = 0.038 ± 0.002).

Figure 8.

The relative contributions of individual sites of plasticity in isolation. (A) Contribution of each site of plasticity to pattern rate and bias (n = 20 for each group). Notably, changes in the input conductance of B51 contributed a strong bias toward iBMPs at the cost of a small reduction in overall rate, whereas changes in the input conductance of B63 contributed a large increase in pattern rate at the cost of a decrease in bias toward iBMPs. (B) Contribution of each site to BMP regularity (n = 15 for each group). The input conductance of B30 and each of the two coupling ratios all contributed to a decrease in CV. Simulations used to calculate the CV (but not BMP rate) for B51 and B65 input conductances were extended from 6 to 10 min so as to obtain 10 BMPs for CV calculation in each simulation.

These findings suggest that: (1) each site of plasticity in the engram contributed only a portion of the overall changes in network output; (2) part of the effects of individual sites could be opposite in direction to the effects of memory; (3) combined sites of plasticity could modulate each other's contributions and display synergism; and (4) memory emerged from the combined effects of all sites of plasticity.

Taken together, these investigations suggested that intrinsic and synaptic plasticity made differential contributions to operant memory through differential reconfigurations of network activity. Notably, these contributions appeared to consist of unique combinations of changes to multiple behavioral features simultaneously, each of which could be in the same or opposite direction as those induced by the complete engram. Moreover, the contributions of intrinsic and synaptic plasticity seemed to be mutually dependent—that is, one form of plasticity could modulate the contributions of the other.

Discussion

Comparison with previous models of the feeding CPG

Several previous studies described models of the Aplysia feeding neural network that can simulate aspects of feeding. Some of these previous models were abstract in design (Kupfermann et al. 1992; Deodhar and Kupfermann 2000), whereas others adhered to biological constraints (Ziv et al. 1994; Susswein et al. 2002; Cataldo et al. 2006; Hurwitz et al. 2008). The present model was an extension of our previous studies of the Aplysia feeding circuit (Kabotyanski et al. 1994; Ziv et al. 1994; Susswein et al. 2002; Cataldo et al. 2006, 2009; Baxter et al. 2010), and represents the most comprehensive model to date. Several aspects of the model, however, are still incomplete. For example, the present model did not include all of the neurons known to be elements of the CPG nor the possibility of additional sites of plasticity in identified neurons and yet unidentified neurons. Moreover, empirical studies have yet to provide complete details of the biophysical properties of many neurons and synaptic connections. Nevertheless, simulations indicated that the present model network captured salient features of rhythmic activity in the feeding neural network. For example, the subset of neurons in the model could simulate the genesis of BMPs and the switch between different types of patterns.

Roles of individual and combined sites of plasticity

Simulations in this study suggested that sites of intrinsic and synaptic plasticity can make unique contributions to memory. These contributions were mediated by differential reconfiguration of network activity, and consisted of unique sets of concurrent changes to multiple behavioral features. In addition, each change could be in the same or in the opposite direction of the changes induced by the complete engram. Furthermore, changes induced by intrinsic and synaptic plasticity were found to be mutually dependent. These findings support a view of memory expression as a distributed process that emerges from a complex set of reciprocally modulating sites of intrinsic and synaptic plasticity.

The model indicated that the currently identified sites of one type of learning-induced plasticity could, at least in principle, account for observed changes in the behavioral features of feeding examined here. However, this finding does not imply that the known sites constitute the complete engram for several reasons. First, it is possible for multiple sites of plasticity to contribute to the same behavioral effect. Indeed, in the model, each of the three examined behavioral features received contributions, ranging from large to small, from several individual sites of plasticity. Second, multiple loci may have redundant roles, which would contribute robustness to the memory trace. Decreasing the input conductance of B30 or increasing the B63–B30 coupling yielded similar changes to all three features of network activity in the model, suggesting that redundant components could be present in the engram. Third, not all sites recruited as part of an engram need to contribute to the memory. Findings from the model suggested that each site of plasticity led to a unique set of changes in BMP features, often including effects that were in the opposite direction of changes induced by the complete engram. Moreover, one site, the decrease in input conductance of B65, appeared to lead solely to such opposite effects when examined in isolation. It is possible that B65 could contribute to other features of memory that were not examined.

In operant learning behavioral features can be selected regardless of whether they are useful. Indeed, the training paradigm used in operant reward learning in Aplysia feeding merely makes reinforcement contingent upon ingestion responses. Thus, the highest number of rewards would be obtained by a high rate of ingestions, regardless of how regular are the intervals between feeding responses. Nevertheless, feeding responses become highly regular following training (Nargeot et al. 2007). Why, then, is regularity selected at the behavioral level? One possibility is that regularity is enhanced due to some unknown behavioral relationship. Alternatively, it could be a by-product of memory allocation at the neuronal level.

Results from the model indicate that individual sites of plasticity make unique sets of contributions to multiple behavioral features. Moreover, further contributions can arise from the interactions among multiple loci, such as synergistic enhancement of a feature. Thus, the engram appears to be a combinatoric code of plasticity loci in which addition or removal of a single locus can have effects that exceed its own direct contributions. One potential mechanism of memory allocation involves the occurrence of activity in a neuron at a given critical time (Athalye et al. 2018; Josselyn and Frankland 2018). This mechanism requires some form of coincidence detection between activity and a critical event. Such a mechanism has been shown to mediate the changes induced by operant learning in B51, which require activity-dependent Ca2+ influx into the cell and concurrent activation of dopamine receptors (Lorenzetti et al. 2008). A similar motif in other cells may explain why neurons such as B65, which does not appear to make any direct contribution to memory, are recruited into the engram. Neurons B30, B63 and B65 display some level of activity in most motor patterns. Thus, their recruitment into the engram may be automatic following operant training. If that is the case, their synergistic effect on regularity would explain why this behavioral feature emerges after training. Several predictions can be made from this interpretation. First, a training paradigm designed to exclusively reward irregular feeding responses should fail to enhance irregularity and still lead to highly regular ingestions. Second, hyperpolarizing or buffering Ca2+ in B30 and B65 during training should prevent their recruitment into the engram and lead to a high rate of irregular BMPs. Third, several other neurons in the CPG that are reliably active during motor patterns should also be recruited as part of the engram.

Roles of intrinsic and synaptic plasticity

Plasticity of the intrinsic properties of neurons is at times portrayed as secondary in the more dominant view of synaptic plasticity as the primary site of memory, in which case the roles of intrinsic plasticity in modulating synaptic plasticity are emphasized. However, evidence from invertebrates and vertebrates supports intrinsic plasticity as a site of memory in its own right (Hansel et al. 2001; Mozzachiodi and Byrne 2010; Titley et al. 2017). The findings reported herein suggest not only that intrinsic and synaptic plasticity make differential and complementary contributions to memory, but also that these forms of plasticity may reciprocally modulate each other's contributions. Thus, expression of memory may require combinatoric engagement of intrinsic and synaptic plasticity, without which certain effects may not be achievable.

A recent study empirically evaluated the contributions of changes to the input conductance of B30, B63, and B65 or their electrical coupling by artificially adding currents that mimic these changes to cells from a naïve preparation (Sieling et al. 2014). Consistent with the results reported here, decreases to input conductance led to an increase in the rate of motor patterns, and increases in coupling led to an increase in BMP regularity. However, the study found that decreases in input conductance had no effect on the regularity with which motor patterns occurred, in contrast to the effects of the equivalent manipulation in the computational model. The analysis in Figure 7E shows that there are regions in parameter space in which the effects of input conductance on regularity are greatly diminished. One possibility, then, is that the biological system occupies such a region, which the model did not reproduce under control parameters. A second possibility is that this difference may be explained by the fact that the current iteration of the model does not fully capture irregularity in the system. Variability in which neuron among B30, B63, and B65 has the earliest onset of activity for a given motor pattern is present in naïve animals and is reduced by operant conditioning (Nargeot et al. 2009). Although these three neurons are capable of initiating BMPs in the model, all simulations in the present study had motor patterns evoked by the same tonic input to CBI-2, and B63 was consistently the first activated neuron among the three cells in both control and operant memory simulations. Thus, expansion of the model to better reflect sources of irregularity in the system could improve the ability to simulate learning-induced changes throughout the neuronal population.

In intact, freely behaving animals, the absolute rate of feeding is a function of several variables, including motivational state (e.g., food deprivation), chemical and tactile stimulation. Some of these variables converge on neurons such as CBI-2 and other CBIs, and neurons in the buccal ganglia. Indeed, Rosen et al. (1991) show that both rate and regularity of motor patterns can be modified by varying the intensity of CBI-2 stimulus. This is also true of the model—higher magnitude of CBI-2 current injection leads to higher frequency and regularity (data not shown). The overall rate observed in our simulations under control conditions approximated that observed in naïve isolated ganglia in operant conditioning experiments (e.g., Sieling et al. 2014). One limitation of this study is that all patterns were elicited by tonic stimulation of CBI-2 with a fixed intensity. In the biological system, motor patterns can be initiated by many other neurons, including other CBIs and pattern initiating neurons in the buccal ganglia (e.g., B31, B63). This study did not explore whether the same relationships among sites of plasticity would be maintained when patterns are initiated by other means. Expansion of the model to include other CBIs and pattern initiating neurons will allow for similar questions to be addressed.

Relationship with biomechanical models

In addition to neuronal models describing the CPG, such as the one reported herein, biomechanical models of the Aplysia feeding apparatus have also been developed (Drushel et al. 1998, 2002; Lyttle et al. 2017). These models can describe the positions and forces exerted by the muscles in the animal's buccal mass, thereby using motor output to illuminate function (e.g., Novakovic et al. 2006). The conductance-based CPG model described in the present study could be combined with such a biomechanical model. This would make it possible to quantify relationships between cellular and synaptic properties of the network and ultimate motor output and function. Examining the contributions of individual sites of plasticity in this combined model would allow for a more detailed perspective of the behavioral role of each plasticity locus. Moreover, if implementation of the currently known engram proves insufficient to produce all motor changes observed experimentally, novel potential sites of plasticity could be predicted by the model.

Limitations of the model

The present model was designed in part to determine whether it is in principle possible for the currently known sites of operant memory to reproduce the main features of operant learning. Although the answer is affirmative, whether additional sites are necessary in vivo, or with different parameter sets, is still unknown. Furthermore, how robust these changes are in the face of potential interferences by other memories or perturbations remains unclear. Moreover, there is no current experimental evidence on whether there is redundancy in sites of plasticity playing similar roles, as can be expected from a robust system. Whereas many other pattern-initiating and pattern-selecting neurons would be in a position to play a role in operant memory, only a few neurons have been examined thus far. Indeed, dopamine, which mediates operant conditioning (Nargeot et al. 1999c), can modulate the activity of many neurons throughout the network (Neveu et al. 2017). Therefore, as empirical studies identify new sites of plasticity, the present model can be expanded to include these new sites and test their specific contributions to network function and reconfiguration.

The model did not include descriptions of modulatory processes or the molecular processes that underlie the induction of operant learning-related changes in cellular properties. Empirically, the learning-induced sites of plasticity are characterized post training, and thus represent the expression of the memory engram. Currently, few data are available regarding the underlying molecular mechanisms and dynamics of inducing memory during operant learning (Lorenzetti et al. 2008; Mozzachiodi et al. 2008). Thus, the present model mimicked the established sites of plasticity that are known to be present following learning (see Fig. 4). A future goal is to analyze induction processes and incorporate them into the model.

Materials and Methods

Developing the model

The feeding neural network has been analyzed extensively and is known to consist of at least 50 identified neurons including sensory neurons, interneurons, and motor neurons in the buccal and cerebral ganglia (Susswein and Byrne 1988; Church and Lloyd 1994; Cropper et al. 2004, 2019). Additional neurons likely are as yet unidentified. Our goal was to develop a reduced model of the CPG component of the network that was sufficient to capture the salient features of the genesis of BMPs for ingestion and rejection of food. The reduced model included 13 cells: CBI-2, B4, B8, B20, B30, B31, B34, B40, B51, B52, B63, B64, and B65 (Fig. 1). CBI-2 is a command-like neuron that elicits fictive feeding (Rosen et al. 1991; Hurwitz et al. 2005). B8 is a closure motor neuron (Morton and Chiel 1993b). The remaining cells are elements of the CPG (Susswein and Byrne 1988; Plummer and Kirk 1990; Hurwitz and Susswein 1996; Hurwitz et al. 1997, 2003, 2005, 2008; Kabotyanski et al. 1998; Jing and Weiss 2001, 2002, 2005; Susswein et al. 2002; Jing et al. 2003, 2004; Dembrow et al. 2004; Proekt et al. 2007; Wu et al. 2007, 2010; Sasaki et al. 2008, 2009, 2013; Saada et al. 2009; Saada-Madar et al. 2012; Dacks and Weiss 2013). Previous modeling studies predicted that additional excitatory drive onto B64, mediated by a hypothetical “z cell,” was required for patterns to transition from the protraction to the retraction phase (Cataldo et al. 2006). Although at least one promising candidate for this role has now been described (i.e., the complex multifunctional neuron CBI-5; Sasaki et al. 2007, 2008), modeling studies have yet to explore its necessity and sufficiency for phase transitions. Therefore, a hypothetical excitatory synaptic connection from B63 to B64 was included in our reduced model as a parsimonious alternative.

The model was developed using SNNAP (Simulator for Neural Network and Action Potentials, version 8.1). The operation and capabilities of SNNAP are described elsewhere (Ziv et al. 1994; Av-Ron et al. 2006, 2008; Baxter and Byrne 2007). Specific details of the equations and parameters that were used in the present study are included in the SNNAP input files, which are available at the ModelDB website (http://modeldb.yale.edu/261489). These input files are annotated ASCII text files, and can be viewed using a text editor.

Empirical data provide some biological constraints for the reduced model such as the pattern of monosynaptic connections and electrical coupling among the cells (the connectome), homosynaptic plasticity, the firing properties of individual cells in response to stimuli, cellular patterns of activity during BMPs, and sites of learning-induced plasticity. In general, each cell was represented by a single compartment. However, cells B31, B51, and B64, which express plateau-like potentials, were exceptions. These three cells were modeled with two compartments: a somatic compartment and an axonal compartment. The somatic compartment included slow conductances that mediated plateau potentials. The axonal compartment included fast conductances that mediated spiking. This strategy of separating the slow and fast conductances followed an approach described by Vavoulis et al. (2007).

Because the neurons constituting the reduced model are a subset of the complete CPG circuit, we needed to adjust the parameters of individual neurons and synaptic connections to obtain realistic firing properties and BMP generation using physiologically reasonable, but not necessarily exact, values. Parameters were adjusted by trial and error. Briefly, in each of the 13 cells, action potentials were mediated via a fast Na+ conductance and a delayed K+ conductance. The spike threshold was adjusted by either hyper- or depolarizing shifts in the inflection point of the Boltzman functions that defined the voltage-dependence of the membrane conductances (see SNNAP input files for details). To reproduce the unique spiking properties of specific neurons, additional conductances, such as an A-type K+ conductance or an H-type current, were included as needed. Similarly, the properties of synaptic conductances were designed to match the unique attributes of each monosynaptic connection. Thus, the conductance-based model approximated the unique firing properties of identified cells (e.g., tonic spiking, delayed spiking in response to stimuli, frequency adaptation, plateau potentials, rebound excitation), the pattern of electrical coupling among cells, and the properties of each monosynaptic connection (e.g., homosynaptic plasticity, fast and slow postsynaptic potentials, increasing versus decreasing synaptic conductances, voltage-dependent synaptic conductances). Parameters were further adjusted so that the simulated network generated ingestion and rejection BMPs that resembled those recorded physiologically. Given that the model could produce biologically relevant patterns of activity (see Fig. 2), we considered it sufficient to investigate the general issue of how multiple sites of learning-induced plasticity contribute to changes in network activity.

Because of inherent variability in synaptic and membrane conductions, models whose function relies heavily on precise values of parameters may be unrealistic (Prinz et al. 2004; Grashow et al. 2010; Marder 2011; Gutierrez et al. 2013). As an attempt to address this issue, stochastic fluctuations were included in all simulations. In SNNAP, any conductance, including membrane, synaptic and electrical-coupling conductances, can include noise. For example, a general equation for a membrane current is where is the maximum conductance, A and B are Boltzman functions for activation and inactivation, respectively, V is membrane potential, E is the equilibrium potential and is the product of various modulator functions (f[REG]). Noise is introduced by R, which is a random number from a Gaussian distribution of mean zero and SD proportional to . Equations for synaptic current and current due to electrical coupling also included the R variable. Noiseless simulations were only used to obtain exact input values for certain analytical solutions (see last subsection of Materials and Methods), all other simulations included equal amounts of noise. All conductances in the model included noise with an SD equal to 25%, and R was updated every 2 msec. By including noise, the model was shown to be robust to fluctuations in parameter values.

Classifying simulated BMPs as either fictive ingestion (iBMPs) or other patterns

To classify simulated BMPs as being either iBMPs or other BMP types, the present study followed criteria that were established previously (Nargeot et al. 1997). The protraction phase was defined as spike activity in the axonal compartment of B31 (B31a), and the beginning of the retraction phase was defined as the termination of spike activity in B31a (see Fig. 3). Spike activity in B8 defined closure. A simulated BMP was classified as being an iBMP if ≥50% of the activity in B8 occurred after the termination of spike activity in B31a (i.e., during the retraction phase). BMPs that did not meet this criterion were classified as “others,” which included (1) intermediate patterns, in which B8 activity partially overlapped with both the protraction and retraction phases, but the majority of activity in B8 occurred during the protraction phase; (2) fictive rejection (rBMP), in which activity in B8 and B31a co-terminated; and (3) incomplete patterns, in which one of the three phases was lacking (e.g., B8 was not active at any point).

Measuring behavioral features of operant reward memory

Empirical studies describe at least three outcome measures of operant learning. First, learning increases the frequency of biting in vivo, and increases the frequency of BMPs in vitro (Brembs et al. 2002; Nargeot et al. 2007, 2009; see also Nargeot et al. 1997; Mozzachiodi et al. 2008). Second, learning biases activity in the CPG in vitro toward the expression of iBMPs (Nargeot et al. 1997, 1999b; Brembs et al. 2002; Mozzachiodi et al. 2008). Third, learning increases the regularity of biting in vivo and of BMPs in vitro (Nargeot et al. 2007, 2009). These outcome measures were used in the present study.

Characterizing the regularity of simulated rhythmic activity

Empirical studies indicate that in naïve animals, biting behavior occurs sporadically (Brembs et al. 2002; Nargeot et al. 2007, 2009). Similarly, in naïve isolated buccal ganglia preparations, the expression of BMPs is sporadic (Nargeot et al. 1997, 2007, 2009). One of the outcomes of learning is the regularization of biting in vivo and BMPs in vitro (Nargeot et al. 2007, 2009). In the present study, regularization of simulated BMPs was characterized by calculating the coefficient of variation (CV) of the inter-burst interval, where CV = SD/mean (see Figs. 3, 5). This method has been used previously to characterize the regularity of behaviors and spiking and/or bursting activities (e.g., Büschges et al. 2000; Moortgat et al. 2000a,b; Prut and Perlmutter 2003; Johnston et al. 2013).

Creating activity maps of simulated rhythmic activity

Activity maps were used to characterize the level and timing of spike activity for cells in the network. Activity maps were generated for a period of time ±10 sec relative to the transition from protraction to retraction in a BMP (see Fig. 3). This 20 sec epoch was divided into 20, 1 sec, nonoverlapping bins. The number of spikes during each bin was counted in each cell and averaged over 100 BMPs. This measure represented the mean firing frequency (Hz) of each cell during a BMP. These activity maps were used to analyze differences in activity in control simulations versus simulations that included learning-induced plasticity (see Figs. 5, 7). All analyses of simulated activity were performed by custom routines in MATLAB.

Implementing learning-induced plasticity

As was the case with all conductances (see above), conductances implicated in learning were not constrained by specific values reported in the literature (either before or after learning). Instead, modifications to specific conductances were made by trial and error until simulations reproduced the main features of operant reward learning, namely changes in motor pattern rate, regularity, and bias toward iBMPs (Nargeot et al. 1997, 1999a, 2007, 2009; Brembs et al. 2002; Mozzachiodi et al. 2008; Sieling et al. 2014). Previous empirical studies characterized at least three sets of correlates of operant reward learning: (i) a decrease in the input conductances and increase in the excitability of cells B30, B63, and B65 (Nargeot et al. 2009; Sieling et al. 2014); (ii) a decrease in the input conductance and increase in the excitability of cell B51 (Nargeot et al. 1999b; Brembs et al. 2002; Lorenzetti et al. 2008; Mozzachiodi et al. 2008); (iii) an increase in the B63–B30 and B63–B65 electrical coupling coefficients and conductances (Nargeot et al. 2009; Sieling et al. 2014). In the present study, the learning-induced changes in the input conductance and excitability of B30, B63, and B65 were simulated by decreases in leakage conductance (see Fig. 4A,C). The leakage conductances were decreased from control values of 100 nS in B30, 40 nS in B63, and 80 nS in B65 to 80, 20, and 60 nS, respectively. Learning-induced changes in B51 were likewise simulated by decreasing the leakage conductance of the cell from control values of 100 nS in the somatic compartment (B51s) and 80 nS in the axonal compartment (B51a) to 80 nS in B51s and 64 nS in B51a (see Fig. 4D). Decreasing the leakage conductance was a parsimonious method that proved sufficient for decreasing input conductance and increasing excitability in all four cells. Finally, learning-induced changes in electrical coupling were simulated by increasing the coupling conductances between B63 and B30 from 6 to 7 nS, and between B63 and B65 from 6 to 8 nS (see Fig. 4B).

Implementing individual sites of plasticity in isolation and combination

Modifying the leakage conductance of B51 did not lead to changes in any other known sites of plasticity. Thus, the excitability of B51 could readily be manipulated in isolation. Conversely, because B30, B63, and B65 are electrically coupled, changing the leakage or coupling conductances in any of these cells can affect multiple sites of plasticity. Therefore, manipulation of any single site or set of sites of plasticity without affecting others required searching the parameter space of leakage and coupling conductances for a combination of parameters that met this goal. This was achieved by solving the equations for all input conductances and coupling ratios in an idealized circuit consisting of only the three cells with their respective coupling conductances and leakage conductances (Fig. 6A). This approach allowed for all input conductances and coupling ratios to be characterized as a function of all leakage conductances and coupling conductances (Fig. 6B–F). The three-neuron circuit was made to approximate the passive properties of B30, B63, and B65 when they are embedded in the complete network by including an additional conductance in each cell, termed g. This conductance represented the net conductance resulting from voltage-dependent conductances and electrical coupling to all other neurons omitted in the three-neuron circuit (Fig. 6A). Thus, the input conductances (G) and coupling ratios (CR) in the three-neuron circuit were given by the following equations: where g represents the leakage conductance in one cell and g represents the coupling conductance among two cells. The value of g for each neuron was estimated by performing a series of simulations of the complete network in which all three neurons were voltage clamped to the same values over a range from −100 mV to 40 mV. Because B30, B63, and B65 were clamped to the same voltage, there was no current due to coupling among these cells. Current due to the intrinsic leakage conductance was then subtracted from the total current in each cell, yielding a net current resulting exclusively from conductances that were to be omitted in the three-neuron circuit. Finally, a linear fit of the current–voltage relationship in hyperpolarized values was used to obtain g for each neuron. Simulations used to estimate g did not include noise, so as to obtain exact values. Input conductances and coupling ratios resulting from analytical solutions using these values were systematically compared with direct measurements from hyperpolarizing current injection simulations. Observed discrepancies were minimal, and these discrepancies were further shown to be due to small contributions by nonlinear voltage-dependent conductances.

To obtain specific target combinations of coupling ratios and input conductances, Eqs. 2–6 were solved iteratively. Initially, all parameters were set to control values (i.e., g = 100, g = 6, g = 40, g = 80, g = 6, all values are in nS). Each iteration consisted of three steps. First, Eqs. 2 and 3 were solved for g and g. Second, these values were used to solve Eq. 4 for g. Third, values obtained in the previous steps were used to solve Eqs. 5 and 6 for g and g. This process was repeated until input conductances and coupling ratios converged on the target values. This occurred in 6–8 iterations in most cases, and deviations were generally less than 2.5 × 10−12% from the target. Estimation of g and solutions to all equations were computed by custom routines in MATLAB.