## Abstract

Many physiological systems are regulated by complex networks of modulatory actions. Here we use mathematical modeling and complementary experiments to study the dynamic behavior of such a network in the accessory radula closer (ARC) neuromuscular system of *Aplysia*. The ARC muscle participates in several types of rhythmic consummatory feeding behavior. The muscle's motor neurons release acetylcholine to produce basal contractions, but also modulatory peptide cotransmitters that, through multiple cellular effects, shape the contractions to meet behavioral demands. We construct a dynamic model of the modulatory network and examine its operation as the motor neurons fire in realistic patterns that change gradually over an hour-long meal and abruptly with switches between the different feeding behaviors. The modulatory effects have very disparate dynamical time scales. Some react to the motor neuron firing only over many cycles of the behavior, but one key effect is fast enough to respond to each individual cycle. Switches between the behaviors are therefore followed by rapid relaxations along some modulatory dimensions but not others. The trajectory of the modulatory state is a transient throughout the meal, ranging widely over regions of the modulatory space not accessible in the steady state. There is a pronounced history-dependency: the modulatory state associated with a cycle of a particular behavior depends on when that cycle occurs and what behaviors preceded it. On average, nevertheless, each behavior is associated with a different modulatory state. In the following companion study, we add a model of the neuromuscular transform to reconstruct and evaluate the actual modulated contraction shapes.

## INTRODUCTION

In recent years an enormous amount has been learned about how different transmitters and modulators regulate cellular functions. It is now apparent, however, that many physiological systems are regulated, not by single transmitters, but by multiple transmitters all present simultaneously in various combinations and affecting multiple aspects of the system's activity. It is then the overall, integrated result of this complex network of regulatory effects, rather than the effect of any one transmitter alone, that is functionally important (Brezina and Weiss 1997). By the same token, however, the overall complex effect is that much more difficult to measure, understand, and predict.

