JN Fuel your research with LabChart
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
 QUICK SEARCH:   [advanced]


     


J Neurophysiol 90: 2592-2612, 2003. First published July 9, 2003; doi:10.1152/jn.01091.2002
0022-3077/03 $5.00
This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow All Versions of this Article:
90/4/2592    most recent
01091.2002v1
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Right arrow Citation Map
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (15)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Brezina, V.
Right arrow Articles by Weiss, K. R.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Brezina, V.
Right arrow Articles by Weiss, K. R.

Neuromuscular Modulation in Aplysia. I. Dynamic Model

Vladimir Brezina, Irina V. Orekhova and Klaudiusz R. Weiss

Department of Physiology and Biophysics and Fishberg Research Center for Neurobiology, Mount Sinai School of Medicine, New York, New York 10029

Submitted 6 December 2002; accepted in final form 5 July 2003


    ABSTRACT
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 DISCLOSURES
 ACKNOWLEDGMENTS
 REFERENCES
 
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
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 DISCLOSURES
 ACKNOWLEDGMENTS
 REFERENCES
 
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 1997Go). 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. 1997Go; Weiss et al. 1993Go) 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 1974Go). 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. 1978Go; Kozak et al. 1996Go). 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. 1987aGo, 1990aGo; Lloyd et al. 1984Go; Vilim et al. 1996aGo; Whim and Lloyd 1989Go), and B16 releases the myomodulins (MMs) (Brezina et al. 1995Go; Cropper et al. 1987bGo, 1991Go; Vilim et al. 2000Go). 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 Ca2+ essential for contraction (Brezina and Weiss 1995Go; Brezina et al. 1994aGo, 1995Go; Cropper et al. 1987aGo,bGo, 1991Go; Lloyd et al. 1984Go; Whim and Lloyd 1990Go;); 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 Ca2+ influx (Brezina et al. 1994bGo, 1995Go; Cropper et al. 1987bGo, 1991Go; Orekhova et al. 2003Go); and 3) they increase the relaxation rate of the contractions probably by modulating the muscle's contractile machinery (Brezina et al. 1995Go; Cropper et al. 1990aGo; Heierhorst et al. 1995Go; Probst et al. 1994Go; Whim and Lloyd 1990Go). 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. 2000aGo). 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. 2A, "SCP" and "MM."



View larger version (47K):
[in this window]
[in a new window]
 
FIG. 1. Scope of our model of modulation in B15/B16-accessory radula closer (ARC) neuromuscular system. For details see INTRODUCTION. Dashed line indicates limits of previous static model of Brezina et al. (1996Go).

 


View larger version (42K):
[in this window]
[in a new window]
 
FIG. 2. Combinatorial modulation of ARC muscle contraction size and relaxation rate by small cardioactive peptides (SCPs) and myomodulins (MMs). A: representative experimental records. ARC muscle contractions were elicited by bursts of firing of motor neuron B16 ("B16") every 30 s (see METHODS). Shown are superimposed contractions under control conditions and in steady state after application of 10–8, 10–7, and 10–6 M ("c," "–8," "–7," and "–6", respectively) SCPB (left), MMA (middle), and of 10–7 M SCPB, 10–6 M MMA, or both (right; relaxation phases only). B: final output of previous static model: plot of combinatorial mapping of space of SCP and MM concentrations, C, to space of steady-state effects on size "S" and relaxation rate "R" of ARC muscle contractions. SCP and MM concentrations range from 10–10 to 10–4 M on log scales; large dots are plotted at intervals of 1 log unit, small dots of 0.1 log unit. Pairs of letters "a"–"d" show where corners of modulator space map to in effect space. Effects of SCP and MM alone are also shown. Red lines show how spaces are traversed as SCP concentration varies in presence of constant MM; blue lines, the converse. Space of effects was constructed by plotting against each other Figs. 8B and 9C, shown later in this study.

 



View larger version (57K):
[in this window]
[in a new window]
 
FIG. 8. Modulation of size of ARC muscle contractions by SCP and MM: model variable S. A: relation governing combination of C and K to give S. Plot of Eq. 6 in METHODS. B: modeled steady-state S predicted for all combinations of SCP and MM from 10–10 to 10–4 M, computed by passing relations plotted in Figs. 6C and 7B through that in A. Both A and B also show steady-state S in response to SCP and MM alone ("SSCP,{infty}", "SMM,{infty}"), computed similarly.

 


View larger version (39K):
[in this window]
[in a new window]
 
FIG. 9. Increase of relaxation rate of ARC muscle contractions by SCP and MM: model variable R. A: representative experiment. ARC muscle contractions were elicited by bursts of firing of motor neuron B16 ("B16") every 30 s (see METHODS) while 100 nM SCPB was applied. A1 illustrates how relaxation rate was measured, as reciprocal of time taken by contraction to relax to 1/3 of its peak amplitude. Increased relaxation rate could conveniently be visualized by scaling control contraction to same peak amplitude as modulated contraction (A1) or vice versa (A2). B: experimental group data and model fits. Plots "t < L" show measurements of R from experiments such as that in A (modulator was applied from t = 0 on) with 10 nM SCPB ({bullet}; means ± SE, n = 2), 100 nM SCPB ({blacktriangleup}; n = 13), 10 nM MMA or MMB ({circ}; n = 12), 100 nM MMA ({triangleup}; n = 7), and 1 µM MMA or MMB ({square}; n = 13). Both motor neurons B15 and B16 were used. Raw measurements of relaxation rate, as in A1, were converted to R as described in INCREASE OF RELAXATION RATE in METHODS. Plot "t > L" shows data from Fig. 4 of Probst et al. (1994Go) ({blacktriangledown}; means ± SE), tracing recovery of relaxation rate after release of endogenous SCP by firing of motor neuron B15 in pattern dintra = 3.5 s, dinter = 3.5 s, fintra = 12 Hz for 10 min, converted to R similarly. Curves show corresponding output of model, given by Eq. 5j in METHODS with parameter values for R. C: modeled steady-state R in response to SCP alone ("RSCP,{infty}"), MM alone ("RMM,{infty}"), and all combinations of SCP and MM from 10–10 to 10–4 M, given by Eq. 5i in METHODS.

 
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. 1996Go; see also Brezina and Weiss 1997Go; Brezina et al. 2000aGo). This model showed how the space of SCP and MM concentrations maps to the space of the 2 functional effects (Fig. 2B), 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. 2A, "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. 2B 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 1974Go; Susswein et al. 1976Go, 1978Go; Weiss et al. 1982Go). 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. 2003Go; 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. 2000aGo,bGo) to reconstruct, and test experimentally, the actual contraction shapes.


    METHODS
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 DISCLOSURES
 ACKNOWLEDGMENTS
 REFERENCES
 
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. 1994cGo). 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 (1994Go), Brezina et al. (1994aGo)]. The latter technique gave very rapid changes in modulator concentration (see Fig. 7A, "ACh"). All experiments were done at room temperature.



View larger version (27K):
[in this window]
[in a new window]
 
FIG. 7. Activation of K current in ARC muscle by SCP and MM: model variable K. A: experimental group data and model fits. Voltage-clamped, dissociated ARC muscle fibers were held at –30 mV while modulators were applied rapidly from pressure-ejection pipette. Nominal step of modulator concentration is indicated by gray bar "L"; actual speed of modulator arrival can be gauged from acetylcholine "ACh" response. Shown are means ± SE (SE indicated by zone of gray error bars around each mean curve) of 18 "ON" ("t < L") and "OFF" ("t > L") K-current responses from 8 fibers with 10 µM MMA, 5 "ON" responses from 4 fibers with 100 nM MMA, and 4 "ON" responses from 3 fibers with 10 µM SCPB. All K-current amplitudes were normalized to peak K current activated by 10 µM MMA in that particular fiber. Thin black curves show 2 corresponding outputs of model, "ON" response to 100 nM MMA and "OFF" response from 10 µM MMA, given by Eq. 5j in METHODS with parameter values for K. "ACh" values are means ± SE of 10 "ON" responses from 5 fibers to 10 µM ACh applied at –90 mV, inverted and scaled. B: modeled steady-state K in response to SCP alone ("KSCP,{infty}"), MM alone ("KMM,{infty}"), and all combinations of SCP and MM from 10–10 to 10–4 M, given by Eq. 5g in METHODS.

 

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 Ca2+, 55 Mg2+, 602 Cl, and 10 HEPES buffer (pH 7.6). See further Brezina et al. (1994bGo,cGo, 1995Go).

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 Ba2+ instead of Ca2+, 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. (1994aGo,dGo, 1995Go).

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. (1995Go) and Orekhova et al. (2003Go). These experiments all used a standard preparation and methods (see also Cohen et al. 1978Go; Cropper et al. 1987aGo; Vilim et al. 1994Go; Weiss et al. 1978Go). 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. 2A, 9A2). Such patterns of firing release ACh but not the motor neurons' endogenous peptide modulators (Cropper et al. 1990aGo; Vilim et al. 1996bGo). Contractions were monitored with an isotonic transducer. All experiments were done at room temperature.



View larger version (19K):
[in this window]
[in a new window]
 
FIG. 10. Speed of effects in intact ARC muscle. Size (bottom) and relaxation rate (top) of motor neuron B15-elicited ARC muscle contractions, with application of 1 µM MMA. Contraction size first decreases as K current rapidly activates, then begins to increase again, in part, as Ca current becomes progressively enhanced. See INCREASE OF RELAXATION RATE in METHODS.

 

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.


View this table:
[in this window]
[in a new window]
 
TABLE 1. Definitions and values of principal variables, parameters, and symbols

 

View this table:
[in this window]
[in a new window]
 
TABLE 2. Canonical set of behaviorally realistic motor neuron B15 and B16 firing patterns

 

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 1974Go; 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. 1990bGo; Morton and Chiel 1993aGo,bGo). 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 dintra, the interburst interval dinter, and the intraburst firing frequency fintra (see Brezina et al. 1997Go, 2000bGo). We also define the cycle period P = dintra + dinter. Bursting patterns of this simple kind (illustrated, e.g., in Fig. 3B, 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. 1997Go, 2000bGo) is that here the parameters P, dintra, dinter, and fintra 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 dinter is also the latency from t* to the start of the burst in each cycle, in particular in the first cycle.



View larger version (23K):
[in this window]
[in a new window]
 
FIG. 3. Behaviorally realistic motor neuron B15 and B16 firing patterns. A: trajectory of cycle period through meal. Symbols are relevant measurements, of bite or swallow latency or interval, from Fig. 6 of Susswein et al. (1976Go) ({square}; means ± SE, n = 11), Fig. 1 of Susswein et al. (1978Go) ({circ}; means ± SE, n = 10), Fig. 2 of Kupfermann and Weiss (1982Go) ({triangleup}; means only), and Fig. 11A of Rosen et al. (1989Go) ({triangledown}; means ± SE, n = 6), all scaled to L = 1 h. Curve is fit by Eqs. 1c-e in METHODS. B: expansion of indicated region of A, illustrating how model constructs firing patterns of motor neurons B15 (red lines) and B16 (blue lines) in individual cycles of biting, swallowing, and rejection. Upper two panels plot continuous functions P(t), dinter(t), dintra(t) [indicated here simply by difference between P(t) and dinter(t)], and fintra(t). At start of each cycle (i.e., at t = t*), model reads instantaneous values of these functions ({bullet}) and constructs from them (as indicated by thin slanting lines leading to bottom panel) detailed firing patterns of B15 and B16 over next cycle. Functions P(t), dinter(t), dintra(t), and fintra(t) can switch at any time between different behaviors (vertical dotted lines), but these changes will take effect only in firing patterns at next t*.

 

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

We now have to specify the functions P(t), dintra(t), dinter(t), and fintra(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 1992Go; Kupfermann 1974Go; Kupfermann and Weiss 1982Go; Rosen et al. 1982Go, 1989Go; Susswein et al. 1976Go, 1978Go; Weiss et al. 1986Go), 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. (1976Go, 1978Go), Kupfermann and Weiss (1982Go), and Rosen et al. (1989Go)] are replotted in Fig. 3A 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. (1990bGo)]. In rejection, however, it appears that a corresponding reasonable value is about 18 s (Weiss et al. 1986Go). 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. (1990bGo), 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 dintra, dinter, and fintra 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 dintra and fintra throughout the meal: for biting, dintra,B15 = 1.3 s, fintra,B15 = 10 Hz, dintra,B16 = 2 s, and fintra,B16 = 18 Hz; for swallowing, dintra,B15 = 3 s, fintra,B15 = 10 Hz, dintra,B16 = 3.7 s, and fintra,B16 = 18 Hz; and for rejection, dintra,B16 = 5 s and fintra,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 dinter is then left to vary according to the relation dinter(t) = P(t) – dintra.

Figure 3B 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 dintra,B15 = 3.5 s, dinter,B15 = 3.5 s, and fintra,B15 = 10–12 Hz for motor neuron B15, and dintra,B16 = 3.5 s, dinter,B16 = 3.5 s, and fintra,B16 = 20 Hz for motor neuron B16, which are often used in physiological studies (e.g., Vilim et al. 1996aGo,bGo, 2000Go; 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).



View larger version (54K):
[in this window]
[in a new window]
 
FIG. 14. Further analysis of single-behavior meals in (S, R) output space. A: single-behavior meals as in Fig. 13A, except that model was run with P(t) = 5.5 s instead of Eq. 1c–e and SSCP(t) = S0,SCP and SMM(t) = S0,MM instead of Eq. 2c, eliminating satiation and modulator depletion. To produce "Swallowing, slow K" trajectory, model was run in all-swallowing meal with same changes as well as K rate constants kK+ and kK divided by 250 and R rate constants kR+ and kR multiplied by 2 to equalize speeds of C, K, and R. B: all-rejection meal from A, but showing full continuous trajectory [S(t), R(t)] (thin oscillating curve) as well as its discrete read-out [S(t*), R(t*)] at peak of (unmodulated) contraction in each cycle (blue diamonds), same as shown in A. In dynamical steady state (see Brezina et al. 1997Go, 2000bGo), full trajectory evolves to stationary oscillating loop, of which "dynamical steady state, mean of effects" is period-averaged mean (<S>{infty}, <R>{infty}); discrete read-out evolves to "dynamical steady state, effects at peak contraction," ([S(t*)]{infty}, [R(t*)]{infty}), same as shown in A. "True steady state" is (S{infty}, R{infty})|<C>{infty}, steady state that would be obtained with steady application of corresponding mean modulator concentrations <CSCP>{infty} and <CMM>{infty}. See Trajectories through the output space in RESULTS.

 



View larger version (58K):
[in this window]
[in a new window]
 
FIG. 13. Trajectory of modulation in (S, R) output space during single-behavior meals and switches between behaviors. A: trajectories of all-biting, all-swallowing, and all-rejection meals (yellow squares, red circles, and blue diamonds, respectively, small symbols), and meal with biting->swallowing switch at a, swallowing->biting switch at b, and biting->rejection switch at c (corresponding large symbols). Each symbol represents one cycle; plotted are S(t*) and R(t*), values at start of cycle where (unmodulated) contraction would peak (see Paper II). Gray region is steady-state region from Fig. 2B, right. B: details of input firing patterns (bottom) and behavior of variables C, K, S, and R during switches in A. Complete model was run with input behavior fixed to biting, swallowing, or rejection, or switching between behaviors as shown.

 
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. (1996aGo,bGo). 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. 2000cGo), 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; kp+ and kp 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 Q10 = 5.3 x 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 kp+,B15 = 4.0 x 10–10, kp–,B15 = 3.4 x 10–3 s–1, and y = 3 [Model II of Brezina et al. (2000cGo)]. S0,SCP, the initial size of SSCP before any release, was found to be 541 fmol, without any evidence for replenishment of SSCP 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. (2000Go). The relevant data are reproduced in Fig. 4. Figure 4A 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. 3B, with behaviorally relevant parameter values. With the pattern in Fig. 4A, Vilim et al. obtained the MM and BUC outflow in Fig. 4B, and with systematic variations of the parameters dinter, fintra, and L, the data in Fig. 4C.



View larger version (27K):
[in this window]
[in a new window]
 
FIG. 4. Modulator release from motor neuron B16. A: standard reference firing pattern used by Vilim et al. (2000Go): burst duration dintra = 3.5 s, interburst interval dinter = 3.5 s, intraburst firing frequency fintra = 20 Hz, and so cycle period P = 7 s, mean (period-averaged) firing frequency < f > = 10 Hz; total stimulation length L = 1 h. B: time course of MM ({bullet}) and buccalin (BUC) ({circ}) outflow o, when motor neuron B16 was stimulated to fire as in A. Data from Fig. 3B1 of Vilim et al. (2000Go); means ± SE, n = 4. Vilim et al. actually measured outflow integrated over 2.5-min intervals every 5 min, alternately for MM and BUC, but this has been converted to outflow min–1. Plot was scaled correctly using absolute amounts measured (F. S. Vilim, personal communication). Smooth curves show best single-exponential fit to both MM and BUC values (with different scaling for MM and BUC) in interval 10 < t < 60 min [Eq. 17 of Brezina et al. (2000cGo): see MODULATOR RELEASE in METHODS]. C: total MM ({bullet}) and BUC ({circ}) outflow, O{infty}, resulting from whole block of firing of length L. C1: O{infty} as function of mean firing frequency< f >, with fixed L = 5 min. Data from Fig. 4, A3 and B3 of Vilim et al. (2000Go); means ± SE, n = 4. These experiments varied either dinter or fintra (as indicated here by 2 triplets of connected data points for each modulator), but because outflow appeared to depend simply on < f >, data have been combined. Vilim et al. actually presented data normalized per spike ({propto} O{infty}/< f >), but this has been converted to O{infty} again. Plot was scaled correctly using absolute amounts measured (F. S. Vilim, personal communication). C2: O{infty} as function of L with fixed < f > = 10 Hz, combining data from C1 and B. In both C1 and C2, curves are final best fits of Eq. 15 of Brezina et al. (2000cGo) as described in MODULATOR RELEASE in METHODS.

 

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. (2000cGo) (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. 2000Go), 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. (2000cGo) to the data in the interval 10 < t < 60 min in Fig. 4B, 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 -> {infty} and integrating the area under the curve (the experimental curve for 0 < t < 10 min, the fitted curve later), S0,MM = 2,690 fmol and S0,BUC = 697 fmol (indicated in Fig. 4C2). Next, with these values of S0, we fitted Eq. 16 of Brezina et al. (2000cGo) to the data in Fig. 4C1, 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. 4B yielded, by the argument in Brezina et al. (2000cGo), the "dissociation constant" Kp {equiv} kp/kp+ {approx} 8.9 x 107 s–1. Finally, with these values of S0, y, and Kp, fitting Eq. 15 of Brezina et al. (2000cGo) simultaneously to the data in both Fig. 4, C1 and C2 yielded kp+,B16 = 7.4 x 10–11 and kp–,B16 = 6.6 x 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 kp+ and kp 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. 2000cGo; Karhunen et al. 2001Go; Vilim et al. 1996aGo, 2000Go). 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 S0,SCP and S0,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 1991Go). 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. 4C1 here with Fig. 4A of Brezina et al. (2000cGo)], 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. 4B here with Fig. 2B of Brezina et al. (2000cGo)].

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 kC 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 kC, 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. 1994aGo; Cropper et al. 1990aGo; Hooper et al. 1994aGo; Lloyd et al. 1984Go; Whim and Lloyd 1989Go, 1990Go) and activate cAMP-dependent protein kinase (PKA) (Hooper et al. 1994aGo,bGo) 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.



View larger version (25K):
[in this window]
[in a new window]
 
FIG. 5. Effective concentrations of released SCP and MM, modeled through their elevation of cAMP and activation of protein kinase (PKA) in ARC muscle. A: A{infty} as function of CSCP: steady-state dose–response relation for elevation of cAMP, A, by SCP. cAMP content of ARC muscles (expressed relative to basal, unstimulated content) was measured after 5-min exposure to exogenous SCPB. Data from Fig. 3 of Lloyd et al. (1984Go) ({bullet}; n = 1) and Fig. 2A of Whim and Lloyd (1989Go) ({circ}; means ± SE, n = 3). Curve is best fit of Eq. 4c in METHODS. Effect of realistic B15 firing was established as described in MODULATOR CONCENTRATIONS in METHODS. B: time course of cAMP rise and fall in response to firing of motor neuron B15. Symbols are cAMP measurements [from Figs. 7 and 8 of Whim and Lloyd (1990Go); means ± SE, n = 2–3] made at various times while B15 was fired in pattern dintra = 4 s, dinter = 6 s, fintra = 25 Hz, L > 10 min ("t < L"; {circ}) and at various times after firing was discontinued at L = 3 min ("t > L"; {bullet}), at 22°C. Curves show corresponding output of model. These are numerical solutions of release model, Eqs. 2 in METHODS, with correctly patterned input given by Eq. 1a, for SCP release from B15 at 22°C, then Eqs. 3 and 4a. To simulate fall in cAMP after end of firing, firing was discontinued at correct amplitude of A, rather than after 3 min. C: cAMP elevation with behaviorally realistic B15 firing. cAMP measurements [left, from Fig. 3B of Cropper et al. (1990aGo), means ± SE, n = 6; right, from Fig. 9 of Whim and Lloyd (1990Go), means ± SE, n = 4] made after B15 was fired in patterns indicated for L = 10 min at 16–18°C (left; in running model, we assumed 17°C) or 16°C (right). In each case, corresponding output of model was computed as in B. D: steady-state dose–response relation for activation of PKA by SCP. PKA activity in ARC muscles (here expressed relative to basal, unstimulated activity) was measured after 2-h exposure to exogenous SCPB. Data from Fig. 2 of Hooper et al. (1994bGo) (means ± SE, n = 3–7). Curve is sigmoid like that in A, but with additional slope parameter (Hill coefficient) != 1. To establish effect of realistic B15 firing, B15 was fired by Hooper et al. in standard reference pattern dintra = 3.5 s, dinter = 3.5 s, fintra = 12 Hz, L = 10 min, at 15°C. E: A{infty} as function of CMM: steady-state dose–response relation for elevation of cAMP by MM. cAMP content of ARC muscles was measured after 5-min exposure to exogenous MMA. Data from Fig. 1A of Hooper et al. (1994aGo) (means ± SE, n = 4). Curve is best fit of Eq. 4c, but with additional slope parameter != 1. To establish effect of realistic B16 firing, B16 was fired by Hooper et al. in standard reference pattern dintra = 3.5 s, dinter = 3.5 s, fintra = 20 Hz, L = 10 min, at 15°C.

 

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. 3B of Weiss et al. (1979Go)]. A simple model for A is then

(4a)
where KA 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, kA+ is the maximal rate of cAMP production, and kA 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. (1996aGo,bGo, 2000Go) on which we based the release model was done at 15°C. As both Vilim et al. (1996bGo) and Whim and Lloyd (1990Go) observed, modulator release is extremely sensitive to temperature. For unknown reasons, release decreases with increasing temperature. Using the standard concept of a temperature coefficient Q10, the value of release r at temperature T is related to its value at 15°C by

(4b)
In their Fig. 6B, Vilim et al. (1996bGo) measured SCP and BUC release from motor neuron B15 at several temperatures. Incorporating Eq. 4b into Eq. 3a of Brezina et al. (2000cGo), 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 Q10 = 5.3 x 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 x 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.



View larger version (31K):
[in this window]
[in a new window]
 
FIG. 6. Enhancement of Ca current in ARC muscle by SCP and MM: model variable C. A: representative experiment. Ba currents were elicited by 100-ms voltage steps from –90 to 0 mV every 10 s in voltage-clamped, dissociated ARC muscle fiber while 1 µM MMA was applied. Subtracted, net Ba currents are shown (see METHODS). B: experimental group data and model fits. Measurements of C from experiments such as that in A (modulator was applied from t = 0 on) with 100 nM MMA ({triangleup}; means ± SE, n = 3), 1 µM MMA ({square}; n = 8), and 1 µM SCPB ({blacksquare}; n = 4). Curves show corresponding output of model, given by Eq. 5j in METHODS with parameter values for C. C: modeled steady-state C in response to SCP alone ("CSCP,{infty}"), MM alone ("CMM,{infty}"), and all combinations of SCP and MM from 10–10 to 10–4 M, given by Eq. 5i in METHODS. Relations CSCP,{infty} and CMM,{infty} were computed for zero concentration of other modulator, but for compactness are plotted here at 10–10 M.

 


View larger version (21K):
[in this window]
[in a new window]
 
FIG. 11. Intrinsic dynamics of effects C, K, S, and R have very different time scales. A: dynamics of C, K, and R after instantaneous "ON" and "OFF" steps in MM concentration. Computed using Eq. 5j in METHODS; scaled to same steady-state amplitude for comparison. B: dynamics from A plotted in space (C, K, R). Each dot represents 1 s. C: dynamics from A, with C and K passed through Eq. 6, plotted in output space (S, R).

 

We deal first with SCP released from motor neuron B15, for which most data are available. Figure 5A plots cAMP measurements [of Lloyd et al. (1984Go) and Whim and Lloyd (1989Go)] 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{infty} (see Hooper et al. 1994aGo). From Eq. 4a, when dA(t)/dt = 0, we have

(4c)
Fitting Eq. 4c to the data in Fig. 5A yielded KA,SCP = 2.44 µM and kA+,SCP/kA–,SCP = 254.

Figure 5B plots all the available measurements [of Whim and Lloyd (1990Go)] 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, kC, and kA to fit the modeled A(t) to the data as best as possible. This yielded vSCP = 10 µl, kC,SCP = 0.1 s–1, kA–,SCP = 0.1 s, and so kA+,SCP = 25.4 s–1. As Fig. 5B 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 kC and kA 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. 5B appear to show a biphasic rise in cAMP, first a rapid rise within 60 s to a plateau around A {approx} 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 vSCP = 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. 1997Go; Whim and Lloyd 1989Go), and in the 300- to 600-g animals used in the release studies of Vilim et al. (1996aGo,bGo 2000Go) presumably somewhat larger. This is consistent with the findings of Karhunen et al. (2001Go) 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 dintra = 3.5 s, dinter = 3.5 s, fintra = 10 Hz, at 15°C—for a sufficiently long time, we find that the quasi-steady-state mean (period-averaged) SCP concentration, <CSCP>{infty}, is about 130 nM, and <A>{infty} is about 14 (Fig. 5A). The latter value agrees well with the cAMP elevation actually measured with similar firing by Cropper et al. (1990aGo) and Whim and Lloyd (1990Go). Their data are reproduced in Fig. 5C, left and right, respectively. Firing our modeled B15 in the exact pattern used in either case slightly overestimates the measured cAMP elevation, as Fig. 5C shows. On the other hand, the data of Hooper et al. (1994bGo) on the consequent activation of PKA, again allowing the effect of realistic B15 firing to be calibrated against that of known SCP concentrations (Fig. 5D), confirm a realistic <CSCP>{infty} of the order of 130 nM.

For MM release from motor neuron B16, the data are more limited. Figure 5E, the counterpart to Fig. 5A, plots cAMP measurements [of Hooper et al. (1994aGo)] 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 dintra = 3.5 s, dinter = 3.5 s, fintra = 20 Hz, at 15°C— cAMP was elevated about 1.7-fold. Figure 5E 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, kC, or both, above their B15–SCP values to keep the MM concentration as low as that of SCP. The correct <CMM>{infty} {approx} 50–100 nM was obtained with vMM = 35 µl and kC,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. 1978Go; F. S. Vilim, personal communication).

We now return to the issue of temperature dependency. Although the data in Fig. 5A were obtained at room temperature, Whim and Lloyd (1990Go) performed a separate test of the temperature dependency of the elevation of cAMP by exogenous SCP. They found that 100 nM SCPB 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, KA, kA+, and kA, do not vary with temperature over this range. In further support of the temperature independence of these parameters as well as of kC and v, we obtained the same behaviorally realistic estimate of <CSCP>{infty} from Fig. 5A, at room temperature, as from Fig. 5D, 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, kC and v. Therefore, while we proceed with the best estimates of kC and 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. 1996Go) 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 Xi of each modulator i is described by the schema

(5b)
with corresponding differential equation

(5c)
where is the concentration of the modulator Ci modified, in general, by a Hill coefficient hX,i, Xmax,i is the maximal effect of the modulator, kX+,i and kX–,i are ON and OFF rate constants, and Fji, a constant between 0 and 1, specifies the fraction of the effect of modulator j that occludes the effect of modulator i. Always, Fji,j=i = 1, and, in general, Fji != Fij. 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. 1994aGo,bGo, 1995Go, 1996Go), a simple constant Fji in the preceding formalism provides an adequate description.

In the steady state, when dXi(t)/dt = 0, Eqs. 5a and 5c yield

(5d)
where K<