Previous Article | Next Article 
Volume 16, Number 12,
Issue of June 15, 1996
pp. 4017-4031
Copyright ©1996 Society for Neuroscience
A Dynamic Network Simulation of the Nematode Tap Withdrawal
Circuit: Predictions Concerning Synaptic Function Using Behavioral
Criteria
Stephen R. Wicks1,
Chris J. Roehrig1, and
Catharine H. Rankin2
1 Program in Neuroscience and 2 Department
of Psychology, University of British Columbia, Vancouver, British
Columbia, Canada V6T 1Z4
ABSTRACT
- INTRODUCTION
- MATERIALS AND METHODS
- RESULTS
- DISCUSSION
- FOOTNOTES
- REFERENCES
ABSTRACT 
The nematode tap withdrawal reflex demonstrates several
forms of behavioral plasticity. Although the neural connectivity that
supports this behavior is identified (Integration of mechanosensory
stimuli in Caenorhabditis elegans, Wicks and Rankin, 1995, J
Neurosci 15:2434-2444), t
he neurotransmitter phenotypes, and hence
whether the synapses in the circuit are excitatory or inhibitory,
remain uncharacterized. Here we use a novel strategy to predict the
polarity configuration, i.e., the array of excitatory and inhibitory
connections, of the nematode tap withdrawal circuit using an
anatomically and physiologically justifiable dynamic network simulation
of that circuit. The output of the modeled circuit was optimized to the
behavior of animals, which possessed circuits altered by surgical
ablation by exhaustively enumerating an array of synaptic signs that
constituted the modeled circuit. All possible polarity configurations
were then compared, and a statistical analysis was used to determine
whether, for a given synaptic class, a particular polarity was
associated with a good fit to behavioral data. The results from four
related experiments were used to predict the polarities of seven of the
nine cell classes of the tap withdrawal circuit. In addition, the model
was used to assess possible roles for two novel mechanosensory
integration neurons: DVA and PVD.
Key words:
Caenorhabditis elegans;
mechanosensory;
habituation;
lesion;
inhibition;
reflex
INTRODUCTION
Reflexes have long been recognized as the
substrate of many aspects of organismal behavior in diverse phyla from
invertebrates to mammals, including primates (Sherrington, 1906
; Carew
et al., 1972; du Lac et al., 1995). Reflexive behaviors are often
produced by small neural architectures that are amenable to empirical
and computational investigation (Selverston, 1985). Thus, reflexes can
be used to delineate the mechanisms underlying various forms of
behavioral plasticity, such as adaptation (Corfas and Dudai, 1990),
inhibition (Rankin, 1991; Fischer and Carew, 1993; Vu et al., 1993),
habituation (Groves et al., 1969; Krasne, 1969; Carew et al., 1972;
Wicks and Rankin, 1995b), and associative conditioning (Kandel and
Schwartz, 1982; Steinmetz et al., 1992). One system that is complex
enough to exhibit sophisticated behaviors yet simple enough to
accommodate a precise neural model of a reflexive behavior is the
nematode Caenorhabditis elegans.
The nematode C. elegans has been the subject of intensive
biological study over the past three decades (Brenner, 1974; Wood,
1988) and has become a model system within which the relationship
between genetics and behavior has been extensively studied. The adult
nematode possesses only 302 neurons, each completely described in terms
of location and morphology (Ward et al., 1975; Ware et al., 1975; White
et al., 1986; Hall and Russell, 1991). The synaptic connectivity of the
approximately 5000 chemical synapses, 600 gap junctions, and 2000 neuromuscular junctions has been encoded in machine-readable form
(Achacoso and Yamamoto, 1992). The extensive genetic analysis of the
worm has confirmed that a number of basic biological mechanisms, such
as second-messenger signaling (Gross et al., 1990; Lu et al., 1990;
Mendel et al., 1995; Segalat et al., 1995), synaptic release (Nonet et
al., 1993), and other cell-signaling events (Stern and DeVore, 1994),
are conserved between it and more complex neurobiological systems.
These advantages have been exploited by modeling researchers;
chemotaxic (Lockery et al., 1993) and locomotory behaviors in the worm
have been modeled previously (Niebur and Erdös, 1991, 1993;
Erdös and Niebur, 1993).
The demonstration that the worm is an adaptive system (Rankin et
al., 1990) is interesting because sophisticated nematode neurobiology
and genetics allow the possibility of describing the nature of this
plasticity at a circuit and cellular/molecular level. Many of the
demonstrations of learning in this organism have been made through the
study of the simple tap withdrawal reflex, for which the underlying
neural circuitry is identified (Wicks and Rankin, 1995a). Although the
detail of the anatomical data used to describe this circuit is
considerable, the functional polarities of the neurons of this circuit
are unknown. This report describes an array of putative polarities for
synapses of the tap withdrawal circuit. These predictions were made by
optimizing the behavior of a modeled circuit under conditions of
degradation isomorphic to laser ablation of the circuitry underlying
the tap withdrawal behavior of real animals. Additionally, this model
was used to investigate roles for two novel mechanosensory integration
neurons.
MATERIALS AND METHODS
The model. The model used was a physiologically
motivated one. However, in the absence of detailed physiological data
from C. elegans, it was necessary to make a number of
extrapolations from the related nematode Ascaris
lumbricoides. These assumptions are presented in physiological
rather than mathematical form to ensure that they are realistic.
Furthermore, preliminary investigations suggested that polarity
predictions based on the modeled circuit were more strongly determined
by circuit connectivity than the exact values of parameters used. Thus,
approximate ranges for these parameters rather than precise values were
derived. The effects of varying some of the more uncertain parameters
were then assessed by rerunning the same experiments with values of
these parameters varied over three orders of magnitude.
The circuitry was constructed by extracting connectivity data from
AY's Neuroanatomy for Computation (Achacoso and Yamamoto, 1992). This
data indicated not only the presence or absence of a set of synaptic
contacts between a pair of neurons, referred to here as a synaptic
class, but also incorporated the actual number of documented electrical
and chemical connections within that synaptic class. Each synaptic
contact within a class of synapses was assigned the same reversal
potential and conductance as all other synapses within that class. This
enabled the simple construction of complex circuits in which all
documented synapses (including all bilateral asymmetries) were included
in the model. It was assumed that the functional efficacy of a synaptic
class was correlated with the number of contacts observed within that
synaptic class. Thus, circuits constructed in this way possessed
connections with weights determined by anatomical criteria. These
weights were not varied further in this model; only the reversal
potential, which determined the sign of the connection, was allowed to
vary. The complete connectivity of the modeled tap withdrawal circuit
is shown in Figure 1.
Fig. 1.
The complete connectivity of the tap withdrawal
circuit. The circuit consists of seven sensory neurons (shaded
circles), nine interneurons (unshaded circles), and two
motorneuron pools (not shown), which produce forward and backward
locomotion (triangles). Chemical connections are indicated
by arrows, with the number of synaptic contacts being to the
width of the arrow. Gap junctions are indicated by
dotted lines. Every connection represented in this figure
was also represented in the model. This representation is useful for
identifying connection asymmetries, which might underlie the origin of
oscillations that control locomotion and are hidden in simpler views of
the circuitry.
[View Larger Version of this Image (59K GIF file)]
The model was based on all available physiological and anatomical data
from C. elegans and the related nematode,
Ascaris. The physiological parameters used in the derivation
of the data presented in this report are shown in Tables
1 and 2. The model was implemented in
Objective-C on Intel-486, HP series 9000 and NeXT computers running
NEXTSTEP software, and was integrated using fourth-order Runge-Kutte
(Press et al., 1989) to an accuracy of 0.5%.
Table 1.
Average simulation parameters
| Neuron parameters |
Value |
Units
|
|
| Membrane resistance |
See Table 2 |
Ohms
|
| Membrane capacitance |
See Table 2 |
Farads |
| Membrane
leakage potential |
 0.035 |
Volts |
| Synaptic parameters
|
| EPSP reversal potential |
0.00 |
Volts |
| IPSP reversal
potential |
0.048 |
Volts |
| Synaptic
conductance |
6.00E-10 |
Siemens
|
| VRange |
0.035 |
Volts |
| Gap
junction conductance |
5.00E-09 |
Siemens |
| Tap parameters
|
| Pulse rest |
0 |
Amps |
| Phasic pulse |
1.00E-11 |
Amps
|
| Start time |
0.01 |
Sec |
| Duration |
0.3 |
Sec
|
| Tonic pulse |
2.50E-10 |
Amps |
|
|
The list of physiological parameters used in the four experiments
run in this report are summarized. For a more detailed discussion of
the origin of these values, see Materials and Methods.
|
|
The neuron. The neurons of C. elegans have simple
morphologies, which are preserved across individuals. Many neurons
consist of a single unbranched process, and few have more than two
branches (Chalfie and White, 1988). Electrophysiology on C. elegans cells is still in its infancy (however, see Raizen and
Avery, 1994; Avery et al., 1995), and little is known about the
membrane characteristics of its neurons. However, electrophysiology has
been done on Ascaris, a larger nematode related to C. elegans (Davis and Stretton, 1989a,b). Its dorsal and ventral
nerve cords have been reconstructed and show considerable similarity to
those of C. elegans, and homologs of C. elegans
motor neurons have been found in Ascaris (Chalfie and White,
1988; Stretton et al., 1992). For this model, electrophysiological data
from Ascaris was used to determine model parameters.
Evidence from Ascaris suggests that signal propagation in
C. elegans neurons is likely accomplished electrotonically,
without classical all-or-none action potentials. Intracellular
recordings of Ascaris motorneurons and interneurons show no
evidence of action potentials, and it has not been possible to evoke
them (Davis and Stretton, 1989a). Membrane resistivity in
Ascaris is unusually high (60-300
k
cm2) and is within the range that would
permit signal propagation without action potentials. Niebur and
Erdös (1993) have used Ascaris data to do detailed
computational studies of the electrotonic characteristics of C. elegans neurons and have shown that the integration of C. elegans locomotion can be accounted for by purely electrotonic
signals.
Davis and Stretton (1989a) have measured specific membrane resistances
(
m) and intracellular resistivity
(i) in Ascaris. In four
motorneurons, m varied from 89 to
251 kcm2, and
i from 79 to 314 cm. We assumed
that membrane properties in C. elegans were similar and used
an average of the four measurements:
i = 180 cm and
m = 150 kcm2. We assumed a specific membrane
capacitance of 1 µF/cm2, a standard value for a
lipid bilayer (Rall, 1989). These membrane properties were adapted to
C. elegans anatomy by using the surface area of each cell
(see Table 2). Each neuron's branching morphology is given in Wood
(1988) and White et al. (1986). This, together with measurements of
electron micrographs in White et al. (1976, 1986), were used to
determine average process lengths and diameters. Diameters varied from
0.2 to 1.0 µm, and an average value of d = 0.5 µm was
used.
Process lengths were taken from diagrams in Wood (1988), assuming a
standard worm length of 1 mm. Soma diameters were taken from camera
lucida drawings in Wood (1988). Soma diameters varied from 2 to 10 µm; we used an average diameter of 5 µm. From these data, a total
membrane surface area for each cell was computed, and the resulting
total membrane capacitance and resistance for the entire cell was
derived (see Table 2). For simplicity, it was assumed that cells were
isopotential. Because the length constant (Rall, 1989)
|
(1)
|
was generally longer than the process (data not shown), this was
reasonable (also see Table 2). Thus, a neuron's membrane potential,
V, was governed by the usual single-compartment membrane
equation (Segev et al., 1989):
|
(2)
|
where Cm is the total
membrane capacitance for the cell,
Rm is the total membrane leakage
resistance for the cell, VLEAK is the
leakage potential of the cell, ISYN is the
current attributable to a synaptic input, and
IEXT is any injected current. A value of
35 mV was used for the leakage potential in these cells (R. E. Davis,
personal communication).
Gap junctions. The anatomical reconstruction of the nematode
nervous system allowed the identification of both electrical and
chemical synapses (White et al., 1986). Gap junctions were modeled as
ohmic resistances, where current flowing into cell i from
cell j was given by:
|
(3)
|
where 
ij was the
conductance of the gap junction. Niebur (1988) used a standard
conductance per unit area value of 1 S/cm2
(Bennett, 1972) and used unpublished micrographs to determine the area
of each gap junction. He reported that gap junctions ranged from 0.2 to
2 µm long and were 0.5 µm wide (Niebur, 1988). In our model, a
standard gap junction length of 1 µm was assumed, with a resulting
conductance of 5 nS for all gap junctions. In some experiments, this
value was increased or decreased by a factor of 10 to test the
sensitivity of the model's predictions to the precise value of the
conductance used.
Synapses. Synaptic classes consisted of a number of
individual synaptic contacts. The number of contacts in each class was
extracted from an anatomical database of C. elegans synaptic
connectivity (Achacoso and Yamamoto, 1992). The identification of
chemical synapses from the anatomical reconstruction was done by
identifying presynaptic specializations and inferring postsynaptic
partners based on proximity; no postsynaptic specializations were
evident in the electron microscope reconstructions (White et al.,
1986). Any error associated with this technique would tend to
overestimate the number of chemical synapses in the organism, however,
and any spurious synaptic classes would not be duplicated with high
frequency. Therefore, such synapses, although present in the model,
would not be assigned a large synaptic weight and thus would not have a
large impact on the response of the circuit to stimulation. Each
modeled synapse represented a class of synaptic contacts with total
synaptic conductance
the ``weight''given by the product of the
number of individual contacts within the class and the individual
synaptic conductance.
The synapse model used was based on the graded synapse model used by
Lockery and Sejnowski (1992) in the leech local bending circuit.
However, it was extended to explicitly include the synaptic reversal
potential as well as the conductance. Postsynaptic current was
attributable to gated channels in the postsynaptic membrane with inward
current given by:
|
(4)
|
where g(t) is the synaptic conductance of
the postsynaptic membrane, which was gated according to presynaptic
potential, and ESYN is the reversal
potential for the synaptic conductance, which was assumed to be
constant. For excitatory synapses, a reversal potential of 0 mV was
used, and for inhibitory synapses 45 mV was used (R. E. Davis,
personal communication). For simplicity, it was assumed that all
synapses made by a single cell possessed the same reversal potential.
This amounts to assuming that all synapses made by a given presynaptic
cell were of the same polarity and class. Computationally, this reduced
the number of optimized parameterseach of which possessed two
possible valuesto the number of neurons in the modeled circuit. It
was further assumed that all modeled synapses functioned as fast
ligand-gated channels. It was possible that some anatomically defined
synapses were modulatory and acted via slow second-messenger systems,
or that synaptic function was altered by the milieu
interieur (Harris-Warrick et al., 1992); however, we assumed that
these modulatory effects did not affect a single tap withdrawal
response.
In the leech local bending reflex, Lockery and Sejnowski (1992)
observed multiple time courses in some of the postsynaptic responses,
and to model this they used a slow (10 msec) and a fast (1500 msec)
conductance, each governed by its own first-order equation. Preliminary
versions of this model used a fast (10 msec) synaptic time constant,
but no significant differences were noted in the results of circuits
containing these synapses and simulations, which used an instantaneous
synapse. We therefore used a synaptic conductance that depended only on
the presynaptic membrane potential:
|
(5)
|
where g
represents the
steady-state conductance in response to a presynaptic voltage.
Synaptic activation. No direct recordings have yet been made
in C. elegans to determine properties of synaptic
activation. In recordings made from Ascaris commissural
motorneurons, Davis and Stretton (1989b) demonstrated that synaptic
transmission is graded and transmitter is released tonically between
both excitatory and inhibitory motorneurons and postsynaptic muscle and
motorneurons. They found that changes in postsynaptic potential were
related to presynaptic depolarizing current by a sigmoidally shaped
curve and that the presynaptic resting potential lies approximately in
the middle of the voltage-sensitive range of synaptic transmission.
Dynamic network simulations based on graded synaptic transmission have
been described previously (Lockery and Sejnowski, 1992; DeSchutter et
al., 1993). We assumed that synaptic activation and transmission in
C. elegans was similar to Ascaris; namely, that
it was graded and sigmoidally shaped with presynaptic potential and was
tonically active with the equilibrium potential in the middle of the
voltage-sensitive range. Accordingly, we used a symmetric sigmoidal
function to model the steady-state postsynaptic membrane
conductance:
|
(6)
|
where 
is maximal postsynaptic membrane
conductance for the synapse, VEQ is the
presynaptic equilibrium potential at the middle of the
voltage-sensitive range, and VRANGE is the
presynaptic voltage range over which the synapse activated. We used a
value of
|
(7)
|
so that the conductance changed from 10 to 90% of its maximal
value over a presynaptic voltage range of
VRANGE. Note that because synapses were
tonically active, a cell's equilibrium potential was not defined
solely by its resting membrane potential, but rather was determined
from the fixed point of the entire system of equations governing the
circuit. This was computed before each run in the following way.
Equilibrium potential. The assumption that tonically active
synapses were active in the middle of their voltage-sensitive range at
the equilibrium potential implied that the postsynaptic conductance
g(t) was one-half its maximal value
when the circuit was at equilibrium. Let
Vi denote the membrane potential for
neuron i and similarly for other quantities pertaining to
neuron i. Let Iij denote
the synaptic current flowing into neuron i from neuron
j across a single synapse and let
ij denote the total number of synaptic
connections from neuron j to neuron i. Similarly,
let Îij denote current flow
across a single gap junction in the direction from neuron j
to neuron i and let
ij denote the total number of
gap junctions between neuron j and neuron i.
Finally, let Ii denote the external
current flow into cell i (caused by either sensory
stimulation or external current injection). Then, the entire system is
given by:
|
(8)
|
|
(9)
|
|
(10)
|
|
(11)
|
|
(12)
|
where N is the number of neurons in the circuit and
ij is the synaptic time constant.
At equilibrium, Vi = VEQi, and
dVi/dt, dgij/dt and
Ii are zero. Synaptic conductances
are tonic and are at their half-activation at equilibrium so that
|
(13)
|
After algebraic manipulation, this yields a system of
linear equations that can be solved using standard Gaussian elimination
to find VEQi (Press et
al., 1988):
|
(14)
|
where Aij is the
ith row and jth column of matrix A and
is given by:
|
(15)
|
|
(16)
|
and b is a vector that is given by:
|
(17)
|
The computed equilibrium potential of cells varied
within a physiologically realistic range from 47 to 0 mV, with a mean
of 24 mV and an SD of 13 mV. This did not vary appreciably from cell
to cell, but rather depended on the circuit's polarity configuration
and ablation condition. As an aside, the system (Eq. 14) can also be
inverted to explicitly give VLEAK in terms
of VEQi:
|
(18)
|
Synaptic parameters. To determine values
for VRANGE and , the
synapse model was fitted to published measurements (Davis and Stretton,
1989b) of Ascaris muscle cell postsynaptic response to
presynaptic current injection as detailed below. Thus, values for
VRANGE and for
particular Ascaris synapses were obtained. We assumed that
C. elegans synapses activated over voltage ranges similar to
Ascaris synapses. However, the maximal synaptic conductance,
, needed to be adapted to C. elegans. We
assumed that represented the product of a synaptic
conductance per unit area and a synaptic area. In the case of a synapse
mediated by a single population of ion channels,
would be equivalent to the single-channel conductance of an open
channel multiplied by the total number of channels. To adapt the
value from Ascaris to C. elegans, we assumed that C. elegans synapses had
similar unit-area conductances and accordingly scaled the Ascaris
by a factor to account for the presumed difference
in synaptic areas. We assumed that the total synaptic area between two
cells was proportional to the length of the process, which we estimated
by the ratio of body lengthsapproximately 1/250. As this represents
only a gross approximation, the value of used in
these studies was varied over three orders of magnitude in different
experiments (see Results).
Modeling Ascaris monosynaptic responses. Data
from Davis and Stretton (1989b; their Figs. 13 and 14) show the
postsynaptic response of an Ascaris dorsal muscle cell (DM)
to current injected into a presynaptic excitatory motorneuron, DE1, and
the response of a ventral muscle cell (VM) to current injected into a
presynaptic inhibitory motorneuron, VI. Both of these response profiles
were sigmoidal in shape, were centered approximately at the resting
potential, and had asymptotic saturated postsynaptic responses at the
extremes of positive and negative presynaptic current injection.
Therefore, the sigmoidal tonic response model presented here was well
suited to fitting these synaptic responses.
Davis and Stretton (1989a,b) placed a recording electrode in a muscle
cell within the output zone of the motorneuron, and an injecting
electrode at the ventral end of a commissural process leading to the
synapse. The measured input resistance of the motorneuron was used to
obtain the resulting membrane potential at the point of current
injection, and an infinite cable model was used to determine the
membrane potential at the presynaptic site. Because the recordings that
were used to determine input resistance were made at the same ventral
end of the commissure as was injected for the synaptic response
measurements (R. E. Davis, personal communication), and because the
input resistance was approximately constant over the relevant range of
injected current (Davis and Stretton, 1989a), it is possible to
directly use the measured input resistance to determine the membrane
potential at the point of current injection.
Davis and Stretton (1989a) determined the motorneuron cable properties
by fitting their measurements along the commissure to an infinite cable
model (Rall, 1989) and found that the length constant was unusually
high (
~ 8 mm)approximately the same magnitude as the length of
the process. With such a large length constant, it is possible that the
cable's branching morphology and sealed ends played a significant role
in determining these cable properties (Rall, 1977), suggesting that a
sealed-end cable model might be more appropriate. However, for
consistency we used an infinite cable model with cable constants as
determined by Davis and Stretton (1989a) to reproduce their measured
voltage response along the commissure.
Specifically, the presynaptic depolarization in response to an injected
current is given by:
|
(19)
|
where L is the distance from the point of
current injection to the synapse and RPRE
is the input resistance at the point of injection. For the DE1-DM
synapse, the distance, L, was 5-8 mm; for the VI-VM
synapse, it was 0.5-2.5 mm (R. E. Davis, personal communication). We
used the mean of 7 and 1 mm, respectively. The input resistances for
the DE1 and VI motorneurons were reported to be 6 and 17 M,
respectively (Davis and Stretton, 1989a).
According to the sigmoidal tonic response model, the steady-state
plateau response of the postsynaptic muscle is given by:
|
(20)
|
where VLEAK and
ESYN pertain to the postsynaptic muscle
cell, RPOST is the postsynaptic cell's
input resistance,
g(VPRE)
is the steady-state synaptic conductance, and n is the total
number of synapses between the motorneuron and the muscle cell. The
input resistance of VM and DM was measured to be 0.18-0.50 M, with
a mean of 0.3 M (R. E. Davis, personal communication). Values of
VLEAK = 35 mV and
ESYN = 0 mV for excitatory and
ESYN = 45 mV for inhibitory reversal
potentials were used (R. E. Davis, personal communication). A light
microscope study of dye-injected muscle cells in Ascaris
suggested that the DE1 motorneurons make approximately 5-10 synapses
to each DM and the VI motorneurons make approximately 8-16 synapses to
each VM (J. Donmoyer, personal communication). Values from the high end
of these ranges were adopted for n in Equation 20 because
muscle cells form gap junctions with one another (J. Donmoyer, personal
communication), and thus it was possible that synapses from neighboring
muscle cells played a role in determining a muscle cell's postsynaptic
potential.
Equation 20 can be rearranged to yield values of
VPOST explicitly in terms of
VPRE:
|
(21)
|
where
|
(22)
|
The change in postsynaptic potential was therefore given
by:
|
(23)
|
where VEQ is the equilibrium
potential achieved by the postsynaptic cell under unstimulated
in-circuit tonic synaptic input (see Eqs. 14-18) and is given
by:
|
(24)
|
Equations 19 and 21, 22, 23, 24 define a nonlinear function for
postsynaptic membrane potential in terms of presynaptic injected
current and contain two unknown parameters: and
VRANGE. Levenberg-Marquardt's method
(Press et al., 1988) was used to fit this function to the data from
Davis and Stretton (1989b).
We obtained good results to the fit for the inhibitory VI-VM synapse.
This fit included the reversal potential
ESYN as a fit parameter; this improved the
fit substantially without significantly changing the reversal potential
(48 mV as opposed to 45 mV). The results of this fitting procedure
yielded values of = 150 nS,
VRANGE = 52 mV,
ESYN = 48 mV and were stable under
various initial conditions. The DE1-DM fit was less precise because the
steady-state conductance was not precisely a symmetric sigmoid. We
therefore manually explored a number of parameter values and attempted
to reproduce the approximate range and amplitude of activation. For
this we obtained values of VRANGE = 20 mV
and = 90 nS.
To adapt these values to C. elegans, an average of the two
VRANGE values was used to estimate the
activation range (35 mV), and the from the VI-VM
synapse fit was scaled by the 1/250 ratio of body lengths to yield a
maximal conductance of 0.6 nS for an individual C. elegans
synapse.
The gearbox. This model does not explicitly incorporate
nematode locomotion; these issues have been dealt with adequately
elsewhere (Niebur and Erdös, 1991, 1993; Erdös and Niebur,
1993). Rather, this report concentrates specifically on sensorimotor
integration. However, because a behavioral variable was used to
optimize the modeled output, it was necessary to rigorously define the
relationship between the animal's locomotion and activation of the
circuitry that controls that behavior. This issue was addressed with
simple assumptions, which were consistent with work on the modeling of
nematode locomotion (Niebur and Erdös, 1991, 1993; Stretton et
al., 1992; Erdös and Niebur, 1993) and current theories of tap
withdrawal circuit function (Chalfie et al., 1985; Wicks and Rankin,
1995a). The output of the tap withdrawal circuit was assumed to control
locomotory behavior primarily through the action of the interneurons
AVB and AVA. These two interneurons make electrical connections with
motor neurons all along the ventral cord of the worm. The AVA
interneurons make gap junctions with the motor neurons AS, VA, and DA,
which are presumed to excite backward locomotion; the AVB interneurons
form gap junctions with the motor neurons VB and DB, which are presumed
to excite forward locomotion. Ablation of these cells almost completely
destroys an animal's ability to move forward (in the case of AVB
ablations) or backward (in the case of AVA ablations) (Chalfie et al.,
1985; Wicks and Rankin, 1995a). Thus, it was simply assumed that the
degree to which an animal reversed was proportional to the
depolarization of the AVA interneuron and inversely proportional to the
depolarization of the AVB interneuron. Forward locomotion in response
to tap was also proportional to this value; a lower propensity to
reverse was equivalent to a higher propensity to accelerate. The exact
nature of this proportionality was not defined because in
vivo it will be modulated by a number of neural, hydrostatic, and
physical forces that are beyond the scope of this endeavor. The
gearbox, i.e., the transformation equation that was used to convert
depolarization of AVA and AVB into behavior, was simply:
|
(25)
|
The integration was calculated from the time of the tap
stimulation until either the end of the simulation or until the
integrand changed sign. Additionally, the test for a change of
integrand sign was suppressed for a grace period of 100 msec to allow
for initial transients after the tap.
One consequence of the gearbox assumption is that, because of
uncertainty regarding the exact nature of the proportionality between
the output of the AVA and AVB interneurons and the magnitude of the
evoked behavior, comparisons of model data and empirical data must be
limited to relative changes in response magnitude. Thus,
such comparisons were made between data profiles that had been
normalized about the mean of that polarity configuration's response
level (see below). This measure detected changes in the levels of
responding to tap produced by an ablation series, without being
sensitive to the absolute response magnitude of a particular circuit
configurationinformation which in any case is meaningless in the
context of the gearbox assumption.
Strategy for neuron polarity determination. The animal's
response to a light mechanosensory stimulus has been intensively
studied (Chalfie and Sulston, 1981; Chalfie et al., 1985; Wicks and
Rankin, 1995a). Specifically, Chalfie and colleagues have described the
circuitry that underlies the worm's reflexive response to a light
touch to either its head or tail. This touch circuit was the starting
point for the delineation of the circuitry responsible for the
animal's response to a mechanically delivered stimulus (i.e.,
``tap'') to the side of the substrate on which the animal moves. An
intact worm will generally respond to a tap stimulus with a cessation
of forward motion, a reversal through some distance, and a resumption
of forward locomotion in a new direction. This response has been termed
the tap withdrawal reflex (Rankin et al., 1990). The magnitude of the
tap withdrawal reflex can be quantified by measuring the distance
through which an animal reverses.
The circuitry underlying the tap withdrawal reflex (see Fig. 1) (Wicks
and Rankin, 1995a) has been identified using the laser ablation
technique (Sulston and White, 1980; Avery and Horvitz, 1987, 1989).
Ablations of neurons in this circuit can quantitatively and
qualitatively alter an animal's response to tap. Three neuron classes,
each composed of one or two cells, transduce the tap stimulus and
segment the animal into two mechanosensory fields; anterior
mechanosensory stimuli are transduced by the ALM and AVM cell classes,
and the posterior input is transduced by the PLM cell class. The AVA
and AVB interneurons are required for normal locomotion, and the AVD
and PVC interneurons couple anterior and posterior mechanosensory
input, respectively, onto AVA and AVB.
The ablation strategy previously used to describe the tap withdrawal
circuit resulted in the accumulation of data sets that represented the
worm's response to a tap stimulus when individual neurons in the tap
withdrawal circuit were destroyed with a laser mircobeam (see Fig.
2) (Wicks and Rankin, 1995a). For example, the ablation
of either the ALM sensory neuron class alone or in conjunction with the
AVM sensory neuron resulted in animals that accelerated forward, rather
than reversed, in response to a tap. Furthermore, the magnitude of both
the reversal and acceleration behaviors elicited by a tap stimulus was
modulated by ablation of the tap withdrawal circuit neurons. For
example, the reversal response of animals lacking the PLM sensory
neuron was larger than the reversal response of control animals with an
intact circuit. This difference was presumed to be attributable to the
loss of excitatory electrical input from the PLM sensory neuron to the
PVC interneuron, as well as chemical input to the AVA and AVD
interneurons (see Fig. 1). This chemical input could be either
excitatory or inhibitory, and that variable should, in theory, affect
the magnitude of the reversal response observed. Similarly, the
ablation of both the ALM and AVM sensory neurons resulted in larger
accelerations in response to a tap than did the ablation of ALM
alone.
Fig. 2.
Mean response magnitudes of the tap withdrawal
reflex of seven groups of animals. These data were used to optimize the
array of underlying functional polarities of the modeled circuit. Some
ablations (for example, the removal of PLM) resulted in larger reversal
responses than in control animals; other ablations resulted in
consistent accelerations in response to tap (indicated by a negative
reversal magnitude for the ALM and AVMALM groups). Note that the
acceleration measure is a change in velocity, whereas the reversal
measure is a distance; these are not directly comparable. Thus, the
ordinate represents a number that is proportional to the amount of
forward or backward locomotion. It was assumed that the pattern of
response represented on this figure was dictated partially by the array
of synaptic signs, which constituted the underlying circuitry (see
Materials and Methods, Strategy for neuron polarity determination).
These data are adapted from Wicks and Rankin (1995a).
[View Larger Version of this Image (17K GIF file)]
The connectivity of the worm's nervous system has been well described.
However, the signs or polarities of the synapses, excitatory or
inhibitory, which determine the behavior produced by a given circuit,
are unknown. This set of synaptic signs is referred to in this report
as the polarity configuration of that circuit. The same connectivity
may result in many distinct behavioral outcomes depending on the
particular polarity configuration that the circuit possesses. The data
that represent the animal's tap withdrawal behavior after neuronal
ablationpresented in Figure 2reflect the polarity configuration of
the worm's tap withdrawal circuit; it is a representation of how the
circuit that is defined by that set of polarities responds, in various
states of degradation, to a tap stimulus. Another circuitone with
identical connectivity but a different polarity configurationwould
not necessarily yield the same profile of tap responses as a result of
that series of lesions. Indeed, each possible polarity configuration of
the tap withdrawal circuit might result in a distinctive pattern of
behavior in response to the ablation of individual, or small sets of,
neurons. This profile of behavior (see Fig. 2), representing the real
animal's response to a tap in a variety of ablation conditions, was
used to find the best match within the space of all possible polarity
configurations using a computational model of the tap withdrawal
circuit.
We constructed a model of the tap withdrawal circuit and produced
lesions in that circuit that were isomorphic to the ablation conditions
represented in Figure 2. We then determined the profile of behavior for
all seven ablation conditions predicted by the model for an arbitrary
polarity configuration. The fitness of that polarity configuration was
determined by an error function with three terms (see Eq. 26), each of
which used a least-squares error approach to compare the model data
profile for that configuration and the ablation data profile (see Fig.
3). Then, all possible polarity configurations were
exhaustively enumerated and sorted according to their fitness. The
response profiles (modeled and empirical) were compared by first
separating those conditions that were associated with reversal
responses from those that were associated with accelerations. Then
these two profiles were standardized by expressing each data set as a
series of Z-scores around the mean of each set of responses. This was
done because we were interested in the relative change in
the withdrawal response magnitude as a consequence of ablation, rather
than the absolute value of the tap withdrawal response
magnitude.
Fig. 3.
Two different response profiles from experiment 1one representing a good fit and one representing a poor fitalong
with the empirical data to which they were compared. These data sets
are expressed as three sets of standardized Z-scores used to evaluate
the relative modulation by ablation of three behavioral measures:
REVERSAL MAGNITUDE, ACCELERATION MAGNITUDE, and
RESPONSE TYPE. The first set of Z-scores incorporates those
ablation configurations that result in a reversal response in the
empirical data set shown in Figure 2. It was used to determine the
error associated with the relative modulation of reversal magnitude as
a result of ablation between that data set and each polarity
configuration from the model. The second set of Z-scores similarly
assessed the error associated with the relative modulation of
acceleration magnitude as a result of ablation. A third set of Z-scores
evaluated the qualitative fit between the model and the empirical
profiles with respect to response type. It assessed whether the gearbox
output (Eq. 25) was lower, on average, for the acceleration profile
than for the reversal profile, as was assumed to be the case for the
intact animal. The fitness of a given configuration was calculated by
summing the least-squared error between the model and empirical
profiles for each of these three sets of standardized data. It is clear
that the fitness of the modeled circuit to the behavioral data can be
strongly modified by altering the array of underlying polarities in the
modeled circuit.
[View Larger Version of this Image (23K GIF file)]
The first error term in the fitness function reflected how closely the
ablation-induced modulation of modeled gearbox output matched the
corresponding modulation of reversal behavior in real animals. The
reversal profile comprised those ablation conditions that resulted in
reversals in the real animal: the intact, PLM, PVC, PVD, and
AVM conditions. Figure 3 shows three data sets that have been
standardized to allow comparison. The empirical and modeled reversal
profiles were each standardized as described above and compared. The
reversal error term was simply the sum of the least-squared errors
between the two reversal profiles. The second error term in the fitness
function reflected how well the modulation of the acceleration behavior
of the model produced by ablation matched the corresponding modulation
of acceleration behavior in real animals. The acceleration profile thus
comprised the two ablation conditions, which resulted in animals that
consistently accelerated in response to tap: the ALM and
AVMALM ablations. Again, both the empirical acceleration profile
(see Fig. 3) and the modeled acceleration profile were standardized as
described above, and the least-squared error between the two
acceleration profiles was derived for each configuration.
It was necessary to separate these two sets of comparisons because the
acceleration and reversal behaviors were quantified with distinct
measures (Wicks and Rankin, 1995)an acceleration is a change in
velocity, whereas a reversal is a distance. Although these two
responses were not quantitatively comparable, it was still possible to
qualitatively evaluate the ALM and AVMALM ablation conditions with
respect to the reversal profile in the model. Any output from the
gearbox (see Eq. 25) associated with the acceleration profile should,
to be considered in accord with data from the real animal, be on
average lower than the output of the gearbox associated with the
reversal profile. It is not possible to make any statement with respect
to how much smaller this gearbox output should be because of the
incommensurability of the acceleration and reversal measures, but
accelerations should be associated with relatively lower gearbox
outputs than reversals. Thus, a third error term was derived to
evaluate whether the response type produced by the model was in accord
with the empirical data inasmuch as the gearbox output for these
conditions was on average lower than for those ablation conditions
associated with the reversal profile. The mean gearbox value for the
reversal conditions (control, PLM, PVD, PVC, and AVM) was
calculated for both the empirical data set in Figure 2 and the model
data set. A similar mean acceleration value was calculated for the
acceleration conditions in the two data sets. These two sets of values
were converted to Z-scores, and a least-squares error between the model
and empirical response type profile was calculated. Because there were
only two terms in each of these distributions, this comparison was
strictly qualitative. If the gearbox output for the acceleration
conditions was, on average, lower than the gearbox output for the
reversal conditions, the two normalized distributions would be
identical, and the error contributed by this term would be zero. On the
other hand, if the two acceleration conditions resulted in gearbox
output that was, on average, larger than that produced by the reversal
conditions, the error contributed by this term would be positive and
furthermore, would be insensitive to how much larger the output was.
Examples of two different response profiles from a typical
experimentone representing a good fit and one representing a poor
fitare shown in Figure 3, along with the empirical data to which they
were compared.
These three error terms were summed for each configuration. Thus, the
fitness measure that was used to sort the list of polarity
configurations was given by:
|
(26)
|
where each term represented the least-squared error associated
with the comparison of the model and empirical profiles for that
behavioral measure. After summing the three error terms associated with
each polarity configuration, the entire list of polarity configurations
was sorted according to the resultant fitness function. The top 50 sorted polarity configurations from a single experiment are shown in
Figure 4. The interpretation of this list was
complicated by the fact that it was both exhaustive and complete. The
top n predictions, although sorted according to absolute
fitness, may not have been significantly different from one another
statistically. To address this issue, the polarity of each cell was
analyzed independently. A given neuron was characterized as excitatory
or inhibitory if it was determined that a given polarity for that cell
was clustered at or near the top of the list. For example, a cell that
was predicted to be inhibitory for the top 50 configurations was more
likely to actually be inhibitory than a cell whose predicted polarity
alternated through the top 50 configurations, even if in the case of
the best configuration, it was also inhibitory.
Fig. 4.
Sample polarity configurations. The top 50 polarity configurations sorted according to error from experiment 1 are
shown. This circuit did not include the DVA interneuron, and hence
there were 256 possible configurations (26) in
the complete sorted list. Thus, the top 10% of the list reported in
Table 3A consists of the top 26 polarity configurations shown in this
figure. The PVD sensory neuron class was not externally stimulated
during this run. A polarity consistent with that which resulted from
statistical considerations is shown as a lightly shaded box;
a polarity that is not consistent is shown in an unshaded
box. No polarity predictions were made for the AVA or DVA neurons.
These columns are darkly shaded. In this experiment, the
tenth and sixteenth configurations are entirely consistent with the
consensus configuration predicted in this report.
[View Larger Version of this Image (118K GIF file)]
A sign test (Siegel, 1956) is a nonparametric statistic that can be
used to assess the probability that the two possible polarities (+1 or
1 in this case) appear with equal frequency within a tested fraction
of this sorted list. A significant sign test indicates that one or the
other sign clusters within that fraction of the sorted list at above
chance levels. Thus, several one-sample sign tests were used for
each neuron to determine whether a given polarity for that cell
clustered at a higher-than-chance frequency near the top of the sorted
list of polarity configurations. This analysis was repeated for each of
several fractions of the sorted list of configurations. Because the
sorted list of configurations was complete, i.e., each configuration
differed from all others, and because a significant prediction at one
fraction of the list was indicative of a trend that might also be
detected at other nearby fractions of the list if they were analyzed,
this constitutes a conservative analysis of the data.
This paper reports the results of analyses based on four fractions of
the list of sorted polarity configurations. Three of these (the top
10%, first quartile, and top half) are arbitrary, but informative, and
the fourth (designated ``
'') was based on the shape of the
frequency distribution of the fitness function that was used to sort
the list of polarity configurations. The fraction was defined as
that fraction of the list of polarity configurations whose members
possessed a fitness value that was >1 SD above the mean of the fitness
frequency distribution. Thus, the fraction of the list that was
defined by was dependent on the distribution of the error term used
to sort the list. These fitness distributions are shown in Figure
5.
Fig. 5.
Fitness frequency distributions from four
experiments. The 256 possible configurations in experiments 1 and 3 and
the 512 possible configurations from experiments 2 and 4 were sorted
according to the fitness measure described in Materials and Methods.
The y-axis represents the number of configurations with a
given error in each of 24 error bins. The frequency distributions are
multimodal. The fraction of the sorted list of configurations that
corresponds to the level (>1 SD below the mean of the
distribution) lies to the left of the indicated level in each case.
The mean of each distribution is also indicated.
[View Larger Version of this Image (34K GIF file)]
Modeling the tap stimulus. The tap stimulus was modeled as a
phasic depolarization of the sensory neurons PLM, ALM, and AVM, which
have been shown to mediate the response to tap in the intact animal
(Wicks and Rankin, 1995a). However, it is likely that the tap does not
represent the only input to the circuit in the intact animal. For
example, the neuron PVD has been thought of as a stretch receptor (E. Hedgecock, cited in Way and Chalfie, 1988) or a background
mechanosensory input detector (Wicks and Rankin, 1995a). In addition,
the interneuron DVA receives considerable synaptic input from
chemosensory circuitry, and therefore may be a cell that modulates the
activity of the tap withdrawal circuit according to the chemical milieu
(Wicks and Rankin, 1995a). Therefore, four complete experiments were
run to determine to what extent variations in the way these two cells
were represented altered the nature of the polarity predictions. The
effect of the inclusion of DVA was assessed by comparing two
experiments in which DVA was included (experiments 2 and 4) with two
experiments in which DVA was absent (experiments 1 and 3). Within each
of these two conditions, the effect of the nature of the stimulation
that PVD and DVA (when present) received was assessed. Two stimulation
parameter sets were used. In the tonic condition (experiments 3 and 4),
the PVD and DVA neurons were activated by a low (one-quarter pulse
input magnitude) tonic stimulation that was continuously present and
not correlated with the phasic tap input to the other touch neurons.
These parameters were chosen to mimic the effects of an unspecified
mechanosensory and chemosensory input to PVD and DVA, respectively. In
the phasic condition, these two neuron classes were not explicitly
stimulated (experiments 1 and 2).
Preliminary versions of these results have been reported previously in
abstract form (Roehrig et al., 1994).
RESULTS
The results of the four experiments are summarized in Table
3A-D. The table represents a series of probability
values for sign tests conducted on the indicated percentages of the
sorted lists of polarity configurations in each of four experiments:
(A) the exclusion of DVA with no stimulation of PVD; (B) the inclusion
of DVA in the circuit and no stimulation of DVA and PVD; (C) the
exclusion of DVA with tonic activation of PVD; and (D) the inclusion of
DVA and tonic stimulation of DVA and PVD. Taken together, the results
of the four experiments made significant predictions for the polarities
of seven of the nine cell classes that constitute the tap withdrawal
circuit. The strength of the prediction for a given neuron was roughly
correlated with the density of chemical connections that that neuron
makes with other cells in the circuit (White et al., 1986; see Fig. 1).
Because it was the reversal potential of the chemical connections that
was the critical factor in determining the fitness measure, no
prediction would have been expected for neurons that make sparse
chemical connections. All of the predictions from each of the four
experiments were consistent. For example, if a cell was predicted to be
inhibitory in the first experiment, any further predictions that were
made in either the same or subsequent experiments were also inhibitory.
Because there were no disparities between the four experiments, a
predicted polarity configuration based on consensus results was
derived. This consensus configuration is shown in Figure
6. Because the polarity of two cells (DVA and AVA) were
not reliably predicted at any level of significance, there are multiple
configurations that are consistent with this consensus
configuration.
6 cm2)