We are studying these questions in an advantageous model system, the accessory radula closer (ARC) neuromuscular system of *Aplysia*. This system contains numerous transmitters and modulators, and its simple construction from just a few large, identified cellular elements allows the various effects of the transmitters to be studied in detail but also understood in their functional relations to each other and to the behavior of the animal. Extensive studies over the past 25 years (reviewed by Kupfermann et al. 1997; Weiss et al. 1993) have characterized many aspects of the system to the point where such integrative work has become possible.

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). The muscle is innervated by 2 motor neurons, B15 and B16, which release their classical transmitter, acetylcholine (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, 1990a; Lloyd et al. 1984; Vilim et al. 1996a; 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 3 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. 1990a; 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 gives the final shape of the contraction (see Brezina et al. 2000a). The main elements of this modulatory network are graphically summarized in Fig. 1. Typical examples of contractions modulated by different concentrations of SCP and MM can be seen in Fig. 2*A*, “SCP” and “MM.”

In previous work we drew together the available data to construct a mathematical model, which we then verified experimentally, of how different combinations of SCPs and MMs interact through the network of modulatory actions to produce different ratios of the 2 ultimate, functionally relevant effects on contraction size and relaxation rate (Brezina et al. 1996; see also Brezina and Weiss 1997; Brezina et al. 2000a). This model showed how the space of SCP and MM concentrations maps to the space of the 2 functional effects (Fig. 2*B*), and specifically to a region within that space over which combinatorial changes in SCP and MM concentrations can achieve any combination of the effects on contraction size and relaxation rate, as indeed was observed experimentally (Fig. 2*A*, “SCP + MM”).

This work had several limitations, however. Perhaps most important, the model and its verification were both static: with steady applications of fixed concentrations of SCPs, MMs, or SCP–MM mixtures, the effects were modeled or measured only after they had developed to quasi-steady state. What is really wanted, however, is a dynamic model, which would predict the dynamic trajectory of the system through the spaces in Fig. 2*B* as the 2 motor neurons B15 and B16 fire in different behaviorally relevant patterns, releasing different SCP–MM combinations, that change gradually over the course of an hour-long meal or more abruptly with switches between the different types of feeding behavior (Kupfermann 1974; Susswein et al. 1976, 1978; Weiss et al. 1982). Construction of such a dynamic model entails modeling the dynamics at each stage, from the pattern of the motor neuron firing through the dynamics of the release of the modulators to the dynamics of their various effects. Fortunately, sufficient data on all of these aspects of the system are available in the literature, or from new experiments reported here, and in this study we construct a dynamic model of the modulation and explore some of its predictions. In the following companion study (Brezina et al. 2003; referred to as Paper II), we then join this model to a model of the basal B15/B16-ARC neuromuscular transform (see Brezina et al. 2000a,b) to reconstruct, and test experimentally, the actual contraction shapes.

## METHODS

### Experiments

Methods were as in previous work in the ARC neuromuscular system, briefly described in the following sections.

dissociated arc muscle fibers. Single fibers, dissociated from the ARC muscle by collagenase treatment, were impaled with a single sharp microelectrode and voltage-clamped using the discontinuous single-electrode technique (for details see Brezina et al. 1994c). The fiber was held by the electrode in a fast-flowing stream of solution; modulators were added to the solution stream (in Ca-current experiments) or puffed onto the fiber from an adjacent pressure-ejection pipette [in K-current experiments; see Brezina (1994), Brezina et al. (1994a)]. The latter technique gave very rapid changes in modulator concentration (see Fig. 7*A*, “ACh”). All experiments were done at room temperature.

To study the modulator-activated K current, the fiber was held at a steady voltage, usually –30 mV or in some experiments –40 or –50 mV, in normal artificial sea water (ASW), containing (in mM) 460 Na^{+}, 10 K^{+}, 11 Ca^{2+}, 55 Mg^{2+}, 602 Cl^{–}, and 10 HEPES buffer (pH 7.6). See further Brezina et al. (1994b,c, 1995).

To study the enhancement of Ca-channel current, we recorded Ba current through the Ca channels. Ba currents were elicited by 30- or 100-ms steps from –90 to 0 mV every 10 s in Ba/TEA/4-AP ASW, containing 460 mM tetraethylammonium (TEA^{+}) instead of Na^{+}, 11 mM Ba^{2+} instead of Ca^{2+}, and 10 mM 4-amino pyridine (4-AP). Leakage and capacitive currents that remained when Co/TEA/4-AP ASW was substituted at the end of the experiment were routinely subtracted. See further Brezina et al. (1994a,d, 1995).

intact arc muscle. To evaluate the kinetics of the increase of the muscle's relaxation rate (Figs. 9 and 10), we reanalyzed relevant experiments from several studies done in recent years in our laboratory, unpublished as well as those of Brezina et al. (1995) and Orekhova et al. (2003). These experiments all used a standard preparation and methods (see also Cohen et al. 1978; Cropper et al. 1987a; Vilim et al. 1994; Weiss et al. 1978). Briefly, the preparation consisted of the buccal ganglia, an ARC muscle, and the connecting nerve 3 through which the buccal motor neurons B15 and B16 innervate the muscle. The muscle was placed in a separate subchamber, continuously perfused throughout the experiment, to which the modulators were applied. Either motor neuron B15 or B16 was impaled with a double-barreled microelectrode and stimulated with brief current injections to fire spikes in the desired pattern, usually at 10–15 Hz (B15) or 15–25 Hz (B16) for 1.5 or 2 s every 30 s. These values were chosen so that each spike burst produced a muscle contraction of moderate size (see, e.g., Figs. 2*A*, 9*A2*). Such patterns of firing release ACh but not the motor neurons' endogenous peptide modulators (Cropper et al. 1990a; Vilim et al. 1996b). Contractions were monitored with an isotonic transducer. All experiments were done at room temperature.

### Modeling

Here we give the equations of the model, present the supporting data—in particular, unpublished data assembled or obtained in new experiments specifically for this purpose—and describe in detail how we used the data to evaluate parameters of the model equations. The principal variables, parameters, and symbols used, along with the parameter values obtained from the data, are summarized in Tables 1 and 2.

behaviorally realistic motor neuron firing patterns. We adopt here a simple, but reasonably realistic, picture of *Aplysia* consummatory feeding behavior, consisting of major elements on which there is consensus, as well as some quantitative data, in the literature. We assume that the behavior is cyclical, composed of many relatively short cycles of 3 distinct types— biting, swallowing, and rejection cycles— occurring sequentially, in principle in any order, in a long meal (see Kupfermann 1974; and other references below). In our model we also allow cycles of a fourth, null “behavior” to model cycle-length pauses and longer gaps of inactivity within the meal.

Corresponding to this behavior is a cyclical, meal-long firing pattern of motor neurons such as B15 and B16 (Cropper et al. 1990b; Morton and Chiel 1993a,b). We model this as a pattern of the “instantaneous” firing frequency *f* as a function of time *t*. Because, to a first approximation, the motor neurons fire only one spike burst of relatively constant intraburst frequency in each cycle, the pattern of *f* can be adequately described by just 3 parameters such as the burst duration *d*_{intra}, the interburst interval *d*_{inter}, and the intraburst firing frequency *f*_{intra} (see Brezina et al. 1997, 2000b). We also define the cycle period *P* = *d*_{intra} + *d*_{inter}. Bursting patterns of this simple kind (illustrated, e.g., in Fig. 3*B**, bottom*) can be compactly described by the equation (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 meal which began at *t* = 0. The major departure from previous formulations (Brezina et al. 1997, 2000b) is that here the parameters *P, d*_{intra}, *d*_{inter}, and *f*_{intra} are not fixed but vary depending on, most importantly, the type of the cycle and its location within the meal, that is, on the time *t*. Note also that *Eq. 1a* implies that the interburst interval precedes the burst within each cycle, so that *d*_{inter} is also the latency from *t** to the start of the burst in each cycle, in particular in the first cycle.

For the length of the entire meal, although shorter and longer times have been reported, *L* = 1 h is reasonable (Kupfermann and Carew 1974; Kupfermann and Weiss 1982; Kuslansky et al. 1987; Susswein et al. 1976).

We now have to specify the functions *P*(*t*), *d*_{intra}(*t*), *d*_{inter}(*t*), and *f*_{intra}(*t*). Two kinds of relevant data are available in the literature:

1. cycle period. A number of studies report various interrelated behavioral measurements, for example of the latency to a food-stimulated bite or of the interbite or interswallow interval (e.g., Hurwitz and Susswein 1992; Kupfermann 1974; Kupfermann and Weiss 1982; Rosen et al. 1982, 1989; Susswein et al. 1976, 1978; Weiss et al. 1986), that we take to be more or less direct reflections of the speed of cycling of the behavior, that is, of the cycle period *P*. When, furthermore, such measurements were made at various times during the meal, they provide a measure of *P*(*t*), and, in particular, of its changes during arousal and satiation. These measurements [from Susswein et al. (1976, 1978), Kupfermann and Weiss (1982), and Rosen et al. (1989)] are replotted in Fig. 3*A* and fitted with the function (1c) where (1d) describes the initial rapid decrease of *P*(*t*) in food-induced arousal, and (1e) describes the later slow increase of *P*(*t*) in satiation, with *t* in seconds. *Eqs. 1c-e* give a *P*(*t*) never smaller than 5.5 s and normally, in the middle part of the meal after arousal but before significant satiation, about 6 s, a very reasonable value for biting and swallowing [see the references above, and Cropper et al. (1990b)]. In rejection, however, it appears that a corresponding reasonable value is about 18 s (Weiss et al. 1986). In rejection cycles, therefore, we multiply *P*(*t*) given by *Eqs. 1c-e* by a factor of 3. Because the cycle period *P* is an attribute of the entire behavior, we use *P*(*t*) given by *Eqs. 1c-e*, or 3 times that value in rejection, for both motor neurons B15 and B16.

2. intra-cycle firing parameters. Essentially the only data available come from the work of Cropper et al. (1990b), who used electrodes chronically implanted in the ARC muscle to record the firing patterns of motor neurons B15 and B16 in representative biting, swallowing, and rejection cycles. These recordings give an idea of the reasonable values of *d*_{intra}, *d*_{inter}, and *f*_{intra} in the 3 types of cycles, but leave unclear how, if at all, these parameters change in relation to *P*(*t*) through the meal. As our best model of the data, we adopt fixed values of the burst parameters *d*_{intra} and *f*_{intra} throughout the meal: for biting, *d*_{intra,B15} = 1.3 s, *f*_{intra,B15} = 10 Hz, *d*_{intra,B16} = 2 s, and *f*_{intra,B16} = 18 Hz; for swallowing, *d*_{intra,B15} = 3 s, *f*_{intra,B15} = 10 Hz, *d*_{intra,B16} = 3.7 s, and *f*_{intra,B16} = 18 Hz; and for rejection, *d*_{intra,B16} = 5 s and *f*_{intra,B16} = 20 Hz (only B16 fires in rejection). In all cases, the burst ends at the end of the cycle (thus in biting and swallowing, the B16 burst starts before the B15 burst, but both end simultaneously). The interburst interval or latency *d*_{inter} is then left to vary according to the relation *d*_{inter}(*t*) = *P*(*t*) – *d*_{intra}.

Figure 3*B* illustrates how the model constructs from the raw parameter constants and functions the firing patterns of the motor neurons B15 and B16 in individual biting, swallowing, and rejection cycles. The type of cycle being commanded and its particular parameter values are determined at the beginning of the cycle (i.e., at *t* = *t**). The cycle is then performed in its entirety, essentially as a fixed-action pattern; a switch to a different cycle type takes effect only in the next cycle.

Null cycles within the meal are modeled simply as cycles of period *P*(*t*) in which neither motor neuron B15 nor B16 fires.

As regular equivalents of the variable firing patterns modeled here, the standard reference patterns *d*_{intra,B15} = 3.5 s, *d*_{inter,B15} = 3.5 s, and *f*_{intra,B15} = 10–12 Hz for motor neuron B15, and *d*_{intra,B16} = 3.5 s, *d*_{inter,B16} = 3.5 s, and *f*_{intra,B16} = 20 Hz for motor neuron B16, which are often used in physiological studies (e.g., Vilim et al. 1996a,b, 2000; see further below), seem reasonable, essentially averaging the parameter values of biting and swallowing cycles over the middle part of the meal. For analytical simplicity, we continued to use such regular patterns in certain simulations and experiments here (e.g., in Fig. 14, and Paper II).

The canonical set of behaviorally realistic firing-pattern parameters arrived at in this modeling is summarized in Table 2.

modulator release. A considerable amount of experimental data is available on the release of the modulators, the SCPs and the buccalins (BUCs), from motor neuron B15, from the work of Vilim et al. (1996a,b). Vilim et al. used radioimmunoassay to directly measure the amounts of SCP and BUC appearing in the perfusate of the ARC muscle while motor neuron B15 was stimulated to fire in various patterns; the outflow of the modulators from the muscle was taken as a reflection of their release from the neuron's terminals within the muscle. In a previous study (Brezina et al. 2000c), we have already used this data to construct a model of modulator release from B15, which we simply incorporate here. The basic equations are (2a) (2b) (2c) where *r* is the release; *S* is the size of the releasable pool of modulator; *p* is the probability of release or availability for release of the modulator in *S*; *k _{p}*

_{+}and

*k*

_{p}_{–}are the on and off rate constants of

*p*, respectively;

*y*is an exponent providing nonlinearity of release; and

*f*(

*t*) is the firing pattern given by

*Eq. 1a*. The factor , with

*Q*

_{10}= 5.3 × 10

^{–2}, is included here to model the dependency of the release on the temperature

*T*(see modulator concentrations below); unless specified otherwise,

*T*= 15°C and so . For release from B15, we found

*k*

_{p}_{+,B15}= 4.0 × 10

^{–10},

*k*

_{p}_{–,B15}= 3.4 × 10

^{–3}s

^{–1}, and

*y*= 3 [Model II of Brezina et al. (2000c)].

*S*

_{0,SCP}, the initial size of

*S*

_{SCP}before any release, was found to be 541 fmol, without any evidence for replenishment of

*S*

_{SCP}on the time scale (

*L*≤ 1 h) of the experiments.

For release of the MMs and BUCs from motor neuron B16, similar, albeit less extensive, data exist in the work of Vilim et al. (2000). The relevant data are reproduced in Fig. 4. Figure 4*A* shows, to scale, the standard reference firing pattern used by Vilim et al. for B16: note that it is a regular version of the pattern in *Eq. 1a* and Fig. 3*B*, with behaviorally relevant parameter values. With the pattern in Fig. 4*A*, Vilim et al. obtained the MM and BUC outflow in Fig. 4*B*, and with systematic variations of the parameters *d*_{inter}, *f*_{intra}, and *L*, the data in Fig. 4*C*.

To model these data, we used the same *Eqs. 2* as for B15, and fitted them to the data using a similar strategy, as described in detail by Brezina et al. (2000c) (their Model II). Because there is good evidence that in motor neuron B16 the MMs and BUCs (as in B15 the SCPs and BUCs) are contained in the same vesicles and obligatorily coreleased (Vilim et al. 2000), throughout we pooled MM and BUC data or constrained the fits to yield identical parameter values for both, except absolute amounts. Briefly, we first fitted *Eq. 17* of Brezina et al. (2000c) to the data in the interval 10 < *t* < 60 min in Fig. 4*B*, yielding a single-exponential rate constant of decay of 0.027 min^{–1} (time constant of about 37 min), and, extrapolating the fit to *t* → ∞ and integrating the area under the curve (the experimental curve for 0 < *t* < 10 min, the fitted curve later), *S*_{0,MM} = 2,690 fmol and *S*_{0,BUC} = 697 fmol (indicated in Fig. 4*C2*). Next, with these values of *S*_{0}, we fitted *Eq. 16* of Brezina et al. (2000c) to the data in Fig. 4*C1*, in a first pass (fit not shown) yielding as the best integral value, as for B15, *y* = 3. With this *y*, the rate constant from Fig. 4*B* yielded, by the argument in Brezina et al. (2000c), the “dissociation constant” *K _{p}* ≡

*k*

_{p}_{–}/

*k*

_{p}_{+}≈ 8.9 × 10

^{7}s

^{–1}. Finally, with these values of

*S*

_{0},

*y*, and

*K*, fitting

_{p}*Eq. 15*of Brezina et al. (2000c) simultaneously to the data in both Fig. 4,

*C1*and

*C2*yielded

*k*

_{p+,B16}= 7.4 × 10

^{–11}and

*k*

_{p}_{–,B16}= 6.6 × 10

^{–3}s

^{–1}.

For both motor neurons B15 and B16, according to *Eq. 2*, release depends on the firing pattern in 2 ways. Release is instantaneously gated by the firing pattern, but also reflects the probability of release *p*, which, with the small values of *k _{p}*

_{+}and

*k*

_{p}_{–}that we have found, responds to the firing pattern slowly, over many cycles. The slow rise in

*p*after the motor neurons first start to fire is responsible for the “facilitation” of release observed experimentally (Fig. 4; Brezina et al. 2000c; Karhunen et al. 2001; Vilim et al. 1996a, 2000). Finally, as release continues, it gradually decreases— even if the firing pattern and

*p*remain constant—as the releasable pool of modulator

*S*becomes depleted.

Comparing *S*_{0,SCP} and *S*_{0,MM}, we see that motor neuron B16 contains approximately 5 times more releasable MM than B15 does SCP, reflecting probably the ratios of synthesis (see Church and Lloyd 1991). The ratio of release of SCP and MM is not fixed, however, but depends, in particular, on the relative firing frequencies of B15 and B16. Because release from B16 requires somewhat higher frequency than release from B15 does [compare, e.g., Fig. 4*C1* here with Fig. 4*A* of Brezina et al. (2000c)], when the two motor neurons fire at the same frequency, B15 releases more SCP than B16 does MM. In behavior, however, B16 fires at higher frequency than B15 (see behaviorally realistic motor neuron firing patterns above). Then B16 releases considerably more MM than B15 does SCP [compare Fig. 4*B* here with Fig. 2*B* of Brezina et al. (2000c)].

modulator concentrations. To model the effects of the released modulators, which were all studied experimentally by applying known modulator *concentrations* to the ARC muscle, it was necessary to convert the released modulator amounts into the effective modulator concentrations, *C*. We assume that each modulator is released into some volume *v*, within the muscle, from which it exerts its effects. The simplest model for *C* then is (3) where *k _{C}* is the rate constant of removal of the modulator from

*v*. This formulation lumps together all removal processes, principally enzymatic degradation, diffusion, and removal by blood (hemolymph) flow through the muscle. Furthermore, postulating a single fixed volume, the same for all effects, at all concentrations, is likely to be an oversimplification. Nevertheless,

*Eq. 3*specifies well enough the 2 essential aspects, the general amplitude and the speed of modulator concentration changes.

To model *Eq. 3* with the correct values of *v* and *k _{C}*, we required an effect that had been measured with known concentrations of exogenous modulators but also known motor neuron firing patterns—and thus, using the release model, known amounts of released modulators—so that the two could be compared. Such data are available for elevation of cAMP in the ARC muscle by the modulators. Both the SCPs and the MMs elevate cAMP (Brezina et al. 1994a; Cropper et al. 1990a; Hooper et al. 1994a; Lloyd et al. 1984; Whim and Lloyd 1989, 1990) and activate cAMP-dependent protein kinase (PKA) (Hooper et al. 1994a,b) in the ARC muscle fibers. [Although the cAMP/PKA pathway mediates 2 of the downstream effects that we model below— on the Ca current and on the muscle's relaxation rate (Fig. 1)—we model those effects directly. The cAMP model in this section is used solely to establish the parameters of

*Eq. 3*.] The essential data on the elevation of cAMP are reproduced in Fig. 5.

Most measurements were of the total cAMP content of the ARC muscle, not in absolute terms but relative to the basal (unstimulated) content. We deal throughout with this normalized variable, *A*. It is probable that the SCPs and MMs activate adenylyl cyclase, rather than inhibit phosphodiesterase, and it appears that activation of adenylyl cyclase, at least by another ARC-muscle modulator, serotonin, is reasonably fast [see Fig. 3*B* of Weiss et al. (1979)]. A simple model for *A* is then (4a) where *K _{A}* is the dissociation constant of a binding reaction representing the lumped equilibrium states of all of the fast reactions up to the production of cAMP,

*k*

_{A}_{+}is the maximal rate of cAMP production, and

*k*

_{A}_{–}is the rate constant of cAMP breakdown.

A potentially serious problem in fitting *Eqs. 3* and *4a* to the data is that much of the cAMP data were obtained at room temperature, in many cases 22°C, whereas the work of Vilim et al. (1996a,b, 2000) on which we based the release model was done at 15°C. As both Vilim et al. (1996b) and Whim and Lloyd (1990) observed, modulator release is extremely sensitive to temperature. For unknown reasons, release *decreases* with increasing temperature. Using the standard concept of a temperature coefficient *Q*_{10}, the value of release *r* at temperature *T* is related to its value at 15°C by (4b) In their Fig. 6*B*, Vilim et al. (1996b) measured SCP and BUC release from motor neuron B15 at several temperatures. Incorporating *Eq. 4b* into *Eq. 3a* of Brezina et al. (2000c), obtaining the temperature-sensitive counterparts to *Eqs. 14a, 14c*, and *15* of Brezina et al., and fitting the last to the data of Vilim et al. (not shown) yielded *Q*_{10} = 5.3 × 10^{–2}. Thus release decreases about 19-fold over 10°C, and about 7.8-fold from 15 to 22°C. We inserted (5.3 × 10^{–2})^{(}^{T}^{–15)/10} into *Eq. 2a* so as to be able to model release, and hence modulator concentration and cAMP, in a temperature-dependent manner. We neglect the intrinsic temperature dependency of *C* and *A*, and also of the downstream effects modeled in Figs. 6, 7, 8, 9, 10, 11, because we expect that their probably much more modest temperature dependency will be largely irrelevant next to the extraordinarily large temperature dependency of release. We will return to this issue below.

We deal first with SCP released from motor neuron B15, for which most data are available. Figure 5*A* plots cAMP measurements [of Lloyd et al. (1984) and Whim and Lloyd (1989)] made after exposing ARC muscles to different fixed concentrations of exogenous SCP for 5 min, at which time *A*(*t*) will have risen reasonably close to its steady-state value, *A*_{∞} (see Hooper et al. 1994a). From *Eq. 4a*, when d*A*(*t*)/d*t* = 0, we have (4c) Fitting *Eq. 4c* to the data in Fig. 5*A* yielded *K _{A}*

_{,SCP}= 2.44 μM and

*k*

_{A}_{+,SCP}/

*k*

_{A}_{–,SCP}= 254.

Figure 5*B* plots all the available measurements [of Whim and Lloyd (1990)] of the time course of cAMP rise when motor neuron B15 was stimulated to fire in a particular known pattern (“*t* < *L*”), and its fall after the firing was discontinued early on, while cAMP was still rising (“*t* > *L*”). We simulated the same firing with the release model, *Eq. 2*, then fed the output (SCP release) through *Eqs. 3* and *4a*, and selected, in addition to the 2 values already obtained, values of *v, k _{C}*, and

*k*

_{A}_{–}to fit the modeled

*A*(

*t*) to the data as best as possible. This yielded

*v*

_{SCP}= 10 μl,

*k*

_{C,SCP}= 0.1 s

^{–1},

*k*

_{A–,SCP}= 0.1 s, and so

*k*

_{A+,SCP}= 25.4 s

^{–1}. As Fig. 5

*B*shows, we were able to fit the time course after the firing ended—first a brief continued increase in

*A*(

*t*), then a rapid fall, to basal values within 60 s— quite well. However, the relatively large values of

*k*and

_{C}*k*

_{A}_{–}that were required to make this fall fast enough also made

*A*(

*t*) rise too fast at the beginning of the firing. More exactly, the experimental measurements in Fig. 5

*B*appear to show a biphasic rise in cAMP, first a rapid rise within 60 s to a plateau around

*A*≈ 10 that is maintained for 2–3 min, then another rise within several minutes to

*A*> 200. Our model averages the 2 phases.

The fitted volume *v*_{SCP} = 10 μl is not very much smaller than the volume of the whole ARC muscle, which in 100–200 g *Aplysia* might be 15–40 μl (15–40 mg) (see Scott et al. 1997; Whim and Lloyd 1989), and in the 300- to 600-g animals used in the release studies of Vilim et al. (1996a,b 2000) presumably somewhat larger. This is consistent with the findings of Karhunen et al. (2001) that suggest that motor neuron B15 terminals release the modulators, in contrast to ACh, in a spatially diffuse manner.

When our modeled motor neuron B15 is fired in a behaviorally realistic pattern— e.g., in the standard reference pattern *d*_{intra} = 3.5 s, *d*_{inter} = 3.5 s, *f*_{intra} = 10 Hz, at 15°C—for a sufficiently long time, we find that the quasi-steady-state mean (period-averaged) SCP concentration, 〈*C*_{SCP}〉_{∞}, is about 130 nM, and 〈*A*〉_{∞} is about 14 (Fig. 5*A*). The latter value agrees well with the cAMP elevation actually measured with similar firing by Cropper et al. (1990a) and Whim and Lloyd (1990). Their data are reproduced in Fig. 5*C**, left* and *right*, respectively. Firing our modeled B15 in the exact pattern used in either case slightly overestimates the measured cAMP elevation, as Fig. 5*C* shows. On the other hand, the data of Hooper et al. (1994b) on the consequent activation of PKA, again allowing the effect of realistic B15 firing to be calibrated against that of known SCP concentrations (Fig. 5*D*), confirm a realistic 〈*C*_{SCP}〉_{∞} of the order of 130 nM.

For MM release from motor neuron B16, the data are more limited. Figure 5*E*, the counterpart to Fig. 5*A*, plots cAMP measurements [of Hooper et al. (1994a)] made with exogenous MM concentrations. The data are not sufficient to model the MM-induced elevation of cAMP fully, as we did for SCP. However, when B16 was fired by Hooper et al. in a behaviorally realistic pattern—the standard reference pattern *d*_{intra} = 3.5 s, *d*_{inter} = 3.5 s, *f*_{intra} = 20 Hz, at 15°C— cAMP was elevated about 1.7-fold. Figure 5*E* suggests that this elevation corresponds to a MM concentration of the order of 50–100 nM. Thus when motor neurons B15 and B16 are both firing in realistic patterns—B16 at higher frequency than B15—the concentrations of the released SCP and MM are quite similar. However, as noted in modulator release above, B16 releases considerably more MM than B15 does SCP under these circumstances. Consequently, when we simulated the realistic firing of B16 with our model, we found that we had to increase *v, k _{C}*, or both, above their B15–SCP values to keep the MM concentration as low as that of SCP. The correct 〈

*C*

_{MM}〉

_{∞}≈ 50–100 nM was obtained with

*v*

_{MM}= 35 μl and

*k*

_{C,MM}= 0.35 s

^{–1}. These larger values may be a reflection of the more extensive innervation of the ARC muscle by B16 than by B15 (Cohen et al. 1978; F. S. Vilim, personal communication).

We now return to the issue of temperature dependency. Although the data in Fig. 5*A* were obtained at room temperature, Whim and Lloyd (1990) performed a separate test of the temperature dependency of the elevation of cAMP by exogenous SCP. They found that 100 nM SCP_{B} elevated cAMP 62 ± 20-fold at 22°C and essentially identically, 57 ± 26-fold, at 16°C (means ± SE, *n* = 5). This strongly suggests that the relevant parameters in our model, *K _{A}, k_{A}*

_{+}, and

*k*

_{A}_{–}, do not vary with temperature over this range. In further support of the temperature independence of these parameters as well as of

*k*and

_{C}*v*, we obtained the same behaviorally realistic estimate of 〈

*C*

_{SCP}〉

_{∞}from Fig. 5

*A*, at room temperature, as from Fig. 5

*D*, at 15°C. All of the data in Fig. 5,

*C–E*were obtained at physiologically low temperature, at 15, 16, or 16–18°C.

Nevertheless, some uncertainty clearly remains in our modeling of cAMP (particularly as the available cAMP data are not all mutually consistent) and thus in our estimates of the ultimately important parameters, *k _{C}* and

*v*. Therefore, while we proceed with the best estimates of

*k*and

_{C}*v*that we have determined here, in Paper II we examine the impact of variation of these parameters on some of the main findings of the work.

modeling of the three modulator effects. The previous static model (Brezina et al. 1996) modeled separately the SCP and MM dose–response curves of the 3 effects on the Ca current, K current, and relaxation rate, and the convergence of SCP and MM on each effect. Here we develop a general formalism to cover both aspects as well as the dynamics of the effects.

The total magnitude of an effect **X** at time *t* is the sum of the individual effects of all *n* modulators in the system (5a) The individual effect **X*** _{i}* of each modulator

*i*is described by the schema (5b) with corresponding differential equation (5c) where is the concentration of the modulator

*C*modified, in general, by a Hill coefficient

_{i}*h*

_{X}_{,}

*,*

_{i}**X**

_{max,}

*is the maximal effect of the modulator,*

_{i}*k*

_{X}_{+,}

*and*

_{i}*k*

_{X}_{–,}

*are on and off rate constants, and*

_{i}*F*, a constant between 0 and 1, specifies the fraction of the effect of modulator

_{ji}*j*that occludes the effect of modulator

*i*. Always,

*F*= 1, and, in general,

_{ji,j=i}*F*≠

_{ji}*F*. Because SCP and MM do not appear to interact on any of the 3 effects in any complex way beyond simple additivity at low concentrations and partial or complete occlusion at higher concentrations (Brezina et al. 1994a,b, 1995, 1996), a simple constant

_{ij}*F*in the preceding formalism provides an adequate description.

_{ji}In the steady state, when d**X*** _{i}*(

*t*)/d

*t*= 0,

*Eqs. 5a*and

*5c*yield (5d) where

*K*

_{X}_{,}

*≡*

_{i}*k*

_{X}_{–,}

*/*

_{i}*k*

_{X}_{+,}

*. For completely noninteracting, additive modulators,*

_{i}*F*

_{ji}_{,}

_{j}_{≠1}= 0, and (5e) and for a single modulator (5f) Here, the most relevant case is that of 2 interacting modulators. For 2 modulators, subscripted 1 and 2,

*Eq. 5d*yields (5g) If the 2 modulators have the same maximal effect,

**X**

_{max}, and completely mutually occlude so that

*F*

_{21}=

*F*

_{12}= 1, (5h) Finally, if further

*K*

_{X}_{,1}=

*K*

_{X}_{,2}=

*K*

**(and**

_{X}*h*

_{X}_{,1}=

*h*

_{X}_{,2}=

*h*

**) to give a completely symmetrical situation (5i)**

_{X}In the case of a single modulator, with concentration *C _{i}* held constant for

*t*≥ 0,

*Eq. 5c*has the solution (5j) with

**X**

_{i}_{,∞}given by

*Eq. 5f*and a characteristic time constant (5k)

As can be seen for example in Figs. 6*B* and 9*B*, the simple first-order formalism of *Eqs. 5a–c* did not fully capture the real dynamics of the effects at all modulator concentrations. Specifically, with realistically constrained values of *k*_{X}_{+,}_{i}, k_{X}_{–,}* _{i}*, and

*K*

_{X}_{,}

*predicts a time constant that varies with modulator concentration to a significantly greater degree than was observed. For simplicity, we retained the first-order framework, and chose parameter values so as to model correctly, in the first instance, the dynamics of the effects observed in the behaviorally realistic range of modulator concentrations, below and around 100 nM.*

_{i}, Eq. 5kenhancement of calcium current. This effect is symbolized **C** and measured as percentage enhancement of the basal, unmodulated Ca current.

The static model (Brezina et al. 1996) has already established the steady-state dose–response relations of **C** for SCP and MM, identical in the 2 cases, described by *Eq. 5f* with **C**_{max} = 80%, *h*** _{C}** = 0.75, and

*K*

**≈ 3.0 × 10**

_{C}^{–6}M

^{h}**. These steady-state relations are shown in Fig. 6**

^{C}*C*(“

**C**

_{SCP,∞}”, “

**C**

_{MM,∞}”) together with the steady-state

**C**predicted for all combinations of SCP and MM from 10

^{–10}to 10

^{–4}M, with

*F*

**= 1, by**

_{C}*Eq. 5i*.

To model the dynamics of **C**, we performed a series of new experiments in voltage-clamped, dissociated ARC muscle fibers, examining the relaxation of **C** after a quasi-instantaneous step in modulator concentration to a new, maintained level. Figure 6*A* shows a representative experiment, Fig. 6*B* the group data, for “on” concentration steps from 0 to 100 nM or 1 μM. Because previous work (Brezina et al. 1994a) and preliminary experiments here suggested that, like the steady-state dose–response relations, the dynamics of **C** were likely to be identical for SCP and MM, we did not exhaustively test both SCP and MM at each concentration. Fitting in the first instance the 100 nM data, *Eq. 5j*, with the above value of *K*** _{C}**, yielded

*k*

_{C+}= 750 s

^{–1}M

^{–}

^{h}**and**

^{C}*k*

_{C}_{–}= 2.25 × 10

^{–3}s

^{–1}.

The converse, “off” concentration steps showed that **C** reversed slowly, indeed so slowly that the reversal was difficult to quantify accurately, although it was consistent with the time constant of 7.4 min predicted by *Eq. 5k* with the value of *k*_{C}_{–} already found.

The endogenously released modulators are not single peptides, but peptide families: motor neuron B15 releases 2 forms of SCP, and B16 releases, probably, 7 forms of MM. In most of our experiments here, however, we used just one form of each modulator, SCP_{B} and MM_{A}. This use is justified as follows. It appears very likely that each motor neuron obligatorily coreleases all of its modulator forms in an essentially invariant ratio (Brezina et al. 1995; Vilim et al. 1996a,b, 2000). Probably half of the released SCP mixture is SCP_{B}, and more than half of the released MM mixture MM_{A}. The other SCP in the mixture acts indistinguishably from SCP_{B} (Lloyd 1986); the other MMs act, in most respects, very similarly to MM_{A} (Brezina et al. 1995). Thus it is likely that the effects of a particular concentration of exogenously applied SCP_{B} or MM_{A} well reproduce those of a mixture of endogenously released SCPs or MMs with the same total concentration.

activation of potassium current. This effect is symbolized **K** and measured as percentage of the maximal current activated by, for example, 10 μM MM_{A} (Brezina et al. 1994b, 1995, 1996).

The static model (Brezina et al. 1996) established steady-state dose–response relations of **K** for SCP and MM, described by *Eq. 5f* with **K**_{max,MM} = 100%, **K**_{max,SCP} = 27%, *h*** _{K}** = 0.85, and

*K*

**≈ 1.85 × 10**

_{K}^{–6}M

^{h}**. Thus SCP and MM differ in their maximal activatable**

^{K}**K**. Here, in addition, we modified

*K*

_{K}_{,SCP}to correct the one significant discrepancy found by Brezina et al. (1996) between the predictions of the static model and the real observations made to verify it. Specifically, we set

*K*

_{K,SCP}≈ 1.85 × 10

^{–5}M

^{hK}(while retaining

*K*

_{K,MM}≈ 1.85 × 10

^{–6}M

^{h}

**) to fit the steady-state plot of contraction size shown here in Fig. 8**

^{K}*B*, which depends largely on

**K**, better to the experimental plot in Fig. 4

*A*of Brezina et al. (1996). The final steady-state dose–response relations for SCP and MM are shown in Fig. 7

*B*(“

**K**

_{SCP,∞}”, “

**K**

_{MM,∞}”) together with the steady-state

**K**predicted for all combinations of SCP and MM from 10

^{–10}to 10

^{–4}M, with

*F*

_{K}_{,SCP,MM}= 1 and

*F*

_{K}_{,MM,SCP}= 0.27, by

*Eq. 5g*.

To model the dynamics of **K**, we performed a series of new experiments in voltage-clamped, dissociated ARC muscle fibers, examining the relaxation of **K** after a quasi-instantaneous step in modulator concentration to a new, maintained level. In most experiments the fiber was held at –30 mV, the voltage within the physiological range where the K current is largest and most significantly modulates contraction size (Brezina et al. 1994b; Kozak et al. 1996; Orekhova et al. 2003), but the kinetics of the relaxation were not obviously different at –40 or –50 mV. Figure 7*A* shows the group data—simply averages of all the individual traces, so no single individual trace is shown—for “on” concentration steps from 0 to 100 nM or 10 μM MM_{A} (“*t* < *L*”), as well as “off” concentration steps from 10 μM MM_{A} to 0 (“*t* > *L*”; “off” relaxations from 100 nM had similar kinetics). Figure 7*A* also shows that, at the same concentration of 10 μM, SCP_{B} and MM_{A} elicited relaxations with similar kinetics, even though with different amplitudes. Fitting in the first instance the 100 nM MM_{A} “on” data and the “off” data, *Eq. 5j*, with the above values of *K*** _{K}**, yielded

*k*

_{+}= 1.9 × 10

^{5}s

^{–1}M

^{–}

^{h}**to (for both SCP and MM), 3.5 s**

^{K}^{–1}and

*k*

_{K–,MM}= 0.35 s

^{–1}.

Although the relaxations of the K current in these experiments were fast, they were not limited by the speed of application of the modulator— by the fact that the modulator step, inevitably, was not truly instantaneous. Thus when the same method was used to apply ACh, the directly ligand-gated current activated by ACh in these fibers (Kozak et al. 1996) turned on much faster still (Fig. 7*A*, “ACh”).

In these experiments and in the model we neglected desensitization of the K current, which is seen to any significant extent only with very large activation of the current (Brezina et al. 1994b).

modulation of contraction size. The effects **C** and **K** jointly determine the modulation of contraction size, symbolized **S** and measured as percentage of the basal, unmodulated contraction size. That is, **S** = 100% represents the basal contraction size; 0 ≤ **S** < 100%, depressed contractions; and **S** > 100%, potentiated contractions.

The static model (Brezina et al. 1996) has already provided a relation, which we continue to use here in the dynamic model, governing the combination of **C** and **K** to give **S** (6) where **S**_{max} = 400% is the maximal modulation of contraction size, *K*** _{S}** = 8% is a characteristic constant for the potentiation of contraction size with

**C**, κ = 9% is a characteristic constant for the depression of contraction size with

**K**, and

*a*= 8/3 is a constant necessary to ensure that

**S**(

**C**= 0,

**K**= 0) = 100%.

*Equation 6*is plotted in Fig. 8

*A*. Figure 8

*B*then uses

*Eq. 6*to compute, from Figs. 6

*C*and 7

*B*, the steady-state

**S**predicted for all combinations of SCP and MM from 10

^{–10}to 10

^{–4}M. The steady-state dose–response relations for SCP and MM alone are also shown (“

**S**

_{SCP,∞}”, “

**S**

_{MM,∞}”).

increase of relaxation rate. This effect is symbolized **R** and measured as a percentage of the maximal increase of relaxation rate produced by ≥1 μM SCP or MM.

The static model (Brezina et al. 1996) has already established the steady-state dose–response relations of **R** for SCP and MM, identical in the 2 cases, described by *Eq. 5f* with **R**_{max} = 100%, *h*** _{R}** = 0.7, and

*K*

_{R}≈ 3.64 × 10

^{–6}M

^{h}**. These steady-state relations are shown in Fig. 9**

^{R}*C*(“

**R**

_{SCP,∞}”, “

**R**

_{MM,∞}”) together with the steady-state

**R**predicted for all combinations of SCP and MM from 10

^{–10}to 10

^{–4}M, with

*F*

**= 1, by**

_{R}*Eq. 5i*.

To model the dynamics of **R**, we performed or reanalyzed a number of experiments with motor neuron B15- or B16-elicited contractions of the intact ARC muscle, examining the relaxation of **R** after a quasi-instantaneous step in exogenous modulator concentration to a new, maintained level. (The motor neuron firing patterns were such that endogenous modulators were not released in these experiments.) Figure 9*A2* shows a representative experiment.

Figure 9*A1* illustrates how the relaxation rate was measured, as the reciprocal of the time taken by the contraction to relax to 1/3 of its peak amplitude (see Brezina et al. 1995). The increase of relaxation rate produced by each modulator concentration was then expressed as a percentage of the maximal increase produced by 1 or 10 μM SCP or MM in that particular muscle, to give **R**. In cases where the maximal increase was not available, we used the average values from Fig. 10*B* of Brezina et al. (1995). For typical unmodulated contractions relaxing at 1 s^{–1} (relaxation time 1 s), for example, the average maximal relaxation rate achievable is about 4 s^{–1} (relaxation time 250 ms) (see also Brezina et al. 1996).

Figure 9*B*, “*t* < *L*”, shows the group data for “on” concentration steps from 0 to 10 nM, 100 nM, or 1 μM. In addition to SCP_{B} and MM_{A}, MM_{B} was used in some of these experiments because it activates the K current, and thus depresses contractions, much less than MM_{A}, while increasing the relaxation rate just like MM_{A} (Brezina et al. 1994b, 1995, 1996).

Figure 9*B*, “*t* > *L*”, shows group data for “off” concentration steps taken from Probst et al. (1994). In these experiments exogenous modulator was not applied, but rather endogenous SCP was released by stimulating motor neuron B15 to fire in the standard, behaviorally realistic pattern *d*_{intra} = 3.5 s, *d*_{inter} = 3.5 s, *f*_{intra} = 12 Hz for 10 min. The SCP concentration attained was therefore probably of the order of 100 nM (see modulator concentrations above).

Fitting in the first instance the 10 and 100 nM “on” data and the “off” data, *Eq. 5j*, with the above value of *K*_{R}, yielded *k*_{R}_{+} = 305 s^{–1} M^{–}^{h}** ^{R}** and

*k*

_{R–}= 1.11 ×10

^{–3}s

^{–1}.

In the intact ARC muscle, unlike in the dissociated fibers, the speed of access by exogenously applied modulators probably *does* limit the observed speed of some effects. Figure 10 shows simultaneous measurements of the depression of contraction size and the increase of relaxation rate by 1 μM MM_{A}. Although the depression of peak contraction is largely the result of activation of the K current (Orekhova et al. 2003), the depression in the intact muscle is significantly slower than the activation of K current in dissociated fibers (Fig. 7*A*). However, as Fig. 10 shows, the increase of relaxation rate is much slower still. Thus the data in Fig. 9*B* very probably reflect the true dynamics of **R**.

### Analytical and numerical methods

The complete dynamic model is composed of *Eqs. 1a-e, 2a-c, 3, 5a* and *5c* for **X** ∈ {**C**, **K**, **R**}, and *Eq. 6*, and the parameter values in Tables 1 and 2. Unless specified otherwise, *T* = 15°C.

Formulation and initial analysis of the model equations was carried out, in some cases, with the aid of *Mathematica* (Wolfram Research, Champaign, IL). Fitting to data was done using nonlinear least-squares Marquardt–Levenberg regression implemented within the SigmaPlot scientific graphing program (SPSS, Chicago, IL). Relaxation rates were measured from contraction records with a custom-written routine or using MiniAnalysis (Synaptosoft, Decatur, GA). To compute the output of the model, the model equations were integrated numerically using XPP, a general program for solving systems of differential equations written by G. B. Ermentrout (University of Pittsburgh), or, in later versions of the model, *Mathematica*. For numerical solutions of differential equations, *Mathematica* 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

^{–}

*, where*

^{p}*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

### Model overview

The model, constructed and fitted to experimental data as described in detail in methods, makes dynamic all of the variables *X*(*t*) in Fig. 1. Each of the relations between these variables, represented by the arrows in Fig. 1, is instantiated in one or more differential equations or instantaneous relations within the model. The firing patterns of the motor neurons B15 and B16, modeled as the firing frequencies *f*_{B15} and *f*_{B16} parameterized for cycles of the different types of feeding behavior over the whole hour-long meal as summarized in Table 2, can then be fed as input into the model and the corresponding trajectory of each of the modeled variables, and ultimately of the 2 outputs, the contraction size modulation **S** and the relaxation rate increase **R**, can be followed. Here the model is purely that of the modulation. To interpret **S** and **R** in absolute terms so as to predict the complete contraction shapes, the model has to be joined with that of the B15/B16-ARC neuromuscular transform, as we do in Paper II.

### The modulator effects have dynamics on very different time scales

The SCPs and MMs released by the firing of B15 and B16 exert 3 primary effects that we have modeled: they enhance the muscle's Ca current (model variable **C**), activate a K current (**K**) (**C** and **K** then combine instantaneously to give **S**), and increase the relaxation rate (**R**). Experiments that we performed or reanalyzed here (Figs. 6, 7, 9, and 10) to model the intrinsic dynamics of these 3 effects (their relaxation to steady state when their immediate input, SCP and MM concentration, is held constant) revealed that they have strikingly different time scales. Figure 11 summarizes this with a simulation of one of these experiments in the model. In Fig. 11*A*, when the concentration of MM, *C*_{MM}, is stepped from 0 to 100 nM (or, later, the reverse), **K** relaxes to the new steady state rapidly, on a time scale of the order of 1 s (with 100 nM MM, the modeled intrinsic time constant τ** _{K}** = 1.78 s). In contrast,

**C**and

**R**relax 2 orders of magnitude more slowly (τ

**= 154.6 s and τ**

_{C}**= 202.0 s with 100 nM MM). Figure 11**

_{R}*B*plots these relaxations against each other in the 3-dimensional (

**C**,

**K**,

**R**) space, and Fig. 11

*C*, combining

**C**and

**K**to give

**S**, in the 2-dimensional (

**S**,

**R**) space, the output space that we wish to study.

### Single-behavior meals

As a benchmark for subsequent more realistic meals, we first ran the model in complete, hour-long meals, with the full input trajectory of the cycle period *P*, but consisting of repeated cycles of just one type of feeding behavior. Figure 12 shows, for example, the behavior of the main variables of the model during an all-swallowing meal.

As swallowing begins and accelerates with food-induced arousal (Fig. 12, *P*), the rhythmic bursts of B15 and B16 firing come more often; the mean firing frequencies increase. As a result, SCP and MM release increases, and the modulators accumulate, to peak concentrations around 150 nM (*C*_{SCP}) and 100 nM (*C*_{MM}) during swallowing. Consequently, their 3 effects **C**, **K**, and **R** become strongly activated. However, because the depressing effect **K** is high, **S** is low; indeed, as Fig. 12 shows, the model predicts even some period of net depression of contractions. (Recall that **S** = 100% represents the basal contraction size; **S** < 100%, depressed contractions; and **S** > 100%, potentiated contractions.)

All this takes considerable time—≤10 min—to develop fully after the beginning of the meal. To some extent this time course reflects the gradual acceleration of the input firing pattern, but mainly the slow dynamics of downstream steps. **C** and **R** have slow intrinsic dynamics, of course, which contribute to the particularly slow development of those effects (see, e.g., Fig. 15*A* below). However, the fact that in Fig. 12 also **K**, which has fast intrinsic dynamics, and *C*_{SCP} and *C*_{MM} all develop quite slowly suggests that a common rate-limiting step occurs at the level of *C*_{SCP} and *C*_{MM} or above. *C*_{SCP} and *C*_{MM} have relatively fast intrinsic dynamics, however (see next), and the common component of the slow initial development of all of the variables appears to reflect largely the slow facilitation of release of the modulators, a phenomenon well documented in the real system (Fig. 4; Brezina et al. 2000c; Karhunen et al. 2001; Vilim et al. 1996a, 2000; see modulator release in methods).

Another noteworthy feature in Fig. 12 is that the rhythm of the bursts of B15 and B16 firing propagates through *C*_{SCP} and *C*_{MM}—which it would not do if these had slow intrinsic dynamics—to **K**, and thus to **S**. Thus because the cumulative dynamics of the steps upstream of it are relatively fast, and because its own intrinsic dynamics are fast, **K**, and so **S**, can respond to the input firing pattern on the time scale of the individual cycles. In contrast, **C** and **R**, with their slow intrinsic dynamics, do not sense the fast rhythm of the firing pattern at all.

As the meal progresses, *C*_{SCP} and *C*_{MM} gradually decrease again. This is partly because of the deceleration of the firing pattern with satiation, and partly because of the decay of modulator release that we have modeled as depletion of the releasable modulator pool. In any case, **C**, **K**, and **R** all decrease. However, **K** decreases more rapidly than **C**, and because of the asymmetry of their contributions to **S** (Fig. 8), **S** increases considerably in the later part of the meal.

Altogether, we see that early in the meal at the height of food-induced arousal (timepoint *a* in Fig. 12), when contractions are frequent, their relaxation rate is strongly increased and their size is somewhat depressed, whereas later in the meal during satiation (timepoint *b*), when contractions are less frequent, their relaxation rate is decreasing and their size is potentiated. A qualitatively similar progression can be seen (e.g., in Fig. 13*A*) in all-biting and all-swallowing meals.

### Trajectories through the output space

We next examined the trajectories of the single-behavior meals in the (**S**, **R**) output space (Fig. 13*A*; *curves of small symbols*). (In Fig. 13*A* and subsequent figures, each individual symbol represents one cycle.) It was immediately obvious that the trajectories of the all-biting, all-swallowing, and all-rejection meals are very different, and, strikingly, that they are not at all restricted to the steady-state region (*gray region* in Fig. 13*A*)—the region mapped out by the steady-state effects of all possible combinations of SCP and MM concentrations—found in our previous static analysis (Fig. 2*B*). Rather, the trajectories spend much of the meal moving through the region of lower **S** and **R** below the arch of the steady-state region. Even when they enter the steady-state region, they continue moving through it to the very end of the meal. In other words, the entire course of each trajectory is a transient that, very obviously at certain times, is far from the steady state.

For the presence of significant transients, where the trajectory remains far from the steady state a large part of the time, 2 conditions must be satisfied. First, clearly, the input to the system must change to move the system away from the steady state. In more realistic meals that frequently switch between the different feeding behaviors (below), this happens frequently. Even in the single-behavior meals, however, the input firing patterns are never stationary. At the beginning of the meal, there are large, fast changes in firing pattern. Then in the later part of the meal, the continuing change in *P* as satiation progresses prevents the trajectories from ever coming to rest. Furthermore, the gradual depletion and decay of modulator release means that, even if the firing patterns were stationary, the released modulator concentrations would progressively decrease. To demonstrate the operation of these 2 factors, in Fig. 14*A* we eliminated both satiation and depletion from the model. Now the trajectories do eventually come to rest at particular, fixed steady-state points (*large open circles* in Fig. 14*A*).

The second condition that must be satisfied for the presence of significant transients is that at least some of the dynamics must be slow. The fact that the trajectories readily travel outside of the steady-state region in Figs. 13*A* and 14*A* immediately shows that the requisite slow dynamics must occur at the level of the modulator effects, downstream of the concentrations *C*_{SCP} and *C*_{MM} in the model. For the steady-state region is defined with respect to these concentrations—it is the map of the steady-state effects of all possible combinations (*C*_{SCP}, *C*_{MM}). Because at every moment there obtains some combination (*C*_{SCP}, *C*_{MM}), if the downstream effects were in the steady state, the system would necessarily find itself in the steady-state region. The fact that the system is not always located there shows that at least some of the downstream effects are far from the steady state. With their slow intrinsic dynamics, **C** and **R** are the obvious effects in this respect.

However, this does not yet fully describe the mechanism that enables the system to travel out of the steady-state region. For the system that is far from the steady state corresponding to the current values (*C*_{SCP}, *C*_{MM}) can yet be located within the steady-state region, that is, at a steady-state point for some other (*C*_{SCP}, *C*_{MM}). This is the case, for example, in the later part of the all-swallowing meal after the trajectory has entered the steady-state region, yet is not in the steady state (see following text). And it is possible (as shown next) that the system might *always* thus be restricted to the steady-state region.

The major reason for the ready movement through the non-steady-state space is the great disparity between the dynamical time scales of the modulator effects—the fact that, in addition to the slow effects **C** and **R**, there is also the fast effect **K**. Examination of the model equations suggests that when all effects in a system such as this have the same time scale, trajectories are restricted to the steady-state region. To illustrate these conclusions, we slowed **K** to the speed of **C** and **R** in the model and ran it again in an all-swallowing meal (Fig. 14*A*, “Swallowing, slow **K**”). Now the trajectory, although it eventually reaches the same steady-state point as before, is forced to travel there through the steady-state region, over the arch rather than under it.

Close examination of Fig. 14*A* reveals that the 3 steady-state points are not exactly in the locations that might be expected: in the case of rejection, this is very obvious because the steady-state point is completely outside the steady-state region. Analysis of this discrepancy, in Fig. 14*B*, reveals that the mechanism of disparate dynamical time scales operates not only over multiple cycles of behavior, but even within individual cycles.

The model—and, plausibly, the real system—produces continuous waveforms **S**(*t*) and **R**(*t*), which represent the modulation that would occur at any time *t* if a contraction occurred then. The full trajectory [**S**(*t*), **R**(*t*)] is seen as the *thin oscillating curve* in Fig. 14*B*. The trajectory oscillates dramatically in the vertical, **S**, dimension, and little in the horizontal, **R**, dimension because, as already observed (Fig. 12), **K**, and so **S**, has dynamics fast enough to respond to the oscillations of modulator release and concentration in each cycle, whereas **R** does not. The modulatory state is then read out only at discrete times *t**, when a contraction, more specifically the peak contraction in each cycle, *does* occur. (In Paper II, the modulation is implemented in a more sophisticated but similar manner.) The trajectory of these points [**S**(*t**), **R**(*t**)] is shown as *blue diamonds* in Fig. 14*B*; this is the trajectory seen in *A*. In each cycle, the point of read-out lies very near the bottom of the full oscillation because the peak contraction comes at the end of a burst of motor neuron firing to which **K**, and so **S**, also has time to respond.

In the dynamical steady state, the full trajectory evolves to a stationary, vertically oscillating loop. Its mean does lie in the steady-state region, as Fig. 14*B* shows (“dynamical steady state, mean of effects”). Interestingly, even this mean is not exactly identical to the “true steady state” that would be obtained—in the manner in which the steady-state region was defined in the static analysis—with steady application of the corresponding mean modulator concentrations. The difference is due to, and a measure of, a measure of “pattern-dependency”: the fact that, in a step with a nonlinear input–output relation, a temporally patterned input will in general produce a different “amount” of output than will steady application of the corresponding mean input (see Brezina et al. 1997). More important, however, the discrete trajectory evolves to a steady state that, as throughout the trajectory, lies near the bottom of the full oscillation (“dynamical steady state, effects at peak contraction”), in the case of rejection completely outside the steady-state region. Thus thanks to the fast dynamics of **K**, even in the dynamical steady state the system can range widely over the output space, including into non-steady-state regions which would not be predicted from the steady-state analysis.

Because time is implicit in Figs. 13*A* and 14, these figures do not give a good sense of how far the system actually is at any point during the meal, especially once it enters the steady-state region, from the steady state corresponding to the current modulator concentrations (*C*_{SCP}, *C*_{MM}). To give this sense, Fig. 15*A* shows again the time course of the discrete read-out of the effects **C**, **K**, **S**, and **R** during the all-swallowing meal (*curves of black symbols*), compared with the true steady-state values of the effects computed from the mean modulator concentrations (*curves of blue symbols*), and from the peak modulator concentrations that coincide with the read-out of the effects at the peak of contraction (*curves of red symbols*), in each cycle.

We see in Fig. 15*A* that over the first few minutes of the meal, unable to follow the changes in modulator concentrations set up by the rapid changes in motor neuron firing (Fig. 12), the slow effects **C** and **R** lag far behind the steady state. The fast effect **K** lags comparatively little, of course, and, as already explained, it is this disparity that allows the system to move out through the non-steady-state space below the arch of the steady-state region in Figs. 13*A* and 14*A*. Later in the meal, when the firing patterns and mean modulator concentrations stop changing rapidly (in a single-behavior meal), **C** and **R**, and the system as a whole, catch up and track much closer to the steady state. It is noteworthy, however, that whereas **C** and **R**, which respond to the modulator concentrations only over many cycles, track most closely the steady state corresponding to the mean modulator concentrations, **K** does not. Because **K** is fast enough to react as the concentrations oscillate within individual cycles, it tends to depart from the steady state corresponding to the mean concentrations toward, even though it does not have time to actually reach, the steady state corresponding to the peak modulator concentrations. This is least evident in biting (not shown in Fig. 15*A*; however, see, e.g., timepoint *c* in Fig. 15*B*), more during swallowing (Fig. 15*A*; *b* in Fig. 15*B*), and most of all during rejection (*d* in Fig. 15*B*), as the bursts of firing of motor neuron B16, which release MM, the principal activator of the K current, grow longer.

Thus again because of the disparity between the dynamical time scales of the effects, even over long stretches of constant behavior the system as a whole can be far from any single true steady state corresponding to a single combination of modulator concentrations (*C*_{SCP}, *C*_{MM}). The different effects sense different aspects of the temporal pattern of the modulator concentration profile, essentially different concentration values, and so tend toward different steady states.

### Switches between behaviors

Next, we studied meals with one or more switches between the behaviors. The trajectory of a representative meal is super-imposed on those of the single-behavior meals in Fig. 13*A* (*large symbols*). The behavior of **C**, **K**, **S**, and **R** during the switches is shown in Fig. 13*B*.

The trajectory in Fig. 13*A* exhibits, even more explicitly, the features already described. It spends much of its time in the non-steady-state space. It is very obviously a transient, driven back and forth across the space by the switches in input firing pattern. Each switch is followed by rapid movement in the vertical, **S**, dimension that is largely complete in just one or two cycles (each symbol in Fig. 13*A* represents one cycle), but only much slower movement in the horizontal, **R**, dimension. Again, as Fig. 13*B* shows, this is due to the fast dynamics of **K**, and so **S**, as compared with **R**.

The switching trajectory makes explicit one feature that was already implicit in the single-behavior trajectories: history-dependency. History-dependency is implicit in the existence of a transient, whereby each successive bite, say, is associated with a different modulatory state than the previous bite. Even more obviously with the switching trajectory, however, the modulatory state associated with a cycle of a particular behavior depends heavily on when in the meal that cycle occurs and what sequence of behaviors, with diminishing importance back to the very beginning of the meal, preceded it.

History-dependency is possible to the extent that the system is far from the steady state. In the same way as for the all-swallowing meal in Fig. 15*A*, Fig. 15*B* shows the distance of the system during the switching meal from the true steady states corresponding to the current mean and peak modulator concentrations. The same mechanisms noted in Fig. 15*A* are at work in Fig. 15*B*. However, the system is on average even further from the steady state during the switching meal because a prolonged relaxation of the slow effects **C** and **R**, similar to that occurring at the beginning of the all-swallowing meal, follows each behavioral switch.

### Quasi-realistic meals

Switches like this occur frequently in the real feeding behavior of *Aplysia* (e.g., Kupfermann 1974; Morton and Chiel 1993a), but their precise statistical distribution within the meal, as many other statistical features of the behavior (see discussion in Paper II), remains to be determined. As a reasonable approximation, we simulated meals switching between biting, swallowing, rejection, and a null “behavior,” that is, pauses between the other behaviors (e.g., Hurwitz and Susswein 1992; Susswein et al. 1986), randomly 75 times in each hour-long meal, thus on average every 48 s. The trajectory of a representative meal is shown in Fig. 16*A* and all cycles from 21 meals in Fig. 16*B*. The *large symbols* show the mean for each type of behavior. There are statistically significant differences between the means in the vertical, **S**, dimension (for all pair-wise differences, *P* < 10^{–10} by Student's *t*-test), but not in the horizontal, **R**, dimension. Again the reason is that **K**, and so **S**, conforms rapidly to each new behavior, whereas **R**, reacting more slowly than the average interval between the behavioral switches, senses only the overall mix of behaviors. On average, then, a different state of modulation of contraction size, though not of relaxation rate, is associated with each type of feeding behavior.

## DISCUSSION

When the ARC motor neurons B15 and B16 fire in different patterns in the different feeding behaviors throughout the meal, they release different profiles of ACh to produce different shapes of basal contraction. This process, and its description, we have called the *neuromuscular transform* (Brezina et al. 2000b). At the same time, however, the motor neuron firing patterns release different profiles of SCP and MM to modulate the basal contractions to their final, behaviorally appropriate shape. We have constructed here a dynamic model of this “neuro-modulation” transform.

We have identified and briefly examined the major features of the operation of the neuromodulation transform. The most striking features are all consequences of the transform's combination of disparate dynamical time scales, specifically of the fact that while some of its effects, notably modulation of the relaxation rate, are very slow and react to the motor neuron firing only over many cycles of the behavior, one particular effect, modulation of the K current and thus of contraction size, is fast enough to respond to each individual cycle. The presence of slow effects suffices to make the trajectory of the modulatory state a transient throughout a typical meal, but superimposition of the fast effect allows the transient to range widely over the modulatory space, including large regions of it that are not accessible in the steady state. In our earlier static analysis of the modulation we emphasized the principle of coverage of the space by combinations of SCP and MM concentrations (Brezina et al. 1996); what we see here is therefore a complementary, and perhaps even more important, dynamic mechanism. A further consequence of the transient nature of the modulation is its pronounced history-dependency. All these features are especially prominent with switches between the behaviors, which are followed by rapid relaxations in the modulation of contraction size, but not that of relaxation rate. Finally, in quasi-realistic meals, each type of feeding behavior is associated, on average, with a different state of modulation of contraction size. A different modulatory state is indeed expected on theoretical grounds, as specific optimization of each distinct behavior (Brezina et al. 2000a).

However, to see, and discuss further, the functional significance of these features, we must add to the neuromodulation transform a model of the basal neuromuscular transform so as to be able to reconstruct the actual modulated contraction shapes. We do this in Paper II.

## DISCLOSURES

This work was funded by National Institutes of Health (NIH) Grants MH-50235 to K. R. Weiss and NS-41497 to V. Brezina. Some *Aplysia* were provided by the National Center for Research Resources National Resource for *Aplysia* at the University of Miami under NIH Grant RR-10294.

## Acknowledgments

B. Bank performed some of the experiments that were reanalyzed for Fig. 9*B*. We thank A. Kivenson and W. Bronson for help with some of the experiments in Figs. 6 and 7, F. S. Vilim for providing raw experimental data, G. B. Ermentrout for the program XPP, and two anonymous reviewers for very helpful suggestions.

## Footnotes

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 © 2003 by the American Physiological Society