## Abstract

Recent work in computational neuroethology has emphasized that “the brain has a body”: successful adaptive behavior is not simply commanded by the nervous system, but emerges from interactions of nervous system, body, and environment. Here we continue our study of these issues in the accessory radula closer (ARC) neuromuscular system of *Aplysia.* The ARC muscle participates in the animal's feeding behaviors, a set of cyclical, rhythmic behaviors driven by a central pattern generator (CPG). Patterned firing of the ARC muscle's two motor neurons, B15 and B16, releases not only ACh to elicit the muscle's contractions but also peptide neuromodulators that then shape the contractions through a complex network of actions on the muscle. These actions are dynamically complex: some are fast, but some are slow, so that they are temporally uncoupled from the motor neuron firing pattern in the current cycle. Under these circumstances, how can the nervous system, through just the narrow channel of the firing patterns of the motor neurons, control the contractions, movements, and behavior in the periphery? In two earlier papers, we developed a realistic mathematical model of the B15/B16-ARC neuromuscular system and its modulation. Here we use this model to study the functional performance of the system in a realistic behavioral task. We run the model with two kinds of inputs: a simple set of regular motor neuron firing patterns that allows us to examine the entire space of patterns, and the real firing patterns of B15 and B16 previously recorded in a 21/2-h-long meal of 749 cycles in an intact feeding animal. These real patterns are extremely irregular. Our main conclusions are the following. *1*) The modulation in the periphery is necessary for superior functional performance. *2*) The components of the modulatory network interact in nonlinear, context- and task-dependent combinations for best performance overall, although not necessarily in any particular cycle. *3*) Both the fast and the slow dynamics of the modulatory state make important contributions. *4*) The nervous system controls different components of the periphery to different degrees. To some extent the periphery operates semiautonomously. However, the structure of the peripheral modulatory network ensures robust performance under all circumstances, even with the irregular motor neuron firing patterns and even when the parameters of the functional task are randomly varied from cycle to cycle to simulate a variable feeding environment. In the variable environment, regular firing patterns, which are fine-tuned to one particular task, fail to provide robust performance. We propose that the CPG generates the irregular firing patterns, which nevertheless are guaranteed to give robust performance overall through the actions of the peripheral modulatory network, as part of a trial-and-error feeding strategy in a variable, uncertain environment.

## INTRODUCTION

The ultimate goal of nervous systems is to generate successful adaptive behavior, behavior that enhances the survival and reproduction of the animal. However, no disembodied nervous system can achieve this goal: “the brain has a body” (Chiel and Beer 1997). The motor commands of the nervous system must be such as to engage the properties of the animal's body, and those of its environment, so as to change the body, the environment, or their mutual relationship in the required adaptive way. Adaptive behavior is thus not simply “commanded” by the nervous system, but “emerges from interactions of nervous system, body, and environment” (Beer 1995; Beer et al. 1999; Chiel and Beer 1997).

The body, indeed, can have “a mind of its own” (Raibert and Hodgins 1993). In vertebrates as well as invertebrates, the structural and dynamical complexity of the periphery can be as large as that of the central nervous system, so that, seen more abstractly, the computational capability of the periphery rivals that of the nervous system that is attempting to control it. The center and periphery are connected through a channel, the neuromuscular transform (NMT) (Brezina et al. 2000b), which can be very narrow relative to the rich computational capabilities at either end (Brezina and Weiss 2000; Brezina et al. 2000c). The semiautonomous periphery may then facilitate and simplify, but can also considerably complicate and constrain, motor control (Beer 1990; Brezina and Weiss 2000; Chiel and Beer 1997; Dickinson et al. 2000). A fundamental question is how the two halves of the system are matched in their properties, and coordinate in their dynamics, so as to ensure the production of functional behavior.

Here we continue our long-standing study of this question in an advantageous model system, the accessory radula closer (ARC) neuromuscular system of *Aplysia.* The ARC muscle is controlled by just two motor neurons. In contrast to this central simplicity, there is great peripheral complexity, generated by a dynamically complex network of actions of numerous endogenous modulatory neurotransmitters and neuropeptides that shape the contractions of the muscle. For any central motor command, this modulatory network multiplies the range of possible peripheral outputs of the NMT. Extensive studies over the past 25 years (reviewed by Kupfermann et al. 1997; Weiss et al. 1993; see also Brezina et al. 2003a) have characterized many aspects of the ARC system, in particular the various modulatory actions in quantitative terms, to the point where integrative understanding has become possible. We now briefly summarize the details essential for the work in this paper.

The ARC muscle is one of the muscles of the buccal mass, a complex structure that produces rhythmic, cyclical movements of the animal's food-grasping organ, the radula, in several types of consummatory feeding behavior such as biting, swallowing, and rejection of unsuitable food (Kupfermann 1974). All of these behaviors are driven by a single, multitasking central pattern generator (CPG) (for recent review see Elliott and Susswein 2002) whose patterned activity is then distributed to the motor neurons of the various buccal muscles. The ARC muscle's two motor neurons, B15 and B16, release their classical transmitter, ACh, to depolarize and contract the muscle (Cohen et al. 1978; Kozak et al. 1996). In addition, however, both motor neurons release modulatory peptide cotransmitters that shape the basal ACh-induced contractions. B15 releases, in particular, the small cardioactive peptides (SCPs) (Cropper et al. 1987a, 1990b; Lloyd et al. 1984; Vilim et al. 1996b; Whim and Lloyd 1989), and B16 releases the myomodulins (MMs) (Brezina et al. 1995; Cropper et al. 1987b, 1991; Vilim et al. 2000). The SCPs and MMs shape contractions through three main actions on the muscle: *1*) they potentiate the contractions by enhancing the muscle's depolarization-activated Ca current that supplies Ca^{2+} essential for contraction (Brezina and Weiss 1995; Brezina et al. 1994a, 1995; Cropper et al. 1987a, b, 1991; Lloyd et al. 1984; Whim and Lloyd 1990); *2*) they depress the contractions by activating in the muscle a K current that opposes the ACh-induced depolarization, activation of the Ca current, and Ca^{2+} influx (Brezina et al. 1994b, 1995; Cropper et al. 1987b, 1991; Orekhova et al. 2003); and *3*) they increase the relaxation rate of the contractions probably by modulating the muscle's contractile machinery (Brezina et al. 1995; Cropper et al. 1990b; Heierhorst et al. 1995; Probst et al. 1994; Whim and Lloyd 1990). The balance of the competing potentiating and depressing actions determines the net modulation of contraction size, which, together with the modulation of the relaxation rate, then produces the final shape of the contraction. [It should be mentioned that all this happens against the background of yet other, although probably less critical, actions of these and other modulators that are present in or act on the ARC system (see Brezina et al. 2003b).] The main elements of this modulatory network are graphically summarized within the gray box in Fig. 1.

Our understanding of how the ARC system operates has over the years been embodied in a series of models, at first purely conceptual (e.g., Weiss et al. 1993) but increasingly, with more quantitative measurements made in the real system, more realistic models (Brezina and Weiss 2000; Brezina et al. 1996, 2000b, c). Most recently, in Parts I and II of the present study, we have constructed a fully quantitative dynamic model of the network of modulation in the ARC system and of its dependency on the firing of the two motor neurons B15 and B16—essentially, a model of the dynamic interdependence of all of the variables *X*(*t*) down the *right-hand side* of Fig. 1 (Brezina et al. 2003a). Joining this modeled “neuro-modulation transform” to a model of the basal, ACh-mediated neuromuscular transform—down the *left-hand side* of Fig. 1—gave a complete model of the modulated B15/B16-ARC NMT (Brezina et al. 2003b). For any input pattern of firing of the two motor neurons (*top* of Fig. 1), the model outputs the resulting modulated contraction shape (*bottom right* of Fig. 1). [This input–output relationship is, in essence, the definition of the NMT (Brezina et al. 2000b, c).] In particular, the model makes it possible to follow the modulated contraction waveform—and, of course, all of the internal modeled variables—as the motor neurons B15 and B16 fire in different patterns, as they do in the intact feeding animal, over hundreds of cycles that constitute a typical meal (Kupfermann 1974; Susswein et al. 1976, 1978; Weiss et al. 1982).

For adaptive behavior, however, what is important is not the contraction waveform per se, but rather how well it satisfies the functional requirements of the desired behavior (Brezina and Weiss 2000). As a final step, therefore, we will in our work here pass the contraction waveform through a realistically modeled behavioral task (schematized in Fig. 1, *bottom right*, and explained further in results) to evaluate its functional performance.

Running the model will immediately answer basic (although, given the complexity of the system, by no means trivial) questions such as: Given a firing pattern, what contractions and functional performance will it produce? Conversely, given a functional requirement, what firing patterns are optimal? Beyond this, however, the model now allows us to answer some general questions of motor control in this representative brain–body system. Some of these questions were raised in our previous work (see Brezina and Weiss 2000; Brezina et al. 2000b, c) but could then be addressed only conceptually or under unrealistically restricted conditions. Now, for the first time, we have the tools to answer them in a reasonably realistic fashion. Furthermore, for the first time, we can feed as input into the model an extended dataset of the real firing patterns of the motor neurons B15 and B16, recorded with chronic electrodes implanted in freely feeding animals by Horn et al. (2004). With these tools, we will in this paper focus, in particular, on the following questions.

*1)* Why is the additional layer of peripheral modulation necessary? Why cannot the motor neurons fire in such patterns as to immediately give the required contraction shapes through the basal NMT alone? The conceptual answer given before (Brezina et al. 2000c) is that, with fixed properties, the channel of the NMT is too narrow to allow even the entire set of firing patterns to create the broad range of contraction shapes that is required by the range of feeding behaviors in which the neuromuscular plant participates. Is this still true in the realistic model, where, for example, the modulators are no longer exogenously superimposed but rather are released endogenously by the activity of the system? Note that an affirmative answer requires demonstrating not merely that modulation improves the performance of a specific firing pattern, in particular the real pattern recorded in the feeding animal—although that is certainly a prerequisite—but furthermore that no other pattern exists that, without modulation, could have given a better performance.

*2)* The previous work (e.g., Brezina and Weiss 1997; Brezina et al. 1996, 2000b, 2003a, b) has suggested that the various components of the system—the two motor neurons, their modulators, and the different modulator effects—act in interdependence and mutual complementarity, one or another playing a key role depending on the circumstances. Furthermore, there is strong nonlinear coupling between the components, such that the functional role of each component depends on the context of the others. To what extent is this important for performance in the realistic model, with real motor neuron firing patterns?

*3)* The dynamics of the three main modulator effects have very different timescales. The modulation of K current is fast enough to react to each cycle of behavior, whereas the effects on Ca current and relaxation rate are slowly integrated over many cycles (Brezina et al. 2003a, b). Accordingly, it was proposed that whereas the former provides local adjustments, the persistence of the latter effects provides a peripheral memory of past behavior, which then prepares the system for future behavior of the same kind. However, these proposed roles—and other previous ideas—must be reexamined in light of the finding of Horn et al. (2004) that the real firing patterns of motor neurons B15 and B16, along with practically all other parameters of *Aplysia* feeding behavior, are not highly stereotyped as was assumed, but show very large, quasi-random variability from cycle to cycle.

*4)* A fundamental question then is: With stochastic central motor commands and slow, history-dependent dynamics in the periphery, to what extent is the system actually under control of the nervous system? When the nervous system sends motor neuron firing patterns through the NMT, to what extent can it predict the contractions and performance that will result? To what extent can it, through those patterns, control the various desired behaviors?

To answer these questions, we proceed in this paper in two complementary ways. First, as in the previous work, we run the model with a set of simple, regular motor neuron firing patterns. Although idealized, these patterns allow us to examine the complete space of firing pattern parameters as we must do if we are to demonstrate a global optimization of performance. Then, we use a dataset of the real, irregular firing patterns recorded by Horn et al. (2004). These patterns, conversely, are fully physiological but their parameters cannot easily be manipulated. Although each approach has its deficiencies, in combination they provide, as will be seen, considerable insight into motor control in the *Aplysia* B15/B16-ARC neuromuscular system and, we believe, other complex brain–body systems.

## METHODS

### Model

The equations and parameter values of the model were taken from Brezina et al. (2003a, b; also 1996, 2000a), and used without modification, except where explicitly stated. The major components of the model are shown in Fig. 1. Below, after an initial overview, each section gives first the model and its implementation as used with the ideal set of regular firing patterns in Figs. 2–9, then any variations used with the real dataset of irregular firing patterns in Figs. 10–19. The meanings of certain symbols—principally those that are often used in the results—are repeated in Table 1.

##### OVERVIEW.

The input to the model is the firing pattern—the spike patterns or, more precisely, the instantaneous firing frequency waveforms *f*_{B15}(*t*) and *f*_{B16}(*t*)—of the motor neurons B15 and B16 (Fig. 1, *top*). This input sets in motion two cascades of action that are modeled as initially independent of each other but ultimately joining.

First, the motor neuron firing pattern drives the “neuro-modulation” transform, down the *right-hand side* of Fig. 1. The firing patterns *f*(*t*) elicit release of the modulators SCP and MM, *r*(*t*). The released modulator amounts are converted to effective concentrations, *C*(*t*). The concentrations elicit the three main modulator effects: the enhancement of Ca current, **C**(*t*), the activation of K current, **K**(*t*), and the increase of the muscle's relaxation rate, **R**(*t*). The Ca-current and K-current effects then jointly determine the modulation of contraction size, **S**(*t*). All of these variables and their relationships are explicitly modeled as a function of the time *t.*

At the same time, the motor neuron firing pattern drives the basal NMT, down the *left-* *hand side* of Fig. 1. This is modeled without any internal details, simply as the direct conversion of the firing patterns *f*(*t*) to a basal muscle contraction waveform, *c*_{basal}(*t*) (Fig. 1, *bottom left*).

Finally, the outputs of the “neuro-modulation” transform, that is, **S**(*t*) and **R**(*t*), are applied to convert *c*_{basal}(*t*) to the modulated contraction waveform, *c*(*t*) (Fig. 1, *bottom right*). Altogether, for any input motor neuron firing pattern *f*_{B15}(*t*) and *f*_{B16}(*t*), the model outputs the corresponding modulated contraction *c*(*t*). It is thus a model of the complete modulated B15/B16-ARC NMT.

##### FIRING PATTERNS OF MOTOR NEURONS B15 AND B16.

As in Brezina et al. (2003a, b), we use here simple cyclical bursting patterns definable by the burst duration *d*_{intra}, the interburst interval *d*_{inter}, and the intraburst firing frequency *f*_{intra}, with cycle period *P* = *d*_{intra} + *d*_{inter}. We also define the duty cycle *D* ≡ *d*_{intra}/*P.* The instantaneous firing frequency *f*, as a function of time *t*, is then given by (1a) where *t**, the start time of the current cycle, is the largest number ≤*t* generated by the sequence (1b) where *N* is the cycle number, and *L* is the length of the entire sequence of firing which began at *t* = 0.

With the regular firing patterns, the parameter values *P, d*_{intra,B15} (and so, by *d*_{inter} = *P* − *d*_{intra}, *d*_{inter,B15}), *f*_{intra,B15}, *d*_{intra,B16} (and so *d*_{inter,B16}), and *f*_{intra,B16} were set and varied as described in results and figure legends. Specification of these parameters gave two particular functions of the form of *Eq. 1*, *f*_{B15}(*t*) and *f*_{B16}(*t*). With a common cycle period *P* and a common sequence of cycle start times *t**, the bursts of B15 and B16 firing always ended together at the end of each cycle, although, with different durations *d*_{intra,B15} and *d*_{intra,B16}, they could start at different times [see Fig. 3*B* of Brezina et al. (2003a)]. *L* was set as described below.

The dataset of irregular firing patterns was derived from a single in vivo chronic recording of electrical activity in the ARC muscle during a spontaneous, approximately 2½—h—long meal (see results and Horn et al. 2004), shown in Fig. 10*A*. Allocation of each spike in this extracellular record to either B15 or B16 (see results), followed by conversion of the spike times to firing frequencies by assigning to each timepoint in each interspike interval the reciprocal of the duration of that interspike interval, gave the instantaneous firing frequency functions *f*_{B15}(*t*) and *f*_{B16}(*t*) shown in Fig. 10*B*. These replaced *Eq. 1* as the primary input to the model in Figs. 12–19. For the purposes of analysis, we further partitioned the record of all spike times (before their allocation to B15 or B16) into individual cycles. As in *Eq. 1*, we assumed that there was one burst of spikes per cycle. Based on an examination of the distribution of interspike intervals in the record, we defined a burst as a continuous sequence of ≥10 spikes with interspike intervals of <1 s. Automatic implementation of this criterion identified 749 bursts in the 2½-h-long record. Note that application of this procedure before the allocation of B15 and B16 spikes made the B15 and B16 bursts in each cycle nominally identical in their timing and duration, in contrast to *Eq. 1* for the regular firing patterns. Between bursts, the interburst firing frequency for B15 and B16, jointly as well as individually, was by definition <1 Hz, and in most simulations was in fact explicitly set to zero—thus altering *f*_{B15}(*t*) and *f*_{B16}(*t*) somewhat—to allow implementation of the modulation of the muscle's relaxation rate (*Eqs. 6c* and *e* below).

*Equation 1a* implies that the interburst interval precedes the burst in each cycle, and we generally applied this rule—with the irregular as with the regular patterns—in computing, for example, cycle periods or duty cycles (Figs. 11, 14, 16). However, because the muscle contraction caused by a spike burst outlasted the burst considerably into the following interburst interval, it was more appropriate to take the burst together with the following interval as the cycle over which to evaluate quantities such as the contraction area, the functional performance of the contraction waveform (see below), and in some simulations (Fig. 17*A*) the corresponding period-averaged modulatory effects. This distinction between the preceding and the following interburst interval only became an issue with the irregular firing patterns in Figs. 13–18; with the regular firing in the steady state throughout Figs. 2–8, the distinction was irrelevant.

##### MODULATOR RELEASE.

Two modulators are considered: SCP released from B15 and MM released from B16. With either the regular or the irregular *f*(*t*) as input, the release *r* is described by (2a) (2b) (2c) where *S* is the size of the releasable pool of modulator, *p* is the probability of release of the modulator in *S*, *k _{p+}* and

*k*are the “on” and “off” rate constants of

_{p−}*p, y*is an exponent providing nonlinearity of release, and

*Q*

_{10}provides for dependency on the temperature

*T.*Here

*T*= 15°C and so

*Q*

_{10}

^{(T−15)/10}= 1. The other parameter values are (Brezina et al. 2000a, 2003a)

*k*

_{p+,B15}= 4.0 × 10

^{−10},

*k*

_{p+,B16}= 7.4 × 10

^{−11},

*k*

_{p−,B15}= 3.4 × 10

^{−3}s

^{−1},

*k*

_{p−,B16}= 6.6 × 10

^{−3}s

^{−1}, and

*y*= 3.

*Equation 2c*describes depletion of

*S*from an initial pool size

*S*

_{0,SCP}≡

*S*

_{SCP}(

*t*= 0) = 541 fmol and

*S*

_{0,MM}= 2690 fmol. For simplicity, in all simulations here we prevented depletion by setting d

*S*(

*t*)/d

*t*= 0 in

*Eq. 2c*and so

*S*(

*t*) =

*S*

_{0}.

##### MODULATOR CONCENTRATIONS.

With *r*(*t*) given by *Eq. 2* as input, the effective concentrations, *C*, of SCP and MM are described by (3) where ν is the effective volume of the released modulator in the muscle and *k _{C}* is the rate constant of removal of the modulator from ν. The parameter values are (Brezina et al. 2003a) ν

_{SCP}= 10 μl, ν

_{MM}= 35 μl,

*k*

_{C,SCP}= 0.1 s

^{−1}, and

*k*

_{C,MM}= 0.35 s

^{−1}.

##### MODULATOR EFFECTS.

Three independent effects are considered: enhancement of Ca current, **C**; activation of K current, **K**; and increase of the muscle's relaxation rate, **R**. With *C*(*t*) given by *Eq. 3* as input, each effect **X** ∈ {**C**, **K**, **R**} is described by (4a) (4b) where *C*′_{i} ≡ *C*_{i}^{hX,i} is the concentration of the modulator *C*_{i} modified by a Hill coefficient *h*_{X,i}, **X**_{max,i} is the maximal effect of the modulator, *k*_{X+,i} and *k*_{X−,i} are “on” and “off” rate constants, and *F*_{ji} specifies the fraction of the effect of modulator *j* that occludes the effect of modulator *i*. The parameter values are (Brezina et al. 2003a) **C**_{max} = 80%, *k*_{C+} = 750 s^{−1} M^{−hC}, *k*_{C−} = 2.25 × 10^{−3} s^{−1}, *h*_{C} = 0.75, *F*_{C} = 1; **K**_{max,SCP} = 27%, **K**_{max,MM} = 100%, *k*_{K+} = 1.9 × 10^{5} s^{−1} M^{−hK}, *k*_{K−,SCP} = 3.5 s^{−1}, *k*_{K−,MM} = 0.35 s^{−1}, *h* _{K} = 0.85, *F*_{K,SCP,MM} = 1, *F*_{K,MM,SCP} = 0.27; **R**_{max} = 100%, *k*_{R+} = 305 s^{−1} M^{−hR}, *k*_{R−} = 1.11 × 10^{−3} s ^{−1}, *h*_{R} = 0.7, *F*_{R} = 1.

In one set of simulations (Fig. 17*A*) we altered the rate constants of the effects, either dividing *k*_{K+} and *k*_{K−} by 375 to make **K** approximately as slow as **C** and **R**, or multiplying *k*_{C+} and *k*_{C−} by 250, or *k*_{R+} and *k*_{R−} by 500, to make **C** or **R** as fast as **K**.

##### MODULATION OF CONTRACTION SIZE.

**C**(*t*) and **K**(*t*) given by *Eq. 4* combine to determine the modulation of contraction size **S**, according to the equation (5) with the parameter values **S**_{max} = 400%, *K*_{S} = 8%, κ = 9%, and *a* = 8/3% (see Brezina et al. 1996, 2003a).

##### THE MODULATED CONTRACTION WAVEFORM.

The instantaneous modulated contraction amplitude *c* is modeled as a basal contraction amplitude *c*_{basal}, then modulated by **S**(*t*) and **R**(*t*), according to the equations (6a) (6b) (6c) (6d) (6e) where *c*_{basal,∞}(*f*), *r*_{contr}(*f*), and *r*_{relax}(*c*) are empirical functions plotted in Fig. 1, *B–D* of Brezina et al. (2003b), *t** is the start time of the current cycle given by *Eq. 1b* with the regular *f*(*t*) or the end of the previous burst with the irregular *f*(*t*), δ_{R} = 0.03, and (6f)

When *f*_{NMT}(*t*) > 32 Hz or *c*(*t*) > 1, the functions *c*_{basal,∞}(*f*); *r*_{contr}(*f*); and *r*_{relax}(*c*) were extrapolated. For further details see Brezina et al. (2003b).

We routinely ran parallel simulations with modulation, as above, and without modulation, setting **S**(*t*) = 100 and **R**(*t*) = 0 in *Eqs. 6d* and *6e*. In some simulations (Figs. 15–17) we disabled one or more of the effects **X** ∈ {**C, K, R**} by setting **X**(*t*) = 0 in *Eq. 4a*, or alternatively **S**(*t*) = 100 or **R**(*t*) = 0 in *Eq. 6d* or *6e*, as appropriate.

### Performance measure

To evaluate the functional significance of the contraction waveform *c*(*t*), we use the performance measure *m*, from Task VI of Brezina and Weiss (2000) (see their Appendix F), except modified for the single ARC muscle. Here we give just the essential implementation details; for the biological motivation of the task and performance measure, see results.

The contraction waveform in any cycle *n*, [*c*(*t*)]_{n}, is required to alternately rise above and fall below a threshold amplitude, *c* = , during particular phases of the cycle (see Fig. 1). We define a template of 4 timepoints *t1* = *t** + φ, *t2* = *t** + φ + δ, *t3* = *t** + φ + λ, and *t4* = *t** + φ + δ + λ, where *t** is the start time of the current cycle and φ, δ, and λ are time intervals ≤*P*. With the regular firing patterns, we set δ = 0.15*P* and λ = 0.5*P* so as to create in each cycle 2 evenly spaced windows of width 0.15*P* (Figs. 1, 2*B*). We then require that *c*(*t*1) < , *c*(*t*2) > , *c*(*t*3) > , and *c*(*t*4) < , and if these basic requirements are satisfied, that *c*(*t*1) and *c*(*t*4) should be as far below and *c*(*t*2) and *c*(*t*3) as far above as possible. A set of equations that incorporates all of these requirements is (7a) (7b) (7c) (7d) (7e) (7f) [See the Appendices of Brezina and Weiss (2000) for the development of these equations. They are simply computationally compact forms of the required conditions. *Equation 7a*, for example, is equivalent to *A*_{1} = {* − **c*(*t*1) if *c*(*t*1) < *, 0 if **c*(*t*1) ≥ *.] In **Eq. 7f*, we vary φ from 0 to *P* to slide the template of the timepoints *t1–t4* along [*c*(*t*)]_{n} so as to maximize *m*′_{n}(φ). This yields the basic performance measure *m*_{n}. In this paper, however, we work throughout with the modified performance measure *M*_{n} ≡ *m*_{n}/*P* 〈*c*〉* _{n}*, where 〈

*c*〉

*is the mean, period-averaged contraction amplitude, and so*

_{n}*P*〈

*c*〉

*is the area under the contraction waveform over*

_{n}*P*, in the cycle

*n.*

In most simulations with the regular firing patterns, *P* ≤ 10 s and the performance algorithm was used exactly as above. When *P* > 10 s (*maps 9* and *10* of Figs. 2*A*, 4*A*, and 8), the relative spacing of the timepoints *t1–t4* was set as for *P* = 10 s. Throughout, we used the threshold amplitude = 0.5.

With the irregular firing patterns, we derived the window interval δ on a cycle-by-cycle basis from the burst duration *d*_{intra}, setting δ = 0.3*d*_{intra}. The other window interval λ was not defined; like φ, it was implicitly found during the maximization procedure, which was broken into two parts. First we varied φ from 0 to *P* to maximize *m*′_{n} (φ) = *A*_{1}*A*_{2}, say at φ = φ*; then we varied λ from φ* to *P* to maximize *m*′_{n} (λ) = *A*_{3}*A*_{4}. Finally (7g) (see text footnote 3 in results). As threshold amplitude, we used = 0.5, = 0.132, or (in Fig. 18*B*) randomly selected in each cycle from the positive part of a Gaussian distribution with mean and standard deviation both equal to 0.132 (see results).

### Simulation logistics

##### STEADY-STATE CONTRACTIONS AND PERFORMANCE WITH REGULAR FIRING PATTERNS.

With motor neuron firing of constant pattern—i.e., *Eq. 1a* with constant *d*_{intra}, *d*_{inter}, *f*_{intra}, and *P*—we ran each simulation for a time, *L* in *Eq. 1a*, that was sufficiently long to allow the contraction waveform to approach close to the steady state. We then measured performance in the last 1–2 cycles.

To obtain such steady-state contractions in the absence of modulation (e.g., maps in the *top rows* of Figs. 2*A* and 4*A*, and *maps 1–4* of Fig. 6), it was necessary to run the simulations for only several cycles [see, e.g., Figs. 7*D* and 8*F* of Brezina et al. (2003b)]. We generally used *L* = 30–100 s, and in any case *L* ≥ 5*P.*

To obtain steady-state contractions with endogenous modulator release (e.g., maps in the *bottom rows* of Figs. 2*A* and 4*A*, *maps 5–8* of Fig. 6, and Fig. 8), we used *L* = 1 h. With constant firing pattern and no modulator depletion, after 1 h the system should be very close to the steady state [cf. Fig. 14 of Brezina et al. (2003a)], and, by monitoring the trajectories of **S**(*t*) and **R**(*t*), we confirmed this for the majority of the simulations. Only in a few simulations with small *f*_{intra} were **S**(*t*) and **R**(*t*) still noticeably changing after 1 h; all these cases were negligible as they had zero performance.

For each of the two motor neurons B15 and B16, we mapped in these simulations the performance as a function of two parameters, *d*_{intra} and *f*_{intra}, at each *P.* When only one of the motor neurons was studied alone, it was possible to completely, albeit coarsely, map the entire 2-dimensional *d*_{intra} × *f*_{intra} space (up to a reasonable value of *f*_{intra}) using a 10 × 10 grid, i.e., 100 separate simulations with *L* = 1 h each (Figs. 2*A*, 4*A*). (With much more computer time, even a 30 × 30 map was possible: Figs. 3*B*, 5*B*.) Studying both motor neurons together, however, it was not realistically possible to exhaustively map the entire 4-dimensional *d*_{intra,B15} × *f*_{intra,B15} × *d*_{intra,B16} × *f*_{intra,B16} space, of 10,000 simulations, in this way. Instead, we ran the simulations, each with *L* = 1 h, only at selected locations to search the space (up to a reasonable value of *f*_{intra}, 30 Hz for B15 and 40 Hz for B16) for the maximal performance. In the first stage of the search, we sampled the space randomly to locate candidate regions of high performance. In the second stage, starting from the point of highest performance found in the first stage, we used gradient ascent, varying each of the 4 parameters in turn by 1% of its range, to locate the point of maximal performance. Note that this strategy was not guaranteed to find the global maximum, although in many cases—given the simple, connected, relatively smooth structure of at least some of the spaces involved (e.g., Figs. 2–7)—it may well have done so. Finally, we mapped 2-dimensional 10 × 10 sections, as for the individual motor neurons, through the point of maximal performance (e.g., maps in the *bottom row* of Fig. 6; Figs. 7*A*, 8).

In simulations to map the steady-state performance with constant exogenous modulator concentrations (Figs. 3*A* and 5*A*), we analytically precomputed the steady-state modulator effects so as to be able to shorten *L.* With constant concentration *C _{i}*(

*t*) =

*C** of a single modulator

_{i}*i*and, in the steady state, d

**X**

*(*

_{i}*t*)/d

*t*= 0,

*Eq. 4b*yields (8) where

*K*

_{x}

*≡*

_{,i}*k*

_{x−,i}/

*k*

_{x+,i}, for each effect

**X**∈ {

**C, K, R**}. Setting

**R**(

*t*) =

**R**

_{i,∞}and (converting

**C**

_{i,}_{∞}and

**K**

_{i,}_{∞}to

**S**

_{i,}_{∞}using

*Eq. 5*)

**S**(

*t*) =

**S**

_{i,}_{∞}in

*Eqs. 6d*and

*6e*allowed the contractions to come to steady state as fast as unmodulated contractions (i.e., within

*L*= 30 s).

##### TRANSIENT CONTRACTIONS AND PERFORMANCE WITH REGULAR FIRING PATTERNS.

In all transient simulations we started the model from the “resting” set of initial conditions (*p*_{B15} = 0, *p*_{B16} = 0, *S*_{SCP} =*S*_{0,SCP}, *S*_{MM} = *S*_{0,MM}, *C*_{SCP} = 0, *C*_{MM} = 0, **C**_{SCP} = 0, **C**_{MM} = 0, **K**_{SCP} = 0, **K**_{MM} = 0, **R**_{SCP} = 0, **R**_{MM} = 0, *c*_{basal} = 0). We fired the motor neurons in a constant pattern (e.g., Fig. 9*A*, records of *f*_{B15} and *f*_{B16} at the *bottom*) while measuring performance (Fig. 9*A*, *top*) in each successive cycle composed of burst plus the following interburst interval (see above). At selected times we halted the simulation at the current values of all model variables, then mapped the performance in the next cycle over the *d*_{intra,B15} × *f*_{intra,B15} × *d*_{intra,B16} × *f*_{intra,B16} space (Fig. 9*B*) using the same strategy of search and mapping through the point of maximal performance as in the steady-state simulations.

##### IRREGULAR FIRING PATTERNS.

In all simulations, the irregular patterns were run through the model starting from the “resting” set of initial conditions. The same was done with the regular patterns, but irregular task, in Fig. 18*B*.

### Analytical and numerical methods

Motor neuron B15 and B16 spike times were automatically tabulated from the 2½-h record of ARC electrical activity using the threshold event detection module of Clampfit 9 (Axon Instruments, Union City, CA). Subsequent processing and analysis of the irregular dataset, and all analysis of regular data, was done in *Mathematica* (Wolfram Research, Champaign, IL) or during the preparation of figures in SigmaPlot (SPSS, Chicago, IL). The Kolmogorov–Smirnov test was run in Clampfit, other statistical tests in *Mathematica* or SigmaStat (SPSS).

In all regular as well as irregular simulations, the model equations were integrated numerically using *Mathematica*, which implements a nonstiff Adams method and a stiff Gear method, based on the LSODE routine (Wolfram 1996). These implementations attempt to constrain the local error in a numerical solution of size *x* to be smaller than 10^{−a} + |*x*| 10^{−p}, where *a* is the number of digits of accuracy and *p* is the number of digits of precision. Standard settings were *a* = 6 and *p* = 6, with a maximal step size of 10 ms. When solutions were expected to be absolutely very small, either the equations were rescaled or *a* was increased so as to maintain the local relative error ≤10^{−5} in most cases. Results with increased *a* and *p* and decreased maximal step size were routinely compared to confirm global accuracy.

## RESULTS

### The model, functional task, and performance measure

Figure 1 summarizes our model of the modulated B15/B16-ARC NMT and the way in which we use it in this paper. The input to the model is the firing pattern—the spike patterns or, more precisely, the instantaneous firing frequency waveforms *f*_{B15}(*t*) and *f*_{B16}(*t*)—of the motor neurons B15 and B16. For any such input, the model outputs the resulting modulated contraction waveform *c*(*t*). (The meanings of frequently used symbols are listed in Table 1 for convenience.) Mathematically, the model consists of a set of coupled differential equations and associated algebraic equations, fully set out in methods. Because of their complexity, these equations have no analytical solution but must be integrated numerically. Throughout the paper, we run the model in this way for various kinds of input under various conditions.

Having obtained the contraction waveform, we pass it through a modeled functional task to evaluate its performance. The task, developed in our earlier work (Brezina and Weiss 2000), is motivated by the basic kinematics of *Aplysia* consummatory feeding movements. Briefly, in an ingestive cycle (biting, swallowing), the open radula protracts from the mouth, closes on food, retracts, and opens again to release the food into the esophagus. In an egestive (rejection) cycle, similar movements occur but the radula is closed during protraction rather than retraction, pushing material out rather than pulling it in (Kupfermann 1974; Morton and Chiel 1993). In either case, the radula must alternately close and open in each successful cycle, meaning that the ARC muscle, which contributes to the closure (Cohen et al. 1978; Orekhova et al. 2001), must alternately contract and relax. In our formal definition of the task, we require the contraction waveform *c*(*t*) to alternately rise above and fall below some threshold amplitude . The closure and opening of the radula must, furthermore, be coordinated with its protraction and retraction in each cycle. The phase relationship of the two sets of movements is opposite in ingestive and egestive cycles, but what is common to all cycles, from the perspective of the ARC muscle, is that its contraction and relaxation must occur at definite times during the cycle. Formally, we require each of the two transitions across the threshold to occur within a relatively narrow window in each cycle (e.g., first between *t1* and *t2*, and then between *t3* and *t4*, in Fig. 1, *bottom right*). Finally we assume that, provided all of these requirements are satisfied, the larger the movements are, the better. Formally, we combine measurements of the contraction amplitude *c*, relative to the threshold , at the beginning and end of each transition window (two-headed vertical arrows at *t1–t4* in Fig. 1, *bottom right*). In the methods we give a set of equations that embodies all of these functional requirements. These equations yield in each cycle a single-valued basic performance measure, *m.* Throughout this paper, however, we present *m* already in the form of the higher-order performance measure *M* ≡ *m/P*〈*c*〉, where *P* is the cycle period, 〈*c*〉 is the mean, period-averaged contraction amplitude, and *P*〈*c*〉 is therefore the area under the contraction waveform over *P.* This modification attempts to include in the performance measure also the fact that more frequent repetition of the same movement in a given time interval represents better performance, as does a lower, less energetically costly contraction waveform. *M* can be thought of as an elementary measure of functional efficiency (Brezina and Weiss 2000). *M* is zero if the basic requirements of the task are not met at all; once they are, increasing *M* indicates increasing performance or efficiency.

In this task, the threshold amplitude is relatively undetermined. Different values of define, in effect, a family of tasks with different degrees of difficulty. We begin by setting = 0.5, a moderate value relative to the actual contraction amplitude *c*, which varies from 0 to >2 (see, e.g., Figs. 2*B* and 12 below). Later we examine some of the consequences of altering .

### First analysis with regular motor neuron firing patterns

We now ask, does modulation in fact improve performance? And, if it does, how do the various components of the system contribute to this result? To develop an initial understanding, we proceed stepwise. First, we analyze how the endogenous release of SCP from motor neuron B15 modifies performance achieved by B15 alone. Then we analyze how the release of MM from B16 modifies the performance of B16, and finally, how the two motor neurons act together. We use a simple set of regular firing patterns—idealized regular bursts of spikes—that are completely defined by just three parameters, such as their burst duration *d*_{intra}, intraburst firing frequency *f*_{intra}, and cycle period *P* (Brezina et al. 1997, 2000b). To begin with, we run each pattern through the model until the model reaches steady state,^{1} and compute, therefore, the steady-state performance *M*_{∞}.

### SCP release from motor neuron B15 improves performance

In Fig. 2*A* we have constructed maps of *M*_{∞} over the space of parameters defining the regular firing patterns of motor neuron B15. The complete parameter space is 3-dimensional: *d*_{intra,B15} × *f*_{intra,B15} × *P.* Each map covers, on a 10 × 10 grid, two of these dimensions, *d*_{intra,B15} × *f*_{intra,B15}, and the third dimension, *P*, varies across the maps from *left* to *right.* On the *x*-axis of each map, *d*_{intra,B15} ranges from 0 to *P;* on the *y*-axis, *f*_{intra,B15} ranges from 0 to 30 Hz, which is more than double the firing frequency that B15 is likely to sustain in behavior (Church and Lloyd 1994; Cropper et al. 1990a; see Fig. 11 below). The values of *P* chosen cover the range from 1 to 10 s, which includes the values typically recorded within the relatively rapid bouts into which the feeding cycles of a meal are usually grouped (e.g., Hurwitz and Susswein 1992; Kupfermann 1974; Susswein et al. 1976, 1978; Weiss et al. 1986; see Fig. 11). *P* = 100 s has been added to represent the intervals between bouts (consequently, this case is treated as that of *P* = 10 s, but with a 90-s gap, in setting the range of *d*_{intra} and the time windows of the functional task). Altogether, then, these maps cover, albeit in a discrete and relatively coarse fashion, the complete space of parameters that is likely to be functionally relevant. Within each map, each pixel represents a separate simulation, whose steady-state performance *M*_{∞} is indicated by the color of the pixel, ranging through a spectrum from blue (*M*_{∞} = 0) to red (*M*_{∞} = 0.121, the highest value in this set of simulations). In the *top row* of maps in Fig. 2*A*, we ran the model without modulation, whereas in the *bottom row* we allowed motor neuron B15 to release its endogenous SCP.

We see that nonzero performance is achieved only over a small subset of parameter values. In each map these values are found in a rough arc (seen better in Fig. 3 where the maps have a finer, 30 × 30 grid), extending from relatively small *d*_{intra} but large *f*_{intra} to small *f*_{intra} but large *d*_{intra}. Elsewhere in the space the requirements of the functional task are not even minimally satisfied.

The limited region of nonzero performance, and its shaping into an arc by a trade-off between *d*_{intra} and *f*_{intra}, are both consequences of the fact that, for success in this task, the phasic oscillation in the contraction waveform must straddle the threshold amplitude . This can be appreciated in Fig. 2*B*, which shows six representative contraction waveforms, *a–f*, from the locations *a–f* in Fig. 2*A*, at *P* = 3 s. (All of this argument applies, in a more self-evident fashion, also to nonsummating contractions at larger values of *P.*) Roughly, the mean contraction amplitude 〈*c*〉, around which the waveform oscillates, must be comparable to [this indeed becomes rigorously true as *P* becomes very small (Brezina and Weiss 2000)]. Thus, waveforms *b* and *d* in Fig. 2*B* are functional, whereas waveform *a* oscillates too low, and waveforms *c, e*, and *f* oscillate too high, to be functional. The overall height of the oscillation—and, more rigorously, 〈*c*〉 (Brezina et al. 2000b)—increases with both *d*_{intra} (Fig. 2*B*) and *f*_{intra} (see Fig. 3*C* below). If one of these parameters makes the oscillation too high or too low, the other can be adjusted in the opposite direction to compensate and return the system to the arc of functional performance. Above and to the right of the arc in each map, at large *d*_{intra} and *f*_{intra}, the oscillation is too high; below and to the left of the arc, at small *d*_{intra} and *f*_{intra}, the oscillation is too low. The arc attenuates with decreasing *P* as the oscillation gets faster, smaller, and less likely to precisely straddle (see Brezina and Weiss 2000; Brezina et al. 2000b).

When motor neuron B15 fires to elicit basal contractions but does not release SCP to modulate them, performance is generally low (*top row* of maps in Fig. 2*A*). When the motor neuron also releases SCP, performance improves significantly (*bottom row* of maps). The effect is smallest at *P* = 100 s, that is, with infrequent bursts of firing that release little SCP in the model (Brezina et al. 2000a, 2003b) as in reality (Cropper et al. 1990b; Vilim et al. 1996a; Whim and Lloyd 1989, 1990). With frequent bursts at *P* < 10 s, however, performance improves 2-fold or more.

This improvement in performance can be explained by considering the effects of SCP on the muscle. Of the three main modulator effects described in the introduction, SCP exerts principally two: it enhances the muscle's Ca current to potentiate the size of the contractions, and it increases their relaxation rate, in the model (Brezina et al. 2003a, b) as in reality (e.g., Brezina et al. 1994a; Cropper et al. 1987a; Lloyd et al. 1984; Probst et al. 1994; Whim and Lloyd 1990). The first effect can also be viewed as an increase in the contraction rate of the muscle (Brezina et al. 2000c) so that both effects accelerate the NMT. The result, as Fig. 2*B* shows, is that the phasic oscillation of the contraction waveform is larger and so, provided that its mean 〈*c*〉 continues to approximate , gives better performance. However, the overall height of the waveform tends to increase at the same time, changing 〈*c*〉 in relation to . Waveforms that previously were of optimal height now become too large (compare waveforms *b* and *e* in Fig. 2*B*), whereas waveforms that were too small may now become optimal (waveforms *a* and *d*). Hence we see the arc of functional performance shift to noticeably smaller *d*_{intra} and *f*_{intra} in the maps in Fig. 2*A*.

In this complete model of the NMT, the modulation is no longer exogenous, superimposed from outside irrespective of the motor neuron firing pattern, but rather is endogenously set in motion by the same firing pattern that elicits the basal contraction waveform, to which it is thereby partially coupled (Brezina and Weiss 1997; Brezina et al. 2000c). This creates quite severe constraints that must be satisfied. The firing pattern must elicit through the basal NMT a basal contraction waveform, and simultaneously through the “neuro-modulation transform” a modulatory state, such that, when the modulatory state then operates on the basal contraction waveform, it converts it into a modulated contraction waveform whose shape satisfies the functional requirements. To see how this happens, in Fig. 3 we dissect the components of this mechanism.

Figure 3*A* shows maps of the steady-state performance *M*_{∞} over the parameter space *d*_{intra,B15} × *f*_{intra,B15} (all at *P* = 5 s) with simulated exogenous application of increasing concentrations of SCP, *C*_{SCP}, which, for simplicity, we here let stand for the complex multidimensional modulatory state (see later). In these maps, we can thus look up what performance would be achieved by any firing pattern (*d*_{intra,B15}, *f*_{intra,B15}) with any SCP concentration, irrespective of whether that concentration can in fact be endogenously attained by that firing pattern. Now, can it in fact be so attained? *Map 7* in Fig. 3*B* shows the mean steady-state SCP concentration, 〈*C*_{SCP}〉_{∞}, that is attained by the firing patterns over the same space *d*_{intra,B15} × *f*_{intra,B15}. We can now combine *map 7* with the maps in Fig. 3*A* to simply look up the performance achieved by the complete, endogenously modulated contraction waveform. This performance is shown in *map 6* of Fig. 3*B*. Consider, for example, the black triplet of connected points in Fig. 3, *A* and *B*. This firing pattern, (*d*_{intra,B15} = 4.2 s, *f*_{intra,B15} = 7 Hz), attains 〈*C*_{SCP}〉_{∞} ≈ 100 nM in *map 7.* When we look the pattern up in the corresponding, 100-nM *map 3*, however, we find that *M*_{∞} = 0. Therefore *M*_{∞} = 0 in *map 6.* Consider, on the other hand, the red triplet of points. This pattern, (*d*_{intra,B15} = 2.2 s, *f*_{intra,B15} = 17 Hz), attains 〈*C*_{SCP}〉_{∞} ≈ 1 μM, and when we look it up in the corresponding 1-μM *map 4*, *M*_{∞} > 0; therefore *M*_{∞} > 0 in *map 6.* Thus the red pattern, unlike the black pattern, succeeds in activating both the basal and the modulatory components of the NMT in such a way that together they achieve functional performance.

Figure 3*C* plots the unmodulated and modulated steady-state mean contraction amplitude 〈*c*〉_{∞} as a function of *f*_{intra,B15} (for one representative value of *d*_{intra,B15}, specifically along the white lines in *maps 1* and *6* of Fig. 3, *A* and *B*). The colored squares mark the points where *M*_{∞} > 0. In further illustration of our analysis above, we see that functional performance is achieved only around the point where the 〈*c*〉_{∞}–*f*_{intra} plot crosses the line of threshold amplitude , that is, where 〈*c*〉_{∞} ≈ . Modulation steepens the slope of the plot and so shifts this point to smaller *f*_{intra}. What is most important, however, is that SCP acts so that it always preserves the positive slope of the plot, so that the plot always crosses the line of at some *f*_{intra}. (Exactly the same argument applies to *d*_{intra}.) This guarantees that, with SCP modulation, solutions to the functional problem—firing patterns that achieve functional performance—always exist. Fundamentally, the components of the SCP-modulated B15 system are mutually stabilizing. Through the basal NMT, increasing *f*_{intra} or *d*_{intra} increases contraction amplitude; through the “neuro-modulation transform,” the increasing *f*_{intra} or *d*_{intra} increases the release of SCP; increasing SCP also increases contraction amplitude. All of these relationships have a positive sign. Hence, if we picture optimization of performance as an iterative, quasi-dynamical process—iteratively increasing or decreasing *f*_{intra} or *d*_{intra} until the best solution is found—the components all converge on the same, stable solution.

### MM release from motor neuron B16 reduces performance

We now carry out the same analysis for motor neuron B16, releasing MM. Figure 4 is the counterpart of Fig. 2, and Fig. 5 of Fig. 3. Because B16 typically fires at higher frequency than B15, we study *f*_{intra,B16} up to 40 Hz, which again is about double the frequency that B16 is likely to sustain in behavior (Church and Lloyd 1994; Cropper et al. 1990a; see Fig. 11).

The unmodulated contractions elicited by motor neuron B16 are essentially identical to those elicited by B15,^{2} as is therefore also their low, but measurable, functional performance (*top row* of maps in Fig. 4*A*; contraction waveforms *a–c* in Fig. 4*B*). When motor neuron B16 is allowed to release its endogenous MM, however, all functional performance (except at *P* = 100 s where little MM is released) completely vanishes (*bottom row* of maps in Fig. 4*A*).

This can be understood by considering the effects of MM on the muscle. Like SCP, MM enhances the muscle's Ca current to potentiate the size of contractions and it increases their relaxation rate. Unlike SCP, however, MM also exerts to a significant degree the third main modulator effect: in the model (Brezina et al. 2003a, b) as in reality (Brezina et al. 1994b, 1995; Cropper et al. 1991; Orekhova et al. 2003), MM activates in the muscle a K current that depresses the size of contractions. As the MM concentration increases above a certain level—in the range between 10 and 100 nM (Fig. 5*A*) —the depressing effect begins to predominate. When contractions are depressed so that they fail even to reach the threshold amplitude (waveforms *d–f* in Fig. 4*B*; Fig. 5*C*), performance vanishes.

Fundamentally, the MM-modulated B16 system fails to satisfy the simultaneous constraints on the basal contraction, the modulatory state, and their interaction. As *map 6* of Fig. 5*B* shows, there appears to be, in the entire set of firing patterns, no pattern at all that can (in the steady state) satisfy all of the constraints simultaneously. Many patterns (e.g., the black triplet of points in Fig. 5, *A* and *B*) raise the mean steady-state MM concentration, 〈*C*_{MM}〉_{∞}, to ≥100 nM where all contractions, elicited by all patterns, are depressed below (*maps 3–5* of Fig. 5*A*). Even those patterns (e.g., the gray triplet of points) that raise 〈*C*_{MM}〉_{∞} to only 10 nM, where some patterns still give functional performance (*map 2*), fail because these patterns that raise 〈*C*_{MM}〉_{∞} to only 10 nM are not the same patterns that give functional performance at 10 nM. The former attain the low 〈*C*_{MM}〉_{∞} by having a small *d*_{intra} or *f*_{intra}, but with this small *d*_{intra} or *f*_{intra} they do not elicit a large enough basal contraction. The latter have large enough *d*_{intra} and *f*_{intra} to give contraction and performance at a low 〈*C*_{MM}〉_{∞}, but, with larger *d*_{intra} and *f*_{intra}, they attain a much higher 〈*C*_{MM}〉_{∞}.

In sum, we can say that the components of the MM-modulated B16 system are mutually destabilizing because their relationships are of opposite signs. As in the B15 system, increasing *d*_{intra} or *f*_{intra} increases basal contraction amplitude; increasing *d*_{intra} or *f*_{intra} increases the release of MM; but increasing MM, above a certain point, *decreases* contraction amplitude. Hence, if *d*_{intra} or *f*_{intra} is altered in an attempt to move the basal contraction amplitude toward , the action of MM opposes that movement, and indeed, this effect predominates. (With increasing *d*_{intra}, this can be seen in the waveforms *d–f* of Fig. 4*B*; with increasing *f*_{intra}, in the negative slope of the modulated 〈*c*〉_{∞}–*f*_{intra} plot in Fig. 5*C*.) Thus the components of the system do not converge on one stable solution, but rather diverge. With = 0.5, indeed, the modulated 〈*c*〉_{∞}–*f*_{intra} plot never crosses the line of (Fig. 5*C*) and no solution at all can be found.

### Simultaneous release of SCP and MM combines to optimize performance

Thus, release of SCP from motor neuron B15 improves functional performance, whereas release of MM from motor neuron B16 reduces, even eliminates, performance. What happens when both modulators are released simultaneously in the complete system?

A difficulty that we encounter at this point is that, with both motor neurons B15 and B16, the complete space of firing patterns, at each single value of *P*, is 4-dimensional: *d*_{intra,B15} × *f*_{intra,B15} × *d*_{intra,B16} × *f*_{intra,B16}. [The phase between the firing of B15 and B16 is fixed (see methods).] This 4-dimensional space is unmanageable graphically. We therefore present appropriate sections through the space. In the *bottom row* of Fig. 6 and in subsequent figures, wherever a pair of maps of *M*_{∞} is enclosed in a box, these are simultaneous 2-dimensional sections through the complete 4-dimensional space. The “B15 with B16” map plots *M*_{∞} over the B15 parameters *d*_{intra,B15} × *f*_{intra,B15}, and the “B16 with B15” map over the B16 parameters *d*_{intra,B16} × *f*_{intra,B16}, as in our analysis of each motor neuron alone except that here the other motor neuron was present with its parameters fixed at the values that, in a search through the complete 4-dimensional space, gave the maximal *M*_{∞} (see methods). In this way, although we do not show and indeed did not compute *M*_{∞} at every point in the complete 4-dimensional space, we show representative sections that, furthermore, include the region of maximal *M*_{∞} as required for any demonstration of global optimization.

Figure 6 compares the steady-state performance *M*_{∞} that is achievable in the different subsets of the B15/B16-ARC system—B15 alone with and without modulation, B16 alone with and without modulation, both B15 and B16 without modulation, and finally in the complete system, both B15 and B16 with modulation—at one representative value of *P, P* = 5 s. Recapitulating the results already presented, when B15 (*map 1*) or B16 (*map 2*) fires alone without modulation, performance is low; performance is somewhat improved by SCP release from B15 (*map 5*) but eliminated by MM release from B16 (*map 6*). Thus when each motor neuron fires alone, adding its modulation does not improve performance greatly. Similarly, when both B15 and B16 fire together but without modulation (*maps 3* and *4*), performance is not greatly improved over that achieved by each motor neuron alone. When, however, the complete system is assembled—both B15 and B16, releasing both SCP and MM—performance is dramatically improved over that achieved by any of the subsets of the system (*maps 7* and *8*). For example, the maximal performance is 4-fold better than the maximal performance achieved with both B15 and B16, but without modulation. Remarkably, in the complete system the release of MM, far from eliminating performance, instead improves it.

Figure 7*A* shows *M*_{∞} in another 2-dimensional section through the 4-dimensional space, that is, the section *f*_{intra,B15} × *f*_{intra,B16}, without modulation (*map 1*) and modulated by both SCP and MM (*map 2*). Without modulation, the diagonal region of functional performance is shaped by the essential equivalence of *f*_{intra,B15} and *f*_{intra,B16} as inputs to the basal NMT (*Eq. 6f* in methods, plotted as the thick diagonal line in Fig. 7*A*). Modulation greatly increases performance, and at the same time strikingly reorients the region where performance is achieved.

Figure 7*B* shows a series of representative modulated contraction waveforms, *a–e*, from the locations *a–e* in Fig. 7*A*. Figure 7*C* then plots the mean steady-state contraction amplitude 〈*c*〉_{∞} along the thin horizontal lines in Fig. 7*A*, unmodulated and modulated, the latter including the locations *a–e*, as *f*_{intra,B15} increases at a constant *f*_{intra,B16}. Figure 7*C* is the counterpart of the earlier Figs. 3*C* and 5*C* and gives an intuitive idea of the mechanism by which the complete modulation, by both SCP and MM, acts to improve performance. As *f*_{intra,B15} increases, 〈*c*〉_{∞} rises toward the threshold amplitude , but then does not merely cross , as with SCP alone (Fig. 3*C*), or decline from again, as with MM alone (Fig. 5*C*), but rather continues roughly along the level of as *f*_{intra,B15} continues to increase. At the same time the contraction changes its shape around 〈*c*〉_{∞} (Fig. 7*B*), becoming faster for improved performance. Thus the potentiating action of SCP and the depressing action of MM balance so as to maintain contraction amplitude relatively constant, allowing the common effect of the two modulators, the acceleration of the NMT, to develop for greatly improved performance over an extended region of the space of firing patterns.

The improvement in performance is seen not just at the single value of *P* = 5 s, but at all behaviorally relevant values of *P* (Fig. 8). With larger values of *P*, performance is improved less but over a larger region of the space; with smaller values of *P*, it is improved more but over a smaller region of the space (the white pixels in Fig. 8 are spots of outstandingly high, essentially offscale, performance).

### Different dynamical time scales uncouple the modulatory state from the basal contraction

So far we have considered only the steady state, when, after sufficiently long repetition of the same input firing pattern, all of the components of the B15/B16-ARC NMT have stabilized at steady values and are coupled to each other in a unique way for each firing pattern. As we observed, this coupling significantly constrains the functional solutions that can be found. We will see shortly, however, that the real firing patterns of the motor neurons B15 and B16 vary considerably from one cycle to the next. Moreover, the components of the NMT on which the patterns act have diverse dynamical time scales that can be much longer than the cycle period. As a result, the system will never reach steady state: the trajectory of the system will be a succession of transients relaxing on different time scales along different dimensions (Brezina et al. 2003a, b). This transient partial uncoupling of the different dimensions will alleviate the functional constraints that we have analyzed. By the same token, however, there will also be a pronounced history dependency (Brezina et al. 2003a, b) that will make the output of the system less predictable and controllable. Before proceeding to the real firing patterns, in Fig. 9 we illustrate these phenomena with the regular firing patterns.

Figure 9*A* shows how the key components of the system evolve from rest toward the steady state upon repetition of one representative firing pattern (indicated by the records of *f*_{B15} and *f*_{B16} at the *bottom*). The colored records show the modulator effects: the enhancement of Ca current (henceforth symbolized **C**), the activation of K current (**K**), the resultant of **C** and **K**, the net modulation of contraction size (**S**), and the increase of relaxation rate (**R**). In the model, based on measurements of the intrinsic dynamics of the effects in real experiments (Brezina et al. 2003a), the effects **C** and **R** are very slow: their intrinsic time constants are of the order of 100 s or more. In contrast, **K** is fast, with an intrinsic time constant of the order of 1 s. As a result, **C** and **R** reflect the firing pattern integrated over many cycles, whereas **K** responds to it even within a single cycle. Thus in Fig. 9*A* we see a smooth trajectory of **C** and **R**, but a cycle-by-cycle ripple in **K**, and consequently also in **S**. Yet another component of the NMT dominates the dynamics during the initial start-up: a slow facilitation of the release of the modulators SCP and MM (Vilim et al. 1996b, 2000; Brezina et al. 2000a, 2003a). This delays the onset of the modulation, so that in the first one or two cycles in Fig. 9*A* we see, essentially, the unmodulated contraction waveform. The contraction then undergoes a transient potentiation, reflecting passage of the modulator concentrations through the low range of values where **C** dominates **S**. [The intrinsic dynamics of the basal NMT also have to be considered, but they are not responsible for the slow potentiation because in the model they are fast, reaching steady state within at most two or three cycles (see, e.g., Figs. 7*D* and 8*F* of Brezina et al. 2003b).] Finally the contraction settles down, with this firing pattern, to a net depression where, a high enough concentration of MM having been attained, **K** dominates **S**. The functional performance *M* reflects this progression, transiently rising above zero only where the contraction is neither too large, nor too small (*top*).

Toward the end of Fig. 9*A*, we interrupted the repetition of the firing pattern with a brief gap. Gaps such as this are a prominent feature of the real firing patterns (see below). Although brief, the gap has a profound impact on the modulatory state, contraction, and performance when the firing resumes. The gap allows the fast effect **K** to decay, and consequently the slow effect **C**, which does not decay, to transiently reveal itself in an enormous increase in **S** and hence contraction amplitude. The performance *M* transiently reappears.

In Fig. 9*B* we stopped the simulation in Fig. 9*A* at the three timepoints marked by the vertical lines—at the very beginning of the simulation, just before the gap, and just after it—and mapped *M* over the B15 parameters *d*_{intra,B15} × *f*_{intra,B15} and the B16 parameters *d*_{intra,B16} × *f*_{intra,B16} simultaneously through the point of maximal *M* in the complete 4-dimensional space, as in previous figures. Here, however, the question is no longer how the modulatory state combines with the contraction elicited at the same time by the same firing pattern, but rather how the modulatory state set up by one pattern combines with contractions elicited subsequently by that pattern and by all other patterns. In Fig. 9*A* we already saw that the modulatory state varies dynamically with the history of the system, and alters subsequent performance by the same pattern. Figure 9*B* now shows how the dynamically varying modulatory state impacts the whole set of patterns. It can reduce as well as improve performance. A sufficiently depressive modulatory state eliminates performance by all patterns (*maps 3* and *4* of Fig. 9*B*), even those that, with their own modulatory states, would give good performance (compare with *maps 5* and *6* of Fig. 8). A potentiating modulatory state, on the other hand (*maps 5* and *6* of Fig. 9*B*), imparts performance to patterns that otherwise would not be capable of it. All of this follows from the fact that the modulatory state outlasts the firing patterns that produced it and persists to combine with the contraction waveform elicited by a new firing pattern. The slow dynamics of the modulator effects thus transiently uncouple the modulatory state from the contraction waveform.

### Real, irregular motor neuron firing patterns

Horn et al. (2004) used chronically implanted electrodes to record electrical activity in the ARC muscles of intact *Aplysia* feeding spontaneously on seaweed (*Gracilaria*). A representative recording of an entire, approximately 2½-h-long meal, consisting of many hundreds of cycles, is shown in Fig. 10*A*, *record 1. Records 2–6* are successive expansions of the indicated portions of the recording. We will use this single internally consistent, high-quality dataset as input to our model throughout the rest of this paper.

In earlier experiments, Cropper et al. (1990a, b) found that such extracellular recordings from the ARC muscle typically contain just two main classes of events, corresponding to the spikes of the motor neurons B15 and B16. Cropper et al. ascertained which class was which by comparing the sizes and shapes of the events recorded during the meal with those elicited by direct stimulation of B15 and B16 after subsequent dissection of the animal. Because Horn et al. (2004) did not do this, it was not immediately clear how to allocate the events in Fig. 10*A* to B15 and B16. At the very end of the meal, however, there were several cycles very different from those occurring previously, with much longer cycle periods and longer bursts composed of high-frequency firing of just the smaller of the two units (*records 3* and *6* in Fig. 10*A*). These are all the hallmarks of cycles of rejection (see, e.g., Church and Lloyd 1994; Cropper et al. 1990a; Weiss et al. 1986) that may tend to occur at the end of the meal as the animal satiates, identifying the smaller unit—as Cropper et al. also concluded in their recordings—as B16. Importantly, these cycles showed that even when B16 fired at high frequency, the corresponding events in the recording did not exceed a certain moderate amplitude (dashed line under *records 5* and *6* in Fig. 10*A*). We therefore allocated, throughout the entire meal, all events smaller than this amplitude to B16 and those larger to B15. This gave us nonrejection bursts (e.g., gray rectangle in *record* *5*) in which the smaller events of B16 came at moderately high frequency throughout the burst, whereas the larger events of B15 (asterisks) came at lower frequency often toward the end of the burst, all very much as in the recordings of Cropper et al. (1990a, b).

From the times of these events we then computed the instantaneous firing frequencies of the motor neurons B15 and B16, *f*_{B15} and *f*_{B16}, shown for the entire 2½-h-long meal in Fig. 10*B*. The same two records are repeated and expanded at the *top* of Fig. 11 to illustrate how the entire meal was composed of the discrete bursts of B15 and B16 firing (*top right*), which were generally well separated by intervals of silence in both B15 and B16 and so could easily be identified (for precise criteria see methods). Each burst, together with the interburst interval, was defined as comprising one cycle. There were 749 cycles in the whole 2½-h-long meal.

Horn et al. (2004) has already documented the extreme irregularity of this and other recordings of electrical activity from the ARC muscle as well as essentially all other aspects of neuromuscular activity in *Aplysia* feeding behavior. The successive expansions in Fig. 10*A*, for example, show well the “fractal” nature of the 2½-h meal—the bursts of activity are grouped into bouts, and bouts of bouts, on multiple time scales—and great variability in burst shape from one cycle to the next (e.g., in *record 4*). Figure 11 summarizes the statistics of this cycle-level variability. As with the regular firing patterns, two basic parameters are the burst duration *d*_{intra} and the cycle period *P.* Within each burst, the precise timing of spikes is irregular, too (*record 5* of Fig. 10*A*; Fig. 11, *top right*), but averaging *f*_{B15} and *f*_{B16} over the duration of each burst yields quantities comparable to *f*_{intra,B15} and *f*_{intra,B16}. For each cycle in the 2½-h meal, aligned with the records of *f*_{B15} and *f*_{B16} at the *top*, these four parameters are plotted down the *left-hand side* of Fig. 11. To the *immediate right* of each of these time series is a marginal histogram (scaled to approximate the probability density function) of the distribution of the values. Because the distributions are skewed, they are summarized using the median and quartiles (horizontal lines). Finally, on the *far right* are the corresponding distributions of cycle-to-cycle differences. These are pairwise differences between the values in successive cycles, a measure that is insensitive to any global progressive trends and emphasizes just the local, cycle-to-cycle variability (see Horn et al. 2004). The differences have been normalized by the mean of all cycles, so that −1 on the horizontal axis of these plots represents a decrease, and 1 an increase, equal in magnitude to the mean. In addition to the median and quartiles (vertical lines), the standard deviation of each distribution, σ (with the normalization, comparable to the coefficient of variation), is given.

The challenge now is to understand how, with these highly irregular, variable firing patterns, the B15/B16-ARC NMT can nevertheless continue to operate successfully.

### Irregular firing patterns produce irregular contractions

Figure 12 shows the trajectories of the key components of the modulated NMT during a simulation of the entire 2½-h-long meal. Because the inputs *f*_{B15} and *f*_{B16} are irregular (*top*), so are all the downstream components: the modulator concentrations *C*_{SCP} and *C*_{MM}, the modulator effects **C, K, S**, and **R**, and the output modulated contraction waveform *c* (*bottom*). However, the components react to the input irregularity on different time scales and so integrate or filter it to different degrees. Thus, as already explained, **C** and **R** are slow, whereas **K** is fast. **S**, the resultant of **C** and **K**, is consequently both slow and fast, exhibiting fast transients riding on a slow envelope.

The modulated contraction waveform is repeated in Fig. 13 (*top left*, in red) together with the statistics of its cycle-level variability presented similarly to those of the firing patterns in Fig. 11. (The actual contraction shapes in representative individual cycles may be seen in Figs. 16*D* and 19*E*.) The two basic parameters are the peak contraction and the contraction area (or, when divided by *P*, the mean contraction 〈*c*〉) in each cycle. In red, Fig. 13 shows these parameters for the modulated contraction waveform. The contraction waveform is no less variable, and in its area, for instance, may be even more variable, than the motor neuron firing patterns.

In blue, Fig. 13 also shows the contraction waveform and its statistics when the motor neurons were prevented from releasing their endogenous SCP and MM. This immediately allows us to settle one question of Horn et al. (2004). Horn et al. likewise observed highly variable contractions of the ARC muscle—in that case, contractions of the real muscle during feeding motor programs triggered by nerve stimulation in a reduced preparation—and found that their variability decreased when the cycle period of the motor programs was regularized. Horn et al. therefore suggested that the irregular cycling of the CPG creates the variability through mutual history-dependent reverberations of slow, cycle-spanning processes. Numerous such processes very likely exist in the CPG itself, but good candidates are already known in the periphery: the slow, cycle-spanning modulatory processes of the NMT. Figure 13 shows, however, that the peripheral modulation does not in any way amplify the basal variability of the NMT. Although the modulated contractions are generally larger than the unmodulated contractions, in their distributions of normalized cycle-to-cycle differences—the measure of variability that Horn et al. focused on and found to be regularized—they are no more variable. Thus, the standard deviations of the unmodulated and modulated distributions in Fig. 13, *right*, are very similar. For a more rigorous demonstration, we compared the entire distributions using the Kolmogorov–Smirnov test, a test that is very sensitive to cumulative differences in the shapes of distributions (e.g., Hoel 1971; Press et al. 2002). The Kolmogorov–Smirnov test found no significant difference (*P* > 0.8) between the unmodulated and modulated distributions of cycle-to-cycle differences of either the contraction peak or area.

### Modulation improves performance even with irregular firing patterns

When we now run the irregular contraction waveform that is elicited by the irregular firing patterns through our functional task, what performance does it achieve? A priori, it is by no means certain that it will necessarily achieve any very good performance.

With these irregular contractions, we set the threshold amplitude of the task, initially, to c̄ = 0.132 (as indicated in Fig. 13, *top right*). This is the mean value over all cycles of the meal of the mean contraction amplitude 〈*c*〉 in each cycle, for the unmodulated contractions. By setting equal to this value, according to our analysis in Figs. 3*C*, 5*C*, and 7*C*, we give the unmodulated contractions the benefit of the doubt—we make it possible for them to achieve their best performance. If modulation then still improves performance substantially, we can take that to be globally significant despite not being able to vary the irregular firing patterns arbitrarily to examine the entire space of patterns.^{3}

Figure 14*A* plots a time series of the performance *M* achieved by the modulated contraction waveform in each of the 749 cycles of the 2½-h meal. Some cycles do indeed fail to give any performance at all (*M* = 0). However, many cycles succeed to lesser or greater extent (*M* > 0). The percentage of cycles with *M* > 0 is 71.1%, and the mean *M* over all cycles, including those with *M* = 0, is 0.092. These two measures will be our standard measures of comparison among various conditions, in Figs. 15–18.

Figure 14, *C* and *D*, shows maps of *M* comparable to those constructed in previous figures for the regular firing patterns. In Fig. 14*C* we have projected all 749 cycles onto the plane spanned by the duty cycle *D* ≡ *d*_{intra}/*P* on the *x*-axis, pooling in this way the cycles with very different combinations of *d*_{intra} and *P*, and *f*_{intra,B15} (*top row* of maps) or *f*_{intraB16} (*bottom row* of maps) on the *y*-axis as before. The *left* pair of maps shows the unmodulated performance, the *right* pair the modulated performance. Clearly, the modulation improves the performance dramatically.

Does the modulation change—shrink, enlarge, or shift—the region of the space where performance is achieved? When all cycles with *M* > 0 are considered (Fig. 14*C*), the answer appears to be no. Individual cycles can be reduced as well as improved in performance (see Fig. 16*A* below), but there appears to be no particular part of the space where only one or the other happens. Therefore the region of *M* > 0 retains roughly the same bounds while *M* increases within it. When, however, only a subset of the cycles with superior performance is considered (e.g., in Fig. 14*D*, those with *M* > 0.15), then the modulation can be said to enlarge the region where this performance is achieved.

Figure 14*B* plots the same set of values of *M* against the cycle period *P.* Modulation improves performance at all practically significant values of *P*, but, as with the regular firing patterns (compare, e.g., Fig. 8), it has the largest effect at the smallest values of *P.*

### Components of the NMT contribute to performance in nonlinear, context- and task-dependent combinations

We ran the same 2½-h-long irregular motor neuron firing patterns *f*_{B15} and *f*_{B16} through the NMT many times while enabling, disabling, or altering various components. Some of the results are presented in Fig. 15.

With = 0.132 (Fig. 15*A*), we have already described the basic result, how the basal unmodulated performance (“No modulation”) improves when the full modulation is enabled (“Endogenous modulators”). The percentage of cycles with *M* > 0 improves from 47.5 to 71.1% (Fig. 15*A*, *top*) and the mean *M* improves from 0.036 to 0.092, by 156% (*bottom*, black circles). This improvement is statistically highly significant. [For clarity, we have not indicated statistical significance in Fig. 15 itself, but have given it in the figure legend. Because of the large number of cycles, standard errors are very small in Fig. 15, *bottom*—often smaller than the mean symbol size—and most differences between means, except those few where the means are virtually identical, are statistically significant.]

Enabling subsets of the modulation has interesting results. Enabling only the enhancement of Ca current (“**C** only”) does not improve the basal, unmodulated performance, but instead reduces it considerably. So does the activation of K current (“**K** only”). But with both of these effects—that is, with the complete modulation of contraction size **S**—performance is not significantly worse than the unmodulated performance (“**C** + **K** = **S**”, black circle), much better than the linear sum of the individual effects **C** and **K** would predict (*gray circle*). (The linear predictions in Fig. 15 are the sums of differences from the basal, unmodulated performance, but it is easy to see that multiplying the ratios, for example, would yield similarly inaccurate predictions.) The only effect that, alone, improves performance substantially is the increase of relaxation rate **R**, although even this then combines in a nonlinear way with **S** to give significantly better performance yet by the complete, fully modulated NMT.

Along the lines of our investigation with the regular firing patterns in Figs. 2–8, we also ran simulations in which we allocated all the spikes detected in Fig. 10*A* to B15, or all to B16, in each case with that motor neuron's modulation only. As Fig. 15*A* shows, allocating all the spikes to B15 gives relatively poor performance (“All B15”). Allocating all the spikes to B16, on the other hand (“All B16”), gives performance that is even better than that normally achieved by the complete, fully modulated NMT. Essentially, the large amount of SCP released in the first case, unopposed by MM, makes the contractions much too large relative to = 0.132, whereas in the second case the release of MM sets about the right contraction level, as in Fig. 7, *B* and *C*.

It might be thought at this point that we have a reasonably clear picture of the relative significance of the various components of the NMT. **R** appears to be by far the most important component of the modulation, in agreement with previous ideas (e.g., Brezina et al. 2000c, 2003b; Weiss et al. 1993). Furthermore, it seems that the system would perform equally well, perhaps even better, with just a single motor neuron, B16. However, this is not so. If we alter the particular quantitative requirements of the task, merely by altering *c̄*, the relationships among the components change dramatically. In Fig. 15*B*, for example, we set = 0.5. Now **C** alone somewhat improves performance, whereas with = 0.132 it reduced performance, and **C** + **K** = **S** improves the performance still more. **R** alone, on the other hand, becomes much less significant than it was with = 0.132. Indeed, now all three individual components have to be present simultaneously for a truly robust improvement in performance. Finally, B16 alone is now much less successful, and it is B15 that alone could do a better job even than the combined system.

Previous work has suggested that the different components of the modulation are likely to be important to different degrees in different types of cycles (Brezina et al. 2003b). The previous work, with regular firing patterns, studied idealized, sharply differentiated “biting,” “swallowing,” and “rejection” cycles. With the real, irregular patterns, the cycles are not so sharply differentiated. We can nevertheless ask the same question: in what kinds of cycles, with what parameter values, do the different components of the modulation make particular contributions to performance? As a first approach to the problem, we simply subtracted the performance *M* in each of the 749 cycles of a meal simulated with one set of modulatory components from *M* in the corresponding cycles of a meal simulated with another set of modulatory components. Selected examples are shown in Fig. 16.

Figure 16*A* shows the differences in *M*, Δ*M*, plotted as a function of the duty cycle, *D*, between the two basic meals: with and without modulation. The black circles are the 749 individual cycles. Cycles with Δ*M* > 0, above the horizontal midline of the plot, are cycles in which modulation improved performance; cycles with Δ*M* < 0, below the midline, cycles in which modulation reduced performance. Neither group shows any clear dependency in its structure on *D.* The red circles are the means of each group, and their horizontal locations are not significantly different. Thus the full modulation improves performance, or occasionally reduces performance, in cycles of all kinds, at least insofar as characterized by their duty cycle *D.*

Figure 16*B*, in contrast, compares the meal with the full modulation and a meal with **R** only. The full modulation tends to perform significantly better in cycles of smaller *D*, whereas **R** only performs better in cycles of larger *D.* The difference between the two meals is the presence of the modulation of contraction size **S**. The effect of **S**—in the context of the background presence of **R**—is illustrated in the representative contraction waveforms *a* and *b* in Fig. 16*D*, from the locations *a* and *b* in Fig. 16*B*. With small *D*, when the contraction occupies a small fraction of the cycle (waveform *a*), the contraction is too small without **S** (blue record), and about the right size, relative to = 0.132, with **S** (red record). On the other hand, with large *D*, when the contraction occupies much of the cycle (waveform *b*), the contraction is too large with **S** (red record) and about the right size without **S** (blue record).

The same mechanism in even more elementary form can be seen at work in Fig. 16*C*, and contraction waveforms *c* and *d* in Fig. 16*D*, where meals with **C** only and **K** only are compared. With small *D*, **C** only performs better than **K** only, which makes the contractions too small; with large *D*, on the other hand, **K** only performs better than **C** only, which makes the contractions too large. However, as we saw in Fig. 15*A*, both **C** only and **K** only actually make the performance worse, on average, than it would be if neither component were present. Only when both **C** and **K**, and furthermore also **R**, are present does the performance truly improve.

Altogether, from the results in this section, it appears that superior performance emerges in a nonlinear fashion through the complementary combination of all the different components of the modulated NMT. The relative contributions of the components vary in different contexts, in different tasks and even from cycle to cycle. In a specific context, one or another component may actually reduce performance: for that particular context, a better-performing, simpler NMT could be devised. Across all contexts and tasks that may arise, however, it is the complete NMT, with all the modulatory components, that guarantees the best overall performance.

### Dynamics of the modulatory state and control of the system

So far, with the irregular firing patterns, we have neglected dynamics. However, the fast dynamics of the modulator effect **K**, and even more the slow dynamics of the effects **C** and **R** that integrate over many cycles, can be expected to affect contractions and performance in the same way as they did with the regular firing patterns in Fig. 9.

Indeed, the dynamics of **C, K**, and **R** are well expressed during the 2½-h irregular meal. Figure 17*B* shows the trajectory of the modulatory state—essentially, the records of **S** and **R** from Fig. 12 plotted against each other—in the (**S**, **R**) plane during the meal. Soon after the beginning of the meal, the trajectory migrates to the general region of moderately elevated **R** and **C**, and so **S**. From this basic location, the trajectory then undergoes rapid vertical oscillations of **K**, and so **S**, in each cycle. These oscillations reach well out of the steady-state region, the gray region traced out by the steady-state effects (**S**_{∞}, **R**_{∞}) of all possible combinations of the SCP and MM concentrations *C*_{SCP} and *C*_{MM} (Brezina et al. 1996, 2003a). As previously analyzed theoretically (Brezina et al. 2003a), such excursions out of the steady-state region immediately indicate that **K**, and indeed **C** and **R**, are not in the steady state, with respect to *C*_{SCP} and *C*_{MM}, during much of the 2½-h irregular meal. Rather, the trajectory of the modulatory state is a transient through most of the meal, opening the possibility of significant history-dependent phenomena such as those in Fig. 9 (see Brezina et al. 2003a).

To explicitly examine how the dynamics of the modulatory state affect performance, we used two approaches. First, building on the idea that in each cycle the basal contraction is acted upon by a conceptually separable modulatory state, we scrambled the set of modulatory states with respect to the set of basal contractions. Essentially, we collected the 749 modulatory states (**C**, **K**, **R**) from the cycles of the 2½-h meal and randomly reassigned them, or their components, to the 749 cycles. The results are shown in Fig. 17*A*. Scrambling the whole modulatory state (**C**, **K**, **R**) (“Scrambled modulation”) impairs performance significantly, compared to the normal, unscrambled, fully modulated level (“Actual modulation”). The effect is fully reproduced by scrambling **K** only, whereas scrambling (**C**, **R**) has no effect. The latter is expected because the slow effects **C** and **R** do not vary much over the 2½-h meal (Fig. 12), so that scrambling them merely replaces the actual modulatory state with another of similar value. The fast effect **K**, in contrast, does vary from cycle to cycle (Figs. 12, 17*B*), and the fact that performance drops when it is scrambled shows that the cycle-by-cycle matching of **K** to the firing pattern and basal contraction is important for performance.^{4}

In a second approach, we altered the speed of the modulator effects. The results continue in Fig. 17*A*. Slowing **K** down to about the speed of **C** and **R** (“Slow **K**”) impairs performance just as scrambling **K** did. Again, the cycle-by-cycle adjustment of **K** appears beneficial. To understand it we can consider, as we did in analyzing Figs. 3 and 5, whether the fast or slow effect stabilizes or destabilizes the contraction mean 〈*c*〉 with respect to the goal, , around which performance is best. A fast **K** is stabilizing: high-frequency or prolonged motor neuron firing in any cycle, producing an excessively large basal contraction, simultaneously increases **K** to depress the contraction; conversely, low-frequency or short motor neuron firing, producing a small basal contraction, negligibly increases **K** and allows the contraction to be potentiated by the underlying persistent effect **C** (Figs. 9 and 12). By the same analysis, a fast **C** would be destabilizing, and a slow **C**, on the contrary, is stabilizing: the persistence of **C** integrated over many cycles ensures that small contractions are potentiated more, and large contractions less, than they would be potentiated within each cycle by a fast **C**. Hence we expect the normal, slow **C** to be beneficial: speeding it up should impair performance. Indeed, this is what we find (“Fast **C**” in Fig. 17*A*). Finally, by the same analysis, a fast **R** should be stabilizing, and indeed, speeding **R** up to about the speed of **K** (“Fast **R**”) improves performance even above that normally achieved by the system. This analysis fails to clarify, therefore, why **R** is normally slow (see discussion).

With these dynamics on multiple time scales interacting in the periphery and significantly influencing performance, to what extent is the behavior of various components of the system, and ultimately performance, controllable, or even predictable? To provide data with which to address this question, in Fig. 17*C* we have constructed a partial diagram of the average correlation strengths between selected components of the system over the 749 cycles of the fully modulated 2½-h meal. A representative correlation plot, from the composite parameter *f*_{intra,B16}*D* to **K**, is shown in the *inset* at *right.* (The inclusion of *D* tended to increase correlation strengths somewhat above those with *f*_{intra,B16} alone.) The coefficient of determination, *R*^{2}, of the best linear fit in this case is 0.50. The main diagram summarizes 25 such values, ranging from <0.01 (no correlation; there was no obvious nonlinear relationship, either, in these cases) to >0.99 (perfect correlation), indicated by the numbers themselves and by the thickness of the corresponding arrows.

The implications of the structure of this diagram for prediction and control will be addressed in the discussion. Here we note only how the diagram bears out some of the features that we have been concerned with so far. In cycle *n*, parameters of the CPG (here taken to include the motor neurons for convenience) such as *f*_{intra,B16}*D* predict reasonably well the current modulator effects **K**, and so **S**, and through them, as well as directly, the contraction mean 〈*c*〉—the fast components of the periphery. 〈*c*〉 in turn predicts, albeit less well, the performance *M.* The parameters of the CPG do not at all predict the current slow modulator effects **C** and **R**. On the other hand, **C** and **R** in the current cycle *n* are almost perfectly predicted by **C** and **R** in the preceding cycle *n* − 1. What is most striking in Fig. 17*C*, however, is how these reasonably strong correlations in the periphery contrast with the almost complete absence of correlations within the CPG itself. In the CPG, *f*_{intra,B15} does not predict *f*_{intra,B16}; *P* does not predict *f*_{intra,B15} or *f*_{intra,B16}; *P* or *f*_{intra,B16} in the preceding cycle does not predict *P* or *f*_{intra,B16} in the current cycle. Thus, as already concluded by Horn et al. (2004), it is the CPG that, by generating a new, essentially random combination of parameters of the motor program for each new cycle of behavior, is the primary source of the irregularity or variability throughout the system.

### With an irregular task, irregular firing patterns are optimal

Why does the CPG generate, almost actively it would seem, very variable motor neuron firing patterns, muscle contractions, and ultimately feeding movements? Horn et al. (2004) suggested that the variability implements a strategy of trial and error that the feeding animal pursues in the face of unknown functional requirements. Essentially, the animal does not know precisely what the task will be in the next cycle; perhaps it cannot know because the task is itself inconstant and irregular (see discussion).

Accepting for the moment that the variable strategy is a suitable one under the circumstances, surely it comes at a price? If the task *were* constant and could be determined, surely there is a regular motor neuron firing pattern that could achieve better performance than does the irregular pattern, which sometimes succeeds but often fails. If in some cycles the irregular pattern produces a contraction that is too large relative to *c̄*, and in other cycles a contraction that is too small, surely some average, intermediate firing pattern, repeated regularly, will produce an intermediate contraction—presumably, as we saw in Figs. 3*C*, 5*C*, and 7*C*, a contraction with mean amplitude 〈*c*〉 ≈ —that will succeed in every cycle.

To address this question, in Fig. 18*A* we compared the performance in the 2½-h irregular meal with that of several regular patterns. First, we constructed a regular pattern from the median parameter values of the irregular meal taken from Fig. 11. Remarkably, this average regular pattern achieves no performance whatever (“Median of irregular, regular”). This can be understood as a form of *pattern dependency* (Brezina et al. 1997): because of their nonlinear dependency on firing frequency (Brezina et al. 2003a, b), the basal contraction and the release of the modulators both respond disproportionately to the peaks in the irregular firing pattern, over groups of cycles as well as within cycles (see further below), and so are reduced considerably when these peaks are averaged out. (Downstream components of the system contribute additional pattern dependency: see, e.g., text footnote 4.) Next, we took the best-performing single cycle from the irregular meal, where it had *M* = 0.447. However, when taken out of its context and repeated over and over, this pattern, too, achieves no performance at all (“Best of irregular, regular”). Finally, we explicitly searched the 4-dimensional space of the motor neuron B15 and B16 parameters *d*_{intra,B15} × *f*_{intra,B15} × *d*_{intra,B16} × *f*_{intra,B16}, as we did in Figs. 6–9, to identify the best-performing regular pattern. This pattern (“Optimal, regular”)^{5} indeed gives better performance than that achieved in the 2½-h irregular meal. However, although the two situations are not completely comparable,^{6} it appears that the regular pattern is only marginally better.

Moreover, this regular pattern is only better in this particular task, with = 0.132. If the true requirements of the task have not been identified correctly, and the pattern therefore does not match the true task, its performance will diminish considerably. What happens, now, if the task is itself irregular, with requirements that vary from cycle to cycle? To model this, we imposed in each cycle a different value of randomly selected from a Gaussian distribution with mean and standard deviation both equal to 0.132. [This is about the degree of variability found across the various parameters of *Aplysia* feeding by Horn et al. (2004); see also Figs. 11 and 13.] The results are shown in Fig. 18*B*. Remarkably, with the irregularly varying *c̄*, the performance of the 2½-h irregular meal is not very different from its performance with the fixed *c̄.* The unmodulated NMT still gives some basal performance (“Irregular, unmodulated”), which the modulation then improves as much as before (“Irregular, modulated”). Both the unmodulated and modulated values of *M* are only about 19% smaller in Fig. 18*B* with the irregularly varying than with the fixed in Figs. 15*A* and 18*A*. The distribution of the performance over the space of the motor neuron firing parameters is not very different (Fig. 18*C*; compare with Fig. 14*C*). In contrast to this good performance of the 2½-h irregular meal, regular repetition of its median pattern again gives very poor performance (“Median of irregular, modulated” in Fig. 18*B*). Finally, as in Fig. 18*A* we explicitly searched the entire space *d*_{intra,B15} × *f*_{intra,B15} × *d*_{intra,B16} × *f*_{intra,B16} to identify the best-performing regular pattern. This pattern is essentially identical (parameters given in Fig. 18 legend) to the best-performing regular pattern with the fixed = 0.132. But now, faced with the irregularly varying *c̄*, the performance of this pattern (“Optimal, modulated”) is 47% lower, significantly worse than that of the 2½-h irregular meal.

Because the two situations are not completely comparable, the absolutely worse performance of the regular as compared to the irregular pattern should not be overinterpreted. However, the trends are significant. The performance of the regular pattern degrades considerably when faced with the irregularly varying *c̄*, whereas that of the irregular pattern is much more robust. To provide insight into how this works, Fig. 19 shows some details of the individual cycles of the three key meals from Fig. 18, that is, the meal with the “optimal” regular pattern but irregular (Fig. 19*A*), the 2½-h irregular meal with the fixed = 0.132 (*B*), and the 2½-h irregular meal with the irregular (*C*). In each case we have plotted the distribution of the cycles (*top*) and the performance *M* (*bottom*) as a function of the ratio of the mean contraction 〈*c*〉 to in each cycle. Fig. 19*D* shows one cycle of the “optimal” regular contraction waveform, and Fig. 19*E* shows nine representative cycles of the irregular contraction waveform of the 2½-h meal. Altogether, examination of Fig. 19 reveals three interrelated factors that, we believe, account for the superiority of the irregular firing pattern and contraction in the irregular task.

First, there is only a single shape of the regular contraction waveform, but many shapes, even with the same mean amplitude 〈*c*〉, of the irregular waveform. [The former can be seen in the unique performance, once the steady state has been reached, at each 〈*c*〉/ in Fig. 19*A*; the latter, in the scatter of performance at each 〈*c*〉/ in Fig. 19*B*.] The shape of the regular contraction waveform is thus necessarily a compromise, not only across cycles but even within each cycle. For example, the rate of rise of the contraction cannot be optimal simply in itself, nor can the rate of relaxation, but each of these features, which are strongly coupled to each other in this kind of task (Brezina and Weiss 2000; Brezina et al. 2000b, c), must be optimized also with regard to the other. In contrast, the irregular contraction waveform is not so restricted. In some cycles it can have a shape that combines a contraction rise that is optimal in itself with a relaxation that, too, is optimal in itself.

This does not yet guarantee better performance overall because cycles that combine well-performing contraction features might be expected to be balanced, on average, by cycles that combine poorly performing features. It appears in Fig. 19, *D* and *E*, however, that the shape of the irregular contraction waveform departs from that of the regular waveform systematically in one direction, a direction that gives better performance. Compared to the smoothly rounded regular contraction, the irregular contraction often has a distinctly square appearance: it rises much more steeply later in the cycle, giving better performance not only because it rises faster but also because it accomplishes this with a smaller mean amplitude 〈*c*〉. (Recall that our performance measure *M* takes 〈*c*〉 into account: the same movement accomplished with a smaller mean contraction amplitude—a smaller expenditure of energy—equates to higher performance.) As a result, whereas in Fig. 19*A* the performance of the regular contraction appears only around 〈*c*〉 ≈ as in our earlier analysis, the performance of the irregular contraction in Fig. 19, *B* and *C*, extends down to 〈*c*〉 where, in some cycles, it can be extremely good. All of this follows again from the fact that the motor neurons fire irregularly not only across cycles but even within each burst, with a tendency to fire more intensely as the burst progresses [see, e.g., Fig. 10*A, record 5*; the same was found by Cropper et al. (1990a)]. The irregular firing patterns thus outperform even the best regular pattern, in part, precisely because they are not members, and so are not bound by the limitations, of the canonical set of idealized, regular patterns.

Finally, because the distribution of the cycles along the dimension of 〈*c*〉/ is inherently very broad in the 2½-h irregular meal (Fig. 19*B*, *top*), it does not become much broader when is also made irregular (Fig. 19*C*, *top*). The inherently broad distribution may cost the irregular meal some performance at the outset (Fig. 18*A*) but it then stabilizes the performance against further degradation. The meal can always be expected to have a similar mixture of well- and poorly performing cycles, essentially irrespective of the distribution of (compare Fig. 19, *B* and *C*, *bottom*). Of course, the well- or poorly performing cycles will be different ones in each new meal—even when, as in the model, each meal replays the same exact sequence of motor neuron firing patterns. Thus, the six cycles marked by the colored circles in Fig. 19*E* can be seen to assume one grouping in the plot in Fig. 19*B*, but a very different grouping in Fig. 19*C*. Cycles that perform well in one meal—for example, *cycles 4* (pink) and *8* (purple) in Fig. 19*B*—may in the next meal be faced with an excessively high (*cycle 4*), or low (*cycle 8*), and perform poorly (Fig. 19, *C* and *E*). However, there will be a comparable number of cycles—*cycle 6* (green), for example—for which the converse is true. Altogether, although performance cannot be guaranteed in any individual cycle, it is almost sure to be robust overall even with completely unknown, randomly varying requirements of the functional task.

## DISCUSSION

The CPG that drives *Aplysia* consummatory feeding behaviors generates complex, multidimensional, variable patterns of activity. The ARC and other buccal-mass muscles that execute the behaviors incorporate complex, dynamically semiautonomous networks of peripheral modulatory actions that shape their contractions.^{7} Yet the CPG is connected to the muscles only through a narrow channel, in the case of the ARC muscle only through the trains of spikes of the two motor neurons B15 and B16. Under these circumstances, how can the central nervous system hope to control the peripheral musculature so as to produce successful behavior? The broad answer given by our work here is that it does not really control it, in a hierarchical, master–slave fashion. Rather, as in other systems (Chiel and Beer 1997), successful behavior “emerges from interactions of nervous system, body, and environment.” As Chiel and Beer (1997) have noted, in such cases “assessing the importance of changes in a motor program can only be done by ‘playing' it through the body, or a model of the body, and observing the behavioral consequences.” That is what we have done here. In emphasizing the interconnections between the nervous system, body, and the performance of a functional task vis-à-vis the environment, we view this work as a contribution to computational neuroethology (Beer 1990; Chiel and Beer 1997).

This work continues a long sequence of construction, analysis, and experimental testing of increasingly realistic models of the modulated B15/B16-ARC NMT (Brezina and Weiss 2000; Brezina et al. 1996, 2000b, c, 2003a, b; Weiss et al. 1993). Some of our questions here were raised before and answered conceptually; here these answers have been largely confirmed, now in a much more realistic setting. In other ways our emphasis here has been different from that in the previous work. Previously, for example, a major issue was how modulation maintains function over an extended range of behavior, as the behavior changes speed or switches between its biting, swallowing, and rejection subtypes (Brezina et al. 2000c, 2003b). In the recording of the 2½-h irregular meal on which we have based much of our work here, however, the subtypes of behavior cannot, after the fact, be unambiguously distinguished; the functional task that we have constructed is therefore appropriately generic. Instead, we have focused in this paper on how the system functions as an integrated whole under realistic conditions, with realistic release of the motor neurons' own endogenous modulators, and, even more important, with the realistic, highly variable motor neuron firing patterns recently documented by Horn et al. (2004). In answer to the four major questions asked in the introduction, we have found the following:

### Peripheral modulation is necessary for superior functional performance

When the motor neurons release their endogenous modulators in the fully functioning system, performance improves robustly, severalfold, under all circumstances that we have examined, with the idealized, regular motor neuron firing patterns (Figs. 6 and 8) as well as the real, irregular patterns (Figs. 14, 15, and 18). But is the modulation *necessary*? Previously we gave a general affirmative answer by showing, conceptually, that the extended range of behavioral tasks that the system will be required to perform is likely to be broader than the range of contraction shapes that an NMT with fixed properties can produce, so that, no matter how these properties are fixed, there will always be some tasks that cannot be performed at all: hence the properties of the NMT absolutely must be modulated (Brezina and Weiss 2000; Brezina et al. 2000c). Here we have not extended the task in the same way and furthermore have set its adjustable parameter so as to maximize the performance of the unmodulated contractions—so as to give the unmodulated system the benefit of the doubt. [The observed severalfold improvement of performance with modulation was thus simply the lower bound on the improvement: if the unmodulated contractions had not been initially favored, the modulation would have improved performance much more (see, e.g., Fig. 14*D*).] Giving the unmodulated system the benefit of the doubt has then allowed us to demonstrate that the contribution of the modulation is absolute in another sense: with the regular firing patterns, at least, it was possible to search the entire space of the patterns and show that it contained no pattern that, without modulation, could have given a better—or even nearly as good—performance as many patterns did with modulation (Figs. 6 and 8). The same may be true for the irregular firing patterns faced with an irregular task (Figs. 18, *B* and *C*, and 19). If the requirements of the task were made more stringent, moreover, performance of it might well require modulation even in the first absolute sense: without modulation, there would be no performance at all (Fig. 14*D*). In all of this, the modulation must be *peripheral.* It must interpose between and change the relationship of the motor neuron firing patterns to muscle contraction shapes—that is, the NMT. Modulating the input to the NMT, the motor neuron firing patterns themselves, is not a substitute.

### Components of the NMT interact in nonlinear, context- and task-dependent combinations for best overall performance

Repeatedly, we have found strong nonlinear coupling between components of the NMT, so that in combination they give very different performance than does each one alone, or even their linear sum. There is coupling between the basal NMT and the SCP “neuro-modulation” transform of motor neuron B15 (Figs. 2 and 3), between the basal NMT and the MM “neuro-modulation” transform of motor neuron B16 (Figs. 4 and 5), and then, overall, between the two motor neurons with their respective modulators (Figs. 6, 7, and 15). This last result complements previous work that, although not explicitly examining performance, showed by a formal argument how two motor neurons with different modulators provide a qualitatively broader coverage of the space of outputs of the system than does each one alone (Brezina et al. 1996). There is strong nonlinear coupling also between the modulator effects **C** and **K** to give **S**, and **S** and **R** to give the full modulation (Fig. 15). Although not analyzed here, finally, the different parameters of the motor neuron firing pattern such as *f*_{intra}, *d*_{intra}, and *P* converge nonlinearly through the NMT such that, in general, all parameters of the firing pattern influence all parameters of the contraction waveform and so functional performance (Brezina and Weiss 2000; Brezina et al. 2000b).

Because of this coupling, in part, the magnitude and even the sign of the contribution that each component makes to performance is heavily context dependent. It depends on the other components that are present, on the functional task, and on the motor neuron firing pattern in each particular cycle. In cycles with small duty cycle *D*, for example, the modulator effect **C** that makes contractions larger improves performance, whereas in cycles with large *D* it actually reduces performance, while **K** and **R**, which make contractions smaller, improve performance (Fig. 16). Thus **C**, or any other component, can improve performance in one cycle but reduce it—because the motor neuron firing parameters vary greatly from cycle to cycle (Fig. 11)—in the very next cycle. In a particular cycle, any particular component can thus actually be counterproductive. Indeed, even on balance over the entire 2½-h meal, **C**, when present alone, tends to reduce the performance (e.g., Fig. 15*A*). (Thus, the naı̈ve idea that larger contractions are necessarily better, and that modulation in a neuromuscular system works by making contractions larger, is not correct.) When present alone, **K**, likewise, reduces performance (Fig. 15, *A* and *B*). When both **C** and **K** are present together, however, performance is reduced much less or improved. Further improvement is achieved when **R**, too, is added to the mix. The relative contributions of these three components change dramatically when simply a quantitative parameter of the task, such as the threshold amplitude , is changed (compare Fig. 15, *A* and *B*). Yet the robust improvement in performance that is achieved by all three components in combination does not change—and this is the key point. We conclude that, even though any component of the NMT may be counterproductive in a specific context—for that specific context, a better-performing, simpler NMT could be devised—together the components constitute a complementary, mutually supporting network that best guarantees robust overall performance across all the contexts and tasks that may arise.

### Slow and fast dynamics of the modulatory state are important for performance

The mechanisms that maintain the robust performance of the peripheral network apparently include stabilizing, “homeostatic” mechanisms based on the dynamics of the modulator effects. As analyzed in results, if **C** were fast enough to respond to each cycle of the variable motor neuron firing pattern, it would amplify the variations and destabilize contraction amplitude; a slow **C**, on the other hand, integrates over the variations and stabilizes contraction amplitude. Indeed, **C** is normally slow, and speeding it up significantly impairs performance (Fig. 17*A*). The converse is true of **K**. A fast **K**, which counteracts the variations of motor neuron firing in each cycle, stabilizes contraction amplitude. Indeed, **K** is normally fast, and slowing it down, or scrambling it so that this component of the modulatory state no longer matches the contraction produced by the basal NMT in each cycle, impairs performance (Fig. 17*A*).

This analysis does not, however, explain the dynamics of **R**. A fast **R** should be stabilizing, and speeding **R** up improves performance (Fig. 17*A*). Yet **R** is normally slow. There are two possible explanations. First, because **R** is probably mediated by a modification of the contractile machinery of the muscle itself (Heierhorst et al. 1995; Probst et al. 1994), it is possible that it simply cannot be altered fast. Second, and more interesting, the dynamics of **R**, and those of **C** and **K** as well, may be understood better when related to the likely more complex statistics of the functional task (see below).

### Degrees of prediction and control by the central nervous system

If all effects in the periphery were fast enough to respond to the motor neuron firing pattern on a cycle-by-cycle basis, or, more generally, if only the steady state mattered (as in Figs. 2–8), there would be no difficulty in prediction and control. Simply from the current parameters of its own activity, the nervous system could predict the events in the periphery, and control them by appropriately setting its own parameters. The desired behavior would, in effect, be given by the solution of a set of instantaneous, algebraic equations, and provided these equations admitted of a simultaneous solution—which they may not (Figs. 4 and 5), but are likely to do in the complete balanced system (Figs. 6–8)—the nervous system could produce the behavior. However, because the slow dynamics in the periphery uncouple the events there from the current motor neuron firing pattern (Fig. 9), the nervous system, just like our model, faces instead a set of differential equations for whose solution it must take into account previous history. To what extent can it then control behavior?

Recall Fig. 17*C*, where we plotted the correlation strengths between components of the system over the course of the 2½-h meal. Although these are simply the elementary pairwise correlations, they give a fair idea of the situation. Interestingly, there *is* the possibility for relatively simple, albeit rough, prediction and control of the periphery by the nervous system. This is because the peripheral dynamics are either fast or *very* slow—there are none in the most troublesome, intermediate range. Accordingly, from the parameters of its own activity—in Fig. 17*C*, this is seen from the perspective of the representative CPG parameter *f*_{intra,B16}*D*—in the current cycle *n*, the nervous system can predict with reasonable accuracy the current value of the fast effect **K**. It can also predict, of course, the fast basal contraction. It cannot at all predict the current values of the slow effects **C** and **R** in this way. On the other hand, these values are exactly the same as they were in the preceding cycle *n* − 1; they reflect the firing of the motor neurons integrated over the past several minutes. So if the nervous system can keep track of its own overall activity—which slow neuronal processes can do—it will have a good representation of the current values of **C** and **R**. In combination all this information allows prediction of contraction parameters such as the mean contraction 〈*c*〉 in the current cycle, and of the performance *M.*

Because *M*, and even 〈*c*〉, is the product of convergence of many weak links, many pieces of information must be assembled, and the prediction is probably not very good. However, Fig. 17*C* also suggests that this does not matter. In striking contrast to the periphery, the CPG contains, apparently, no strong correlations at all. In each new cycle of behavior, the CPG generates a new, essentially random combination of parameters of the motor program. This is not simply a response to a randomly varying environment in the 2½-h meal: Horn et al. (2004) found the same randomly irregular CPG and motor neuron activity, both in intact animals and in reduced neuromuscular preparations without sensory feedback, even under controlled conditions with a constant stimulus. It thus appears that the CPG itself actively generates the irregularity or variability that is then seen throughout the system. (Figure 13 further supports this idea.) As we argue below, this variability is beneficial. Under these circumstances, although the network of events in the periphery is, as we have shown, fundamental for superior performance overall, it can be allowed to operate relatively autonomously: accurate prediction and control of the results in any particular cycle is not critical.

Indeed, the peripheral network can then be viewed as taking over some of the predictive and control functions of the nervous system. The nervous system, for example, no longer needs to integrate its own past activity in some neuronal variable to form an internal representation of what **C** and **R** are doing in the periphery. Rather, **C** and **R** are themselves a representation or memory of the past activity of the nervous system, maintained at the point where this memory can immediately be used to shape the contractions in the next cycle of behavior.

In sum, the nervous system exerts different degrees of control over different aspects of the peripheral function. Through the basal NMT, the nervous system directly controls the actual production of each feeding movement, its timing, and its type, in response to immediate sensory stimuli before and within each cycle (e.g., Evans and Cropper 1998 Kupfermann 1974; Rosen et al. 1982) and its own deterministic (Proekt et al. 2004; Zhurov et al. 2003) as well as stochastic (Horn et al. 2004) dynamics. The nervous system also exerts considerable control over the fast effect **K** to provide local adjustments of contraction shape. Beyond this, however, the nervous system should not be thought of as issuing direct commands but rather “suggestions” (Chiel and Beer 1997; Raibert and Hodgins 1993) for alteration of the trajectories of the semiautonomous dynamics of the periphery (e.g., Fig. 17*B*). These dynamics then perform computations of the feeding movements complementary to those performed by the nervous system.

### Variability as an optimal feeding strategy in an uncertain environment

Why does the CPG generate very variable motor neuron firing patterns, muscle contractions, and feeding movements? Horn et al. (2004) proposed that the variability is an integral part of the feeding animal's strategy for dealing with an uncertain environment. At any point in its meal, the animal may be confronted with any of a wide range of seaweed types and qualities that are best ingested with somewhat different feeding movements—in other words, confronted with a range of functional tasks. It is likely, furthermore, that the precise requirements of the task at any point are unknown. External tactile and chemical cues cannot fully distinguish, for example, the toughness of a piece of seaweed. How the seaweed can best be eaten, indeed whether it can be eaten at all, can be determined, by internal sensory feedback from the buccal mass and esophagus, only once the attempt has actually been made. In such circumstances, a trial-and-error strategy may be optimal. The CPG randomly generates a wide range of movements that efficiently sample the entire space of possibilities. Some of these movements will fail—and in a herbivore the penalty for failure in any particular cycle is relatively low, as indeed it was in our model here—but at least some will succeed. Successful movements can then be reinforced and stabilized in following cycles. [For further discussion of this idea and the evidence that already exists for it, see Horn et al. (2004).]

Is this strategy really “optimal”? To answer this question requires a quantitative cost–benefit analysis that, in elementary form at least, we carried out here. Clearly, given a fixed but unknown task, irregularly varying movements, which sometimes succeed, are better overall that regular movements that systematically miss the requirements of the task and never succeed. What if the fixed task could be known? Then regular movements that are specifically tailored to its particular requirements perform better, but not very much better than the irregular movements that pay no attention to the requirements of the task at all (Fig. 18*A*).

However, the proposal of Horn et al. (2004) is that the task, too, varies irregularly. Then the performance of the regular movements degrades dramatically, whereas—remarkably—that of the irregular movements is hardly affected (Fig. 18, *B* and *C*). In other words, the regular movements, fine-tuned to the requirements of one specific task, break down when faced with a changing mix of tasks, whereas the irregular movements continue to give essentially the same, robust performance no matter what may be the composition of the mix. There is, in that case, no need for the animal to identify the task requirements in any particular cycle precisely.

The robust performance is obtained because the irregular motor neuron firing patterns generate a redundant mix of irregular contraction shapes among which any mix of task requirements finds a sufficient number that give good and even excellent performance, whereas the single regular shape gives, at best, merely average performance (Fig. 19). Technically, the transformation of firing pattern to performance exhibits “positive” pattern dependency (Brezina et al. 1997) whereby a regular pattern constructed from the average parameters of the irregular patterns gives much worse performance (Fig. 18, *A* and *B*). Interestingly, in this regard it is not just the irregularity from cycle to cycle, but apparently even that within each burst of motor neuron firing that is significant (Fig. 19*E*). Previous work has always assumed, in the absence of explicit evidence to the contrary, that bursts could be simplified to their regular equivalents without changing the contractions of this “slow” muscle. This assumption must now be revised.

### Robust networks and trial-and-error strategies in other systems

Complex networks of interacting, mutually reinforcing elements are a common theme throughout biology. Recently studied examples include ecological webs (McCann et al. 1998), networks of ion currents that generate the voltage behavior of neurons (Goldman et al. 2001; Marder and Prinz 2002), and, in particular, networks of gene transcription, intracellular signaling and metabolism, and protein–protein interactions (e.g., Albert et al. 2000; Jeong et al. 2000, 2001; Li et al. 2004). These networks, too, are robust, rather than fine-tuned, in maintaining certain parameters of their output stable in the face of internal noise and external perturbations (Albert et al. 2000; Alon et al. 1999; Barkai and Leibler 1997). Usually such networks incorporate numerous feedback loops that stabilize a small number of discrete attractors in the space of trajectories or outputs of the system (e.g., Li et al. 2004). Our network here, however, does not include any especially strong feedback. Consequently, the network allows the variability generated by the CPG to appear in contractions, movements, and behavior as required by the proposed trial-and-error feeding strategy. What the network stabilizes is not performance in any one cycle, but the overall performance over many cycles in a variety of circumstances.

Strategies of random trial, error, and stabilization as solutions to the problem of uncertainty are also well known. Most concretely, many organisms perform random searches through actual, physical space. Unicellular organisms perform simple biased random walks (Adler 1975; Fenchel 2002; Webre et al. 2003); foraging ants and other insects (Bonabeau 1997; Camazine et al. 2001), and, constructed often in their imitation, robots (e.g., Mobus and Fisher 1999), execute more sophisticated searches that rely, nevertheless, on random variation to sample the space. Indeed, *Aplysia* itself appears to execute a random search to locate seaweed whose presence it has detected through nondirectional chemical cues (Teyke et al. 1992). More metaphorically, some of the genetic and biochemical networks already mentioned express and even amplify variability to sample the “phenotypic space” (Korobkova et al. 2004; McAdams and Arkin 1999). At higher levels of organization, “Darwinian” mechanisms of variation and stabilization are used to form synaptic connections (Hua and Smith 2004) and perhaps construct larger functional modules in the nervous system (Edelman 1993). The active generation of variability by CPGs has been discussed by Rabinovich and Abarbanel (1998) and Mpitsos (2000). A striking example that in some ways parallels our work here has recently been described by Varona et al. (2002) and Levi et al. (2004), who analyzed how a carnivorous mollusk, *Clione*, uses the variability actively generated by its CPG-like statocyst network to trace a quasi-random search path when hunting for its prey.

Other complex networks of modulator and neurotransmitter actions are only beginning to be analyzed from a global network perspective, but we suggest that some of the principles we have discussed here may apply also to other such networks, for example, the well-known networks of neuromodulation in the crustacean stomatogastric nervous system (Marder and Thirumalai 2002; Skiebe 2001), the vertebrate gut (Furness et al. 1992), and certain circuits of vertebrate as well as invertebrate brains (Agnati et al. 1995; Herbert 1993; Nässel 1996).

### Future directions

Here we have considered just two extremes of the functional task: a task that is either completely fixed, or completely random from cycle to cycle. However, it is likely that the real task will fall somewhere between these two extremes. It will be variable, with some possibility of even major change in any cycle, but at the same time with a higher-order temporal structure that is to some extent predictable. Figures 10 and 11, for example, already show such a structure in the cycle periods of the irregular meal: the cycles are grouped into bouts, presumably because multiple swallows are required to ingest each length of seaweed (see Horn et al. 2004), so that a short cycle period is likely to be followed by another short cycle period, although eventually a long cycle period will occur. Under these circumstances, the slow dynamics of the periphery can be understood as a memory of this higher-order temporal structure that then prepares the system for the next cycle, more like the last cycle the longer the bout of similar cycles has continued. Interestingly, the CPG itself has dynamics on a similar slow timescale that govern history-dependent transitions between the different types of feeding motor programs (Proekt et al. 2004). Preliminary work shows that these transitions are indeed reflected in muscle contractions and movement in the periphery (Zhurov et al. 2003). All this makes perfect sense from the brain–body perspective of computational neuroethology (Chiel and Beer 1997), which predicts that there should exist a matching between the dynamics of the nervous system, the dynamics of the periphery, and the temporal structure of the animal's environment. Our work here now makes it possible to study such matching in the *Aplysia* feeding system, both in the model of the modulated B15/B16-ARC NMT and in real, behaving animals.

## GRANTS

This work was funded by grants from the National Institutes of Health.

## Acknowledgments

We thank A. Proekt for helpful suggestions.

## Footnotes

↵1 Strictly, the complete model has no real steady state because it incorporates a slow, progressive decay of modulator release, even with a constant firing pattern, attributed to progressive “depletion” of the releasable modulator pool (Brezina and Weiss 2000). To allow the model to come to steady state, we disabled this depletion (see methods).

↵2 Indeed, by

*Eq. 6f*in methods, the contractions elicited by B15 and B16*are*identical, except that each value of*f*_{intra,B16}is equivalent to a value of*f*_{intra,B15}that is 4 Hz smaller.↵3 The functional task used with the irregular firing patterns and contractions differed from that used with the regular firing patterns in one important detail. Whereas the regular task based its threshold-transition windows a priori on the cycle period, the irregular task derived them a posteriori from the observed burst of motor neuron firing and the contraction waveform itself (see methods). This was necessary because the cycle period in the same functional sense as in the regular task (that is, as the time interval during which the sequence of contraction and relaxation of the muscles and the coordination between the closure and opening of the radula and its protraction and retraction must occur) was not directly observable: there were many very long cycle periods (Fig. 11) with short bursts of firing, best interpreted as consisting of a short “functional” cycle period followed by a gap or pause in the meal. Ultimately, the difference in the formulation of the regular and irregular tasks could be seen as making the coupling requirement between the closure and opening of the radula and its protraction and retraction, like the other requirements of the task, relatively rigid in the regular task but more flexible in the irregular task.

↵4 To enable the scrambling of modulatory states between cycles, the modulatory state first had to be fixed in each cycle, by setting the values of

**C**,**K**, and**R**throughout each cycle to the mean, period-averaged values 〈**C**〉, 〈**K**〉, and 〈**R**〉 computed over that cycle. Interestingly, this in itself reduced performance by about 6% (“Actual modulation” in Fig. 17*A*vs. “Endogenous modulators” in Fig. 15*A*). Presumably, it was again the fixing of**K**to 〈**K**〉 that was responsible. Thus, even the variation of**K**within each cycle has an impact on performance.↵5 Interestingly, here with = 0.132 the optimal regular pattern (parameters given in Fig. 18

*A*legend) was quite different—with much lower*f*_{intra,B15}—from that found for = 0.5 in Figs. 6–8. Thus a relatively small change in the task requirements apparently produced a major migration of the maximal point of the performance landscape over the space of parameters.↵6 Whereas the cycle periods of the irregular firing patterns of course vary, in the search for the best-performing regular pattern we set

*P*= 5 s. This presumably gave the best performance (cf. Fig. 8) yet was still realistic: it was approximately the lower bound on*P*over the 2½-h irregular meal (Fig. 11) as well as in previous observations of*Aplysia*feeding behavior (see Brezina et al. 2003a). In addition to this difference in*P*, the regular and irregular situations differ in the precise formulation of the functional task (text footnote 3).↵7 All other buccal-mass muscles studied so far have been found to be modulated much like the ARC, although the balance of the individual modulatory effects can differ (see, e.g., Church et al. 1993; Evans et al. 1999; Fox and Lloyd 1997; Hurwitz et al. 2000; Scott et al. 1997).

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2005 by the American Physiological Society