Chem. Senses 24: 497-508,
1999
© Oxford University Press 1999
A Kinetic Model of the Transient Phase in the Response of Olfactory Receptor Neurons
Division of Insect Biology, ESPM, 201 Wellman Hall, University of California at Berkeley, CA 94720-3112, USA
Correspondence to be sent to: Wayne M. Getz, Division of Insect Biology, ESPM, 201 Wellman Hall, University of California at Berkeley, CA 94720-3112, USA. e-mail:getz{at}nature.berkeley.edu
| Abstract |
|---|
|
|
|---|
A model is presented that predicts the instantaneous spike rate of an olfactory receptor neuron (ORN) in response to the quality and concentration of an odor stimulus. The model accounts for the chemical kinetics of ligandreceptor binding and activation processes, and implicitly the initiation of second messenger cascades that lead to depolarization and/or hyperpolarization of the ORN membrane. Both of these polarizing processes are included in the most general form of the model, as well as a process that restores the voltage to its negative resting state. The spike rate is assumed to be linearly proportional to the level of voltage depolarization above a critical negative voltage level. The model includes the simplifying assumption that activation of bound ligandreceptor complexes by G-proteins and other enabling molecules follows a Monod function that has the ratio of enabling molecules to bound unactivated ligandreceptor complexes as its argument. Parameters are selected that provide an excellent fit of the model to previously published empirical data on the response of cockroach ORNs to pulsed 1-hexanol stimuli. The sensitivity of model output to various model parameters is investigated and changes to parameters are discussed that would improve the ability of ORNs to follow rapidly pulsed stimuli.
| Introduction |
|---|
|
|
|---|
Olfactory receptor neurons (ORNs) transmit information to the brain of organisms on the presence and concentration of volatile chemicals in the environment. The ORNs are the gatekeepers of olfactory perception, both in the context of which chemicals are detectable and what resolution is possible with respect to temporal variability of odors in turbulent plumes or the concentration gradients of odor fields. The neural firing or spike rate response of an ORN to a temporally and spatially varying odor stimulus results from a concatenation of biochemical, molecular conformational and flux transport processes (Hildebrand and Shepherd, 1997
Olfactory perception in an organism begins with odorant molecules or odor ligands (Hildebrand and Shepherd, 1997
) penetrating the mucous layer of the
olfactory epithelium in vertebrates or entering through the pores of olfactory sensilla in insects,
diffusing or being transported across a liquid medium with the help of odorant binding proteins
or
OBPs; some of these OBPs may also be implicated in deactiving odorant molecules (Breer, 1994
; Kaissling, 1998b
), and then binding with
receptor molecules on the dendritic membranes of ORNs. The individual receptor molecules, of
which there are known to be hundreds or even thousands of types (Buck and Axel,
1991
; Lancet et al., 1993
; Reed, 1994
) belong to the seven transmembrane domain, G-protein coupled family of receptor
proteins. Some of these receptor molecules are known to be highly specific for particular
odorants, typically those implicated in transducing sex pheromone signals in insects (Kaissling, 1987
). Other receptor molecules, however, may well be generalists, able
to
bind to several different members of a group of ligands of the same or similar moieties (Singer and Shepherd, 1994
). In this case, we expect
ligandreceptor affinities or binding constants to vary among similar ligands (Lancet et al., 1993
).
Once a ligand binds to a receptor protein, the resulting ligandreceptor complex
couples with and activates G-proteins (Firestein and Zufall, 1994
). This
coupling in turn initiates a second messenger cascade that ultimately causes ion channels to open
and the membrane to depolarize or, in some cases, hyperpolarize, depending on which channels
have been opened (Ache, 1994
; Trotier, 1994
).
Changes in membrane polarization cause the ORNs to increase or decrease their rate of firing
with
respect to a cell specific background rate, depending on the concentration of the odor stimulus (Fujimura et al., 1991
; Akers and Getz, 1992
, 1993
; Getz and Akers, 1993
, 1997a
).
For all animals, stimulation or inhibition of the ORNs is the first, or sensory transduction,
phase of a multiphase perceptual process that has striking commonalities in all animals, with, of
course, some notable differences among taxons with vastly different brain morphologies. The
archetypal organisms for the model presented here are insects with highly developed olfactory
systems, such as cockroaches and honeybees. In insects, olfactory stimuli are typically
filamentous
plumes of odorant molecules that have a complex spatial and temporal structure (Moore and Atema, 1991
; Murlis et al., 1992
; Dittmer et al., 1995
; Vickers and Baker, 1994
). As the odorant molecules waft around the antennae of an individual organism, the
concentration of ligands available to the ORNs situated in olfactory sensilla of various types will
be highly stochastic and variable over time and at different locations on the same antenna at the
same point in time. The various processes altering the firing rates of the ORNs, as we show here,
smooth out some of the temporal jitter in the concentration, while variations in the concentration
experienced among ORNs are smoothed out by a large ensemble of ORNs converging at the next
level onto orders of magnitude fewer olfactory glomeruli (Shepherd, 1972; Rospars,
1988
; Boeckh et al., 1990
; Masson and
Mustaparta, 1990
; Smith and Getz,
1994
; Hildebrand and Shepherd, 1997
).
Mathematical models of ORNs (or, for that matter, taste receptor neurons) can focus either
on the kinetics of ligandreceptor binding (Ennis, 1989
, 1991
; Malaka et al., 1995
; Kaissling, 1998a
, 1998b
; Lánsky and
Rospars, 1999
) or on the alteration of receptor membrane potential and concomitant
generation of action potentials (Lánsky and Rospars, 1993
; Tuckwell et al., 1996
; Vermeulen et al.,
1996
; Vermeulen and Rospars, 1998
). In the former case, the
assumption can be made that spike rate is an appropriate function of the number of bound
ligandreceptors pairs (Getz and Akers, 1995
) or models
combining both aspects can be considered (Lánsky et al., 1994
; Rospars et al., 1996
). Here, following the ideas of
Lánsky and Rospars (Lánsky and Rospars, 1999
), I
develop a model that combines both ligandreceptor kinetics and membrane potential
dynamics. A central difficulty in this approach is to find a reasonable way to model the second
messenger cascade linking the ligandreceptor dynamics to changes in membrane
potential,
especially when some components of this cascade are not yet well understood. I circumvent this
difficulty by taking an approach inspired by consumer-resource concepts found in population
ecology (Getz, 1993
).
In this study, parameters in the model are selected to simulate observed data from cockroach
ORNs (Lemon and Getz, 1997
). Because the primary motivation for
developing the model is to use it to generate realistic spike rate input for a neural network model
of olfactory coding in the insect antennal lobes (Getz and Lutz, 1999
), I
focus on fitting the model to existing empirical data rather than capturing the details of all the
processes involved with generating the response. In some cases these details are not known,
particularly with regard to second messenger cascades. In other cases, the details are not
necessary to capture the essential dynamic response of ORNs, and the simplest reasonable
description is used; for example, I use one equation to model changes in the membrane voltage
potential rather than the usual four to capture the actual spike profiles (Mascagni and
Sherman, 1998
). Also, if ligandOBP dynamics occur on a much faster time
scale than spike generation, then replacement of a dynamical description of ligandOBP
interaction with an appropriate static (equilibrium) value that, perhaps, includes a small time
delay
may be a reasonable simplification to make. Thus, although the model does provide complete
insight into the processes that determine the response profiles of ORNs, the critical test of the
model is how well it can generate realistic response profiles to be used as input into a network
model of olfactory coding in the insect antennal lobes. As presented below, the model developed
here does provide a good fit to the response of cockroach ORNs (Lemon and Getz,
1997
).
| Materials and methods |
|---|
|
|
|---|
Modeling ligandreceptor binding kinetics
Chemical kinetic models of the rate at which populations of receptor molecules at density R on the membranes of ORNs bind to ligand messenger molecules typically begin with the assumption that monovalent ligands at density L in the perireceptor space (i.e. the space around the receptor neuron membrane where they are sufficiently proximate to interact with membrane receptor molecules) are competing for unbound receptors at density U. This implies that the density B of bound receptors is B = RU, and that the reaction is of the form
![]() |
where k1 and k1 are respectively the
association and disassociation rates of this reaction (Figure 1). Bound
receptors, however, do not initiate second messenger cascades until they have activated a coupled
G-protein (Firestein and Shepherd, 1991
). The rate at which bound
receptors at density B become activated receptors at density A is typically
modeled by the simple kinetic process
|
![]() |
(Lánsky and Rospars, 1999
). One can assume also that the
ligands themselves enter the perireceptor space at a rate determined by an externally determined
flux density k0Lin (Lin is a
density
and k0 is a rate) and leave at the same per-capita (i.e. per molecule) rate they
enter. Further, one can assume, as in Lánsky and Rospars (Lánsky and
Rospars, 1999
), that the ligands become inactivated during the dissociation process
![]() |
where Ld is the density of degraded ligands that are no longer able to
play a role in generating a response. Recent evidence, however, suggests that termination of
olfactory signaling has more to do with direct deactivation of ligand
receptorG-protein complexes, than ligand degradation per se (Breer,
1994
). For want of specifics regarding the mechanisms of receptor inactivation, I
assume that the activated ligandreceptorG-protein complexes are inactivated at a
per-unit rate k2, leaving behind, after a sequence of reactions much
faster than k2, an intact unbound receptor and a degraded ligand
(Figure 1). Thus the inactivation processes can be approximated by the
dissociation process
![]() |
leading to equations (1)(3) in the Appendix.
To model the activation process in a more realistic manner, one may need to take account of
the fact that the activation processes require the presence of resources in the form of ATP,
G-proteins, Ca2+ (Restrepo et al., 1996
) and other
molecules on which the activation rate k2 depends (Figure 1). One might expect the activation rate k2 to decline with
decreases in the density M of these second-messengerrelated molecules, which here I
generically refer to as `enabling molecules`. Thus, at high densities of L,
if
the complex B is produced rapidly, then the activation process may become partially
exhausted through a reduction in the density M, although M will probably not
decline to zero. The reason for this is because these enabling molecules, at density M,
are
maintained by cellular processes that replenish or recycle the molecules involved in the
activation
processes and restore the density of these molecules to a resting level M0
during a refractory period after the ligand input density Lin becomes zero
[in
reality, ligand input density is a function of time Lin(t)].
In the absence of specific details, one might expect the activation rate k2
to be an increasing function of the ratio of the density M of activation enabling
molecules
to the density B of bound receptors, i.e. k2 = k2(M/B). One would expect this function to saturate at some maximal rate k2* when (i.e. when enabling molecules are not limiting). A simple
function
with this saturating property is the Monod function, where the onset of saturation is controlled by
its associated MichaelisMenten `half-saturation` parameter M1/2see equation (4) in Appendix A; for a review of
this function see Smith and Waltman (Smith and Waltman, 1995
). If the
density M of these enabling molecules is drawn down at a rate proportional to the
activation process rate k2B, while simultaneously being restored by
cellular processes to a resting density M0, then the kinetics of M is
modeled by equation (5) in the Appendix.
Modeling neuron spike rates
The relationship between the density of activated ligand receptor complexes and the
spike rate of an olfactory neuron is complicated by the second messenger cascade that ultimately
results in the opening of the ion channels that lead to the production of spikes. Without the
details
of all, or even the major, biochemical pathways involved in this second messenger cascade, I
need
to make some simplifying assumptions regarding the relationship in questions. First, I assume
that
the adaptive processes going on in the cell are captured by the notion of enabling molecules and
the saturating kinetic relationship expressed in equation (4). Second, I
explicitly model the voltage dynamics of the receptor neuron membrane and then make the
standard assumption that spike rate is some appropriate function of membrane voltage (Rospars et al., 1996
).
Consider a membrane voltage depolarization and resting state restoration model based on the simple assumption that competing restorative and depolarizing forces are proportional to differences in membrane voltage respectively with regards to resting (Vrest) and maximum depolarizing (Vdep) levels (note Vrest < Vdep). Also, assume that the rate of membrane depolarization is proportional to the density of activated ligandreceptor complexes A(t). Under these assumptions, changes in membrane voltage V(t) are modeled by the equation (6) in the Appendix. Further, if it is now assumed that the maximum spike rate for a neuron is Smax, and that a neuron will not spike if V(t) < Vcrit, then the simplest model for the instantaneous spike rate incorporating this threshold is a clipped linear time-delay model scaled by the maximum spiking rate with changes in voltage normalized by the maximum possible change maximum equal to (Vdep Vcrit)see equation (7) in the Appendix.
It is not directly possible to measure an instantaneous spike rate, only an average spike rate
(t1,t2) over an arbitrary time
interval [t1,t2] or, of course, interspike intervals
[equation (8) in Appendix]. Here, to be compatible with data used to fit
the model,
we compute the average spike rates over consecutive 50 ms time bins.
| Results |
|---|
|
|
|---|
Baseline dynamics
I first set out to explore the dynamics of the ligandreceptor binding and activation process, modeled by equations (1)(5) in the Appendix. These equations constitute a fourth-order nonlinear system of ordinary differential equations in the `ligand', `ligandreceptor', `activated ligand receptor' and `enabling molecules' density variables L, B, A and M, respectively. These equations contain seven rate parameters (k0, k1, k1, k2*, k2, k3 and k3), an input parameter or, more generally, function of time Lin(t), three density scaling parameters, R, M1/2 and M0, and, of course, four initial conditions, L(0), B(0), A(0) and M(0). Without loss of generality, I scaled the units of all parameters with respect to the total number of receptor molecules R = U(t) + A(t) + B(t) per unit receptor neuron membrane. Note that we have assumed R to be constant over time, even though the relative densities of receptors in the different states A(t) and B(t) change with time. To implement this scaling, we set R = 1 and interpret the other densities in terms of R; for example, Lin = 0.1 means that the ligand density in the perireceptor space is 1/10th the density of the total receptor molecule density R. Without loss of generality, we can also scale time. For notational simplicity, the variable we choose as our basic unit of time is k2* = 1, the maximum rate at which ligandreceptor complexes are activated [which only occurs when densities of M(t) are much greater than M1/2Bsee expression (4)].
The model is capable of producing a complex array of output patterns for different values of the rate and scaling parameters, as well as initial conditions. One way of investigating the effects of these parameters and initial conditions is by systematically studying several simpler special cases and the sensitivity of solution output to these special cases. The particular cases selected here are motivated by the need both to understand the range of possible patterns (analysis of contrasting cases) and to search for output patterns that resemble empirical data. In particular, it would be useful to obtain an understanding of how certain processes influence the temporal profiles of the response of receptors to constant, as well as pulsatile, stimuli.
To begin, I conducted a baseline analysis for the simplified case k0
, k1 = 0, k1 = 1, k2 = 1, k3 = 1, k3 = 1 and Lin(t) = 1 for all t
0 (recall k2* = 1
and R = 1), by first numerically simulating the behavior of the solution to equations (1)(5) when the density scaling
parameters
have the values M1/2 = 1 and M0 = 1 (Figure 2 AD), or M0 = 10 (Figure 2
EH). The initial conditions used are the resting values of the system (i.e. the
equilibrium state of the system in the absence of any external stimuli); i.e. B(0) = 0, A(0) = 0, and M(0) = M0. I then considered the effects of
altering the rate parameters k1, k2, k3 and k3 with respect to the basal simulations of the
equations over the interval t
[0,20]. Note that letting k0
in equation (1) is equivalent to assuming that at all times,
the actual density L(t) is equal to the input density Lin(t). Relaxing this assumption basically sets up a time lag, inversely related to the size of
the
rate constant k0, in which the ligand density follows the input Lin(t) at a distance that depends on how rapidly ligands bind to available
receptors. The dynamics would be a little more complicated than just a pure time lag, however,
and the ORNs would then take on features associated with flux detectors, rather than
concentration detectors, as discussed by Kaissling (Kaissling, 1998a
).
|
The results indicate (Figure 2 A,E) that the rate k1 at which ligand and receptor molecules associate serves to scale the equilibrium proportion of activated ligandreceptor complexes up to a saturating proportion (the value for very large k1) of ~0.35 when M0 = 1 (Figure 2 A) and 0.5 when M0 = 10 (Figure 2 E), given that the other rate parameters are unity. On the other hand, the parameter k1 influences the shape of the transient rather than the equilibrium value when the enabling molecule density recovery rate is a relatively slow k3 = 1/30 (Figure 3 A). The rate k2 at which activated ligand receptor complexes dissociate into available unbound receptors and the density Ld of disabled ligands has an inverse effect on the ligandreceptor complex equilibrium proportion, both when k3 = 1 (Figure 2 B,F) and k3 = 1/30 (Figure 3 B). The enabling molecule recovery rate k3 produces a strong peak-shaped ligandreceptor activation profile when the other parameters have the basal values of unity, but then only when it has a sufficiently small (i.e. slow) value itself (Figure 3 C). This peak disappears when the enabling molecule resting density is increased by an order of magnitude to M0 = 10 (Figure 3 G) and the transient become insensitive to the enabling molecule recovery rate k3. The rate k3 at which at which the enabling molecules are used up when activating ligandreceptor complexes also has an increasing effect on the ligandreceptor complex equilibrium proportion, both when k3 = 1 (Figure 2 D) and k3 = 1/30 (Figure 3 C). At the higher resting density of ten enabling units per receptor (M0 = 10), the equilibrium proportion of activated ligandreceptor complexes hardly increases for corresponding values of the various rate parameters, although a strong peak arises when the rate k3 increases to 30 (Figure 2 H) or even when k3 is only equal to 1 for the case k3 = 1/30 (Figure 3 C).
|
Fitting real receptor neuron spike profiles
Empirical data from both the honey bee Apis mellifera (Akers and Getz,
1992
, 1993
; Getz and Akers, 1993
)
and the American cockroach Periplaneta americana (Lemon and Getz, 1997
) indicate that the average firing rate of a population of olfactory receptor neurons
has a peak-shaped profile similar to that displayed inFig 3A for the cases k1 = 1 and k1 = 30.
These data were obtained from recording the responses of a population of cockroach ORNs
to a pulsed 1.25 Hz square-wave 1-hexanol stimulus (Lemon and Getz, 1997
). The means and standard errors of these data (N = 20 receptors in the data
set) are illustrated inFigure 4.
|
Equations (1)(7) can be used to simulate the data inFigure 4. These seven equations contain one input function Lin(t) and 17 adjustable parameters (R, k0, k1, k1, k2*, k2, k3, k3, M1/2, M0, a0, a1, Vrest, Vdep, Vcrit, Smax and t). Two of these parameters can be set to unity to scale the equations to relative units of time and density. Thus I set k2* = 1 and R = 1. Also, as previously mentioned, I let k0
which allows the
ligand density to be replaced with an input function Lin(t) (note
that
I made this assumption more realistic by imposing a 20 ms time delay on this function). The
remaining 14 parameters represent a large number of degrees of freedom for fitting the output of
a model to the given data. I thus reduced the number of free parameters by further setting Vrest = 50 mV, which is within 10 mV of the resting potential of many
cockroach neurons (Burrows, 1996An informal search over the remaining parameter space produced the following set of parameters that provided an excellent fit to the empirical data (Figure 4): k1 = 5, k1 = 100, k2 = 2, k3 = 100, k3 = 3.5, M1/2 = 0.1, M0 = 10, a0 = 10, a1 = 80, Vcrit = 45 and t = 0.1 (i.e. 20 ms).
From visual inspection, these parameters provide a good fit. Of course, a better fit could be
obtained using least-squared or maximum-likelihood methods (Hilborn and Mangel,
1997
). This level of precision, however, would only be justified if more scenarios
were used to make the fitting procedure more robust (for example, several different
concentrations, input waveforms, etc.), an exercise beyond the scope of this paper and the
available data.
Finally, additional simulations were undertaken to explore the sensitivity of different features in the spike rate profile (for example, existence of a prominent peak, shape of peak, on- and off-set characteristics;Figure 5) and the ability of the spike train to follow a pulsed stimulus (Figure 6). The implications of these results are discussed below.
|
|
| Discussion |
|---|
|
|
|---|
Relationship among parameters
The primary purpose of this paper is to present the model and investigate how well the
processes incorporated in the model are able to reproduce empirically measured output, because
the model will then be used to generate input into a neural network model of the insect antennal
lobes (Getz and Lutz, 1999
). If the model fits a limited data set well, then
one can move to collecting a more demanding set of data (more concentrations, a greater variety
of stimulus patterns, etc.) on which to test the model and refine or extend it. Further, in testing a
model, one's confidence in how well the model captures the essential processes increases
with the number of parameters that are independently estimated; for example, direct
measurement
of Vrest, Vdep and Vcrit, and of
as
many of the rate constants and component densities as possible (Kaissling, 1998b
). Thus, the ultimate test of the model is to see how well it reproduces the output for
a variety of stimuli when all of the model parameters are estimated independently of the output.
Given the paucity of suitable data, I was unable to estimate (i.e. fix) values for the majority of the parameters before selecting the values for the remaining set. Further, the model certainly omits the details of many processes contributing to the production of ORN spike trains in response to a given stimulus. Even so, the model still has considerable value in generating hypotheses regarding relationships among the rate constants of the identified processes, and then testing these to see how well they stand up to scrutiny in real systems.
For example, in trying to find a set of parameters that provide an approximate fit to the data, it is clear fromFigures 2 and3, especiallyFigure 3 A, that one needs to select a set of parameters for which k1 << 1 and k3 >> 1. From a visual inspection of the fit between model output and empirical data (Figure 4), it is apparent that the model captures well the reduction in the second and third peaks of the data compared with the first peak. The model also captures the fast rise in the first peak, but clearly shows that the data indicate a slower rising second peak and an even slower rising third peak that could be due to relatively slow adaptive processes not accounted for in the model. Several important points emerge from a comparative study of how the spike rate changes with the relative values of the parameters used to fit the model (Figure 5).
First, a relatively large (fast) ligandreceptor binding rate, k1, is needed to keep the the observed peak (Figure 4) sharp. For example, if we drop the value of k1 by a factor of five to k1 = 1, no adaptation (i.e. presence of a peak) is evident: the response is cut off at its maximum value by the off-set of the stimulus (Figure 5A). On the other hand, once k1 is sufficiently large, increasing its value by a factor of five to k1 = 25 has a relatively small effect on the shape of the peak (Figure 5 A).
Second, the ligandreceptor dissociation rate, k1, can be relatively large (fast) compared with the ligandreceptor binding rate k1, and the transient peak still retained, but not too large to prevent significant ligandreceptor activation and initiation of a dense second messenger cascade (Figure 5 B). On the other hand, k1 should not be too small (slow), since this leads to a slower decay of spike activity once the stimulus has been removed (Figure 5 B).
Third, generation of a strong peak requires the rate k3 at which enabling molecules catalyze the activity of bound receptor/ligand complexes to be considerably faster than the enabling molecule cellular restoration rate k3 (cf.Figure 2 Cwith k3 = 1/30 andFigure 2 H with k3 = 30). This difference causes enabling molecules to become limiting and adaptation to set in, provided the absolute values of the rates k3 and k3 are appropriately scaled with respect to the resting enabling molecule density M0 and the ligand input density Lin (compareFigure 2 C with 2G and Figure 2 D with 2 H). Further, when the input ligand density is pulsed, if k3 is decreased to the point where it is almost two orders of magnitude smaller (slower) than k3, then the height of the first peak is hardly affected, but the second and third peaks are now a fraction of the height of the first peak (Figure 5 C). By contrast, the second and third peaks become much closer to the height of the first peak if either k3 is increased or k3 is decreased (Figure 5 C,D).
Fourth, the peaked behavior is not evident unless the ligand input and/or the enabling molecule resting concentrations, Lin and M0 respectively, are sufficiently high (Figure 5 E,F). These results, together with the above results, imply that the existence of a transient peak is dependent on the relative relationship of the densities of the variables (ligands, free receptors, enabling molecules) to the values of the rate parameters. Therefore, the absolute values of the rate parameters can only be set by measuring the true densities of the variables involved or vice versa. Note that the half saturation parameter M1/2 has little effect on the shape of the curve for the set of parameters considered here, so the results have not been included inFigure 5.
Fifth, the parameter Vcrit, as expected, affects the absolute height of the peak but, more importantly, it is critical in determining the offset characteristics of the spike train. If the critical voltage at which the cell begins its spike activity, Vcrit, equals the resting voltage Vrest (= 50 mV), then the spike activity follows the process of the restoration of the voltage to its resting state so that the relatively sharp switch of the neuron to its off position is lost (Figure 5 G). On the other hand, since our model is deterministic, no background spike rate is possible when the neuron is switched off. The low background spiking level evident in the empirical data can easily be incorporated into the model by including an appropriate low noise level in the spike generation process.
Finally, the voltage parameters a0 and a1 need to be relatively large for the spike train to reproduce (or follow) steep inclines and declines, as well as sharp peaks, that occur in the density of the activated ligandreceptor complexes (cf. Figures 2 and 3). These peaks, however, will not be effectively followed unless the voltage depolarizing parameter a1 is larger than the voltage restoration parameter a0 (Figure 5 H). Further, the larger this discrepancy, the more easily the receptor is able to rise to its maximum spike frequency of Smax = 200 Hz.
Pulsatile stimuli
Empirical data indicate that odor stimuli often take the form of turbulent odor plumes
containing discrete packets of odor. Thus, as an insect flies through or along an odor plume or an
odor plume is transported over a stationary insect, the packets of odor are perceived as a
sequence of pulses of varying pulse and inter-pulse duration (Murlis et al.,
1992
; Vickers and Baker, 1994
). If the inter-pulse intervals
are short compared with the constants characterizing the response of the ORNs, we can expect
the ORNs to smooth out the pulses, while if the inter-pulse intervals are relatively large, then
each
pulse should be transparent to the system (Lemon and Getz, 1997
). By
comparing model output to empirical data on the response of ORNs to pulsatile stimuli, one can
test whether the time constants in the model are compatible with observed data. The simulation
results presented in Figure 6 provide some insight into how model output
changes for changes in selected rate constants.
The specific set of ORN parameters used to fit the cockroach data presented inFigure 4 smooth out pulsatile square wave stimuli of amplitude 2 oscillating at frequencies of 40 Hz or higher (Figure 6 A). Once these pulsatile stimuli have dropped to frequencies of ~15 Hz, the oscillations are evident in the simulations (Figure 6 B), although in real ORNs noise may mask these oscillations in individual response profiles. At frequencies of ~5 Hz (Figure 6 C), the oscillations are strongly evident and unlikely to be masked by noise, even in individual ORNs.
Real populations of cockroach olfactory receptor neurons appear to be able to follow to some
extent oscillations of 20 Hz and even 40 Hz (Lemon and Getz, 1997
).
One of the reason our simulated receptors are not able to follow frequencies much higher than 10
Hz is in part due to the fact that the offset response (at the end of each square wave) is too slow.
The offset response rate can be increased by decreasing the value of Vcrit.
This, however, would lead to a significant decrease in the spike rate (Figure 5 G), which can be compensated by reducing the value of k3 and the ratio of a0 to a1. Further the absolute
values of both these parameters can be increased so that the voltage can more accurately follow
changes in the densitiy of the activated ligandreceptor complexes. Responses of ORNs
modeled by equations (1)(6), using the
standard parameter set used to generateFigure 6 AC, with the
modifications k3 = 10, a0 = 50, a1 = 1000 and Vcrit = 10, produces output that clearly
follow 4
and 15 Hz pulses (Figure 6 E,F). Populations of these ORNs would even
follow 40 Hz pulses in populations of neurons if the signal-to-noise ratio were not too large
(Figure 6 D).
Synergism and inhibition
Many insects and crustacean chemosensory receptors exhibit synergistic or inhibitory
characteristics when responding to mixtures of several odorants or stimulants (Daniel
and Derby, 1988
, 1991a
, 1991b
; Derby et al., 1991
; Fine-Levy and Derby, 1992
; Getz and Akers, 1995
). Inhibition is particularly widespread
in arthropods and is thought to be associated with IPs second messenger pathways (Ache and Zhainazarov, 1995
; Daniel et al., 1996
). The phenomena of synergism and inhibitions have been defined in several
different
ways, some more useful than others (Getz and Akers, 1995
).
Interactions between odorants in generating a response rate S in an olfactory receptor neuron can occur at several different levels. First, two ligands may compete directly for binding sites on the same receptor molecules. Second, different ligands may bind with different receptor molecules on the dendritic membrane of the same receptor neuron and these bound ligandreceptor complexes may be activated by the same set of enabling molecules; third, by a partially overlapping set of enabling molecules; or, fourth, by a non-overlapping set of enabling molecules. Fifth, even if two different ligands initiate distinct second messenger cascades, the final products of these cascades may open the same ion channels or, sixth, different ion channels in the receptor cell membrane.
In the last case, one ligand may lead to depolarization and the other to hyperpolarization of the receptor neuron membrane. The voltage model developed here is easily extended to include hyperpolarizing pathways. Specifically, the depolarizing and hyperpolarizing pathways can each be modeled by a set of equations [equations (1)(5) in the Appendix] with variables and parameters subscripted (or double subscripted if already subscripted) to indicate the pathway in question. Thus, in this case, for i = 1, 2, we can have Ai(Li,t) representing the density of activated ligand receptor complex at time t when the density of ligands i in the perireceptor space has the time profile Li(t) = Li. If we now assume that depolarizing, restorative and hyperpolarizing ionic channel and pump processes compete to drive the voltage to the respective levels Vdep > Vrest > Vhyp (note that the resting membrane potential Vrest is negative in neurons) at rates that are proportional to the voltage gradients involved and, additionally, the depolarizing and hyperpolarizing rates are respectively proportional to the densities of the activated ligandreceptor complexes A1(L1,t) and A2(L2,t), then the membrane voltage [equation (6)] can be extended to include both pathways [equation (9) in the Appendix]. Unfortunately, no appropriate data exist to explore how well this extension is able to simulate the response of olfactory receptor cells containing both types of pathways, so a test of the extended model remains a tantalizing opportunity for future studies.
Finally, from an olfactory coding perspective, it may be more important for the model to fit
data over some parts of the response profile than others. For example, honey bee ORN data (Getz and Akers, 1992
), which exhibit peak response rates occurring, as
with the cockroach, during the interval 50150 ms after the onset of each stimulus (Figure
4), appear to code more information during this peak response period
than during post-peak period of 150250ms (Getz and Akers, 1997b
). If this situation also applies to cockroach ORN data, then it is more critical that the
model fits the data over the interval 50150 ms after the onset of a stimulus than over
time
intervals much beyond this particular window of time.
| Conclusion |
|---|
|
|
|---|
Chemical signaling processes based on G-protein associated receptor transduction mechanisms incorporate complex second messenger cascades involving many tens of interactions and synthesis of intermediate products, not to mention links to more general cellular processes involving the cycling of such energy resources as ATP, and all-purpose messenger molecules such as cAMP and Ca2+. A precise model for the generation of membrane voltage profiles, and hence spiking behavior, of ORNs would require that all the essential details of the associated second messenger cascades be known. Thus a precise model is not yet feasible, since many details remain to be worked out. Further, as long as a complete description of all cellular processes competing for critical energy and messenger molecules is not included in the model, simplifying assumptions need to be made to obtain workable models of the response of ORNs to stimulation. These assumptions may be reasonable if the processes modeled by equations are orders of magnitude faster or larger than those subsumed under postulates of constant background levels, constant fluxes, or phenomenologically specified relationships, such as the one postulated in equation (4).
Most ORN response models in the literature have been either been applicable to cells that
respond to single odorants (such as pheromone receptors in insects), or to general olfactory cells
stimulated by odorants rather than complex odor stimuli. Some models do consider stimulation
by
complex odor stimuli but, as in Malaka et al. (Malaka et al., 1995
), they only consider the equilibrium situation, they do not model the membrane
voltage and they do not consider nonlinearities arising from the ligandreceptor activation
process being limited by ATP, G-proteins, or other enabling molecules ([f. equation (4)]. Also, many studies focus on the equilibrium (tonic) rather than transitory
(phasic)
properties of response neurons (Rospars et al., 1996
), despite the
fact that the phasic response in general olfactory receptors in insects lasts at least 200 ms (Lemon and Getz, 1997
). This phasic period is not much shorter than the
time it takes worker honey bees, for example, to discriminate different odors (Smith
and Menzel, 1989
).
Although many details of the pathway involved in transducing olfactory signals into electrical responses in insect and other animal olfactory sensory neurons are not known, and a number of simplifying assumptions were made in developing the model presented here, the model captures exceptionally well the spike rate of cockroach olfactory neurons to pulsatile odors. More detailed empirical data on the response of neurons to pure and mixed odor stimuli with various concentration profiles (for example, both square and triangular shaped pulses) are required to provide a richer array of patterns for fitting the model. Also, more details on the various transduction pathways are needed before we can assess how well a nonlinear relationship, such as portrayed in equation (4), is really able to encapsulate the complex array of processes associated with signal transduction by ORNs.
| Appendix |
|---|
|
|
|---|
Glossary of symbols
Independent variable
t time (scaled so that k 2* = 1)
System variables.
The dynamics of each variable is modeled by its own differential equation.
A density of activated ligandreceptor complexes
B density of `bound-but-not-yet-activated' ligandreceptor complexes
L density of free ligands available for binding
M density of `enabling molecules' facilitating the ligand receptor activation process
V membrane voltage of olfactory receptor neuron
Auxiliary variables
R density of membrane receptors normalized to 1 so that all related
densities are in `membrane-receptor-density' units
L in(t )ligand input flux as a function of time t
S (t )spike rate of olfactory receptor neuron as a function of time t
(t 1,t 2)average
spike rate over the interval [t 1,t 2]
U density of unbound receptors, U = R B A
Rate parameters.
a0rate at which membrane voltage decays back to resting state
a1rate at which membrane voltage depolarizes
a2rate at which membrane voltage hyperpolarizes (two pathway model only)
k0rate at which ligands enter the perireceptor space; if it is
assumed that k 0
, this is equivalent to assuming L (t ) = L in(t )
k0rate at which ligands exit the perireceptor space; it is assumed k0 = 0positive k0 is not needed to prevent the pooling of active ligands in the perireceptor space because, after dissociating from receptors, ligands are assumed to be transformed to a degraded state
k1rate at which ligand and receptors associate to form a bound-but-not-yet-activated complex
k1rate at which bound-but-not-yet-activated ligandreceptor complexes disassociate
k2*maximum rate at which bound-but-not-yet-activated ligand receptor complexes become activated; the actual activation rate k2 is a function of the enabling molecule density M [see equation (4) below]; note that time is scaled so that k2* = 1
k2rate at which activated ligandreceptor complexes disassociate to leave behind a degraded ligand molecule
k3rate at which enabling molecules are replenished when their density is very low
k3rate at which enabling molecules are used in ligandreceptor activation process
Scaling parameters.
M0resting or saturation density for enabling molecules
M1/2half-saturation or MichaelisMenten constant in the Monod activation rate function k2(M /B ) (see equation 4)
Smaxmaximum spiking rate (set to 200 spikes/s)
tspike generation time delay
Vrestcell membrane resting voltage (set to 50 mV)
Vdepmaximum voltage to which cell membrane can be depolarized (set to +50 mV)
Vcritdepolarization threshold for spiking to occur
Vhypminimum voltage to which cell membrane can be hyperpolarized
Model equations
From the processes described in the text, it follows that the kinetic equations for L , B and A are (Lánsky and Rospars, 1999
).
![]() | (1) |
![]() | (2) |
![]() | (3) |
For the case where the association rate k2 is a Monod function of the ratio M/B rather the density M itself we have:
![]() | (4) |
[The time argument t on the right-hand side is included to stress that k 2 varies with time as M and B vary with time. Note that ratio M/B does not appear explicitly on the right-hand side of equation (4) because we have multiplied the top and bottom by B].
From the assumptions in the text, the density M of these enabling molecules is modeled by the equation.
![]() | (5) |
The dynamics of the membrane voltage depends on the current voltage state and the number of activated receptors and is modeled by the equation
![]() | (6) |
A clipped linear time-delay relationship between voltage and spiking rate, as described in the text, is given by
![]() | (7) |
Hence the average spiking rates over intervals [t 1, t2] is
![]() | (8) |
In the case where two pathways exist, each satisfying a set of equations (1)(5) that determine the densities of corresponding activated ligandreceptor complexes A1(L1,t) (depolarizing pathway) and A2(L2,t) (hyperpolarizing pathway), then the membrane voltage will be governed by the equation.
![]() | (9) |
| Acknowledgments |
|---|
I thank William Lemon, Petr Lánsky and Jean-Pierre Rospars for comments on the manuscript. This work was support by NSF Grant IBN 98-07938 to W.M.G.
| References |
|---|
|
|
|---|
Ache, B.W. (1994) Towards a common strategy for transducing olfactory information. Cell Biol., 5, 5563.
Ache, B. W. and Zhainazarov, A. (1995) Dual second-messenger pathways in olfactory transduction. Curr. Opin. Neurobiol., 5, 461466.[Web of Science][Medline]
Akers, R.P. and Getz, W.M. (1992) A test
of identified response classes among olfactory receptors in the honeybee worker. Chem. Senses, 17, 191209.
Akers, R.P. and Getz, W.M. (1993) Response of olfactory sensory neurons in honey bees to odorants and their binary mixtures. J. Comp. Physiol., 173, 169185.
Boeckh, J., Distler, P., Ernst, K.D., Hösl, M. and Malun, D. (1990) Olfactory bulb and antennal lobe. In Schild, D. (ed.), Chemosensory Information Processing. Springer, Berlin, pp. 201227.
Breer, H. (1994) Odor recognition and second messenger signaling in olfactory receptor neurons. Cell Biol., 5, 2532.
Buck, L. and Axel, R. (1991) A novel multigene family may encode odorant receptors: a molecular basis for odor recognition.Cell, 65, 175187.[Web of Science][Medline]
Burrows, M. (1996) The Neurobiology of an Insect Brain. Oxford University Press, New York, pp. 682.
Daniel, P.C. and Derby, C.D. (1988) Behavioral olfactory discrimination of mixtures in the spiny lobster( Panulirus argus )based on a habituation paradigm. Chem. Senses, 13,385
395.
Daniel, P.C. and Derby, C.D. (1991a) Chemosensory responses to mixtures: a model based on compostion of receptor cell types. Physiol. Behav., 49, 581589.[Medline]
Daniel, P.C. and Derby, C.D. (1991b) Mixture suppression in behavior: the antennular flick response in the spiny lobster towards binary odorant mixtures. Physiol. Behav., 49, 591601.[Medline]
Daniel, P.C., Burgess, M.F. and Derby, C.D. (1996) Responses of olfactory receptor neurons in the spiny lobster to binary mixtures are predictable using a noncompetitive model that incorporates excitatory and inhibitory transduction pathways. J. Comp. Physiol A., 178, 523536.[Medline]
Derby, C.D., Girardot, M. and Daniel, P.C. (1991) Responses of olfactory receptor cells of spiny lobsters to binary mixtures. II. Pattern
mixture interactions. J. Neurophysiol., 66, 131139.
Dittmer, K., Grasso, F. and Atema, J. (1995) Effects of varying plume turbulence on temporal concentration signals available to orienting lobsters. Biol. Bull., 189, 232233.[Web of Science]
Ennis, D.M. (1989) A receptor model for binary
mixtures applied to the sweetness of fructose and glucose: De Graaf and Frijters revisited. Chem. Senses, 14, 597604.
Ennis, D.M. (1991) Molecular mixture models
based
on competitive and non-competitive agonism. Chem. Senses, 16, 117.
Fine-Levy, J.B. and Derby, C.D. (1992) Behavioral discrimination of binary mixtures and their components: effects of mixture
interactions
on coding of stimulus intensity and quality. Chem. Senses, 17, 307323.
Firestein, S. and Shepherd, G.M. (1991) A kinetic model of the odor response in single olfactory receptor neurons. J. Steroid Biochem. Mol. Biol., 39, 615620.[Web of Science][Medline]
Firestein, S. and Zufall, F. (1994) The cyclic nucleotide gated channel of olfactory receptor neurons. Cell Biol., 5, 3946.
Fujimura, K., Yokohari, F. and Tateda, H. (1991) Classification of antennal olfactory receptors of the cockroach, Periplaneta americana. L. Zool. Sci., 8, 243255.
Getz, W.M. (1993) Metaphysiological and evolutionary dynamics of populations exploiting constant and interactive resources: rK selection revisited. Evol. Ecol., 7, 287305.
Getz, W.M. and Akers, R.P. (1993) Olfactory response characteristics and tuning structure of placodes in the honey bee. Apis mellifera. L. Apidologie, 24, 195217.
Getz, W.M. and Akers, R.P. (1995) Partitioning nonlinearities in the response of honey bee olfactory receptor neurons to binary odors. BioSystems, 34, 2740.[Web of Science][Medline]
Getz, W.M. and Akers, R.P. (1997a) Response of American cockroach (Periplaneta americana) olfactory receptors to selected alcohol odorants and their binary combinations. J. Comp. Physiol. A, 180, 701709.
Getz, W.M. and Akers, R.P. (1997b) Coding properties of peak and average response rates in American cockroach olfactory cells.BioSystems , 40, 5563.[Web of Science][Medline]
Getz, W.M. and Lutz, A. (1999) A neural
network model of general olfactory coding in the insect antennal lobe. Chem. Senses, 24, 351372.
Hilborn, R. and Mangel, M. (1997) The Ecological Detective: Confronting Models with Data. Princeton University Press, Princeton, NJ.
Hildebrand, J.G. and Shepherd, G.M. (1997) Mechanisms of olfactory discrimination: converging evidence for common principles across phyla. Annu. Rev. Neurosci., 20, 595631.[Web of Science][Medline]
Kaissling, K.-E. (1987) In Colbow, K. (ed.), R.H. Wright Lectures on Insect Olfaction. Simon Fraser University, Burnaby, British Columbia, Canada.
Kaissling, K.-E. (1998a) Flux detectors versus
concentration detectors: two types of chemoreceptors. Chem. Senses, 23, 99111.
Kaissling, K.-E. (1998b) Pheromone deactivation
catalyzed by receptor molecules: a quantitative kinetic model. Chem. Senses, 23, 385395.
Koester, J. (1991) Membrane potential. In Kandel, E.R., Schwartz, J.H., Jessell, T.M. (eds), Principles of Neural Science, 3rd edn. Elsevier, New York, pp. 81103.
Lancet, D., Sadovsky, E. and Seidemann, E. (1993) Probability model for molecular recognition in biological receptor repertoires:
significance to the olfactory system. Proc. Natl Acad. Sci. USA, 90, 37153719.
Lánsky, P. and Rospars, J. (1993) Coding of odor intensity. Biosystems, 31, 1538.[Web of Science][Medline]
Lánsky, P. and Rospars, J. (1999) Odorant concentration and receptor potential in olfactory sensory neurons. Chem. Senses, in press.
Lánsky, P., Rospars, J.-P. and Vermeulen, A. (1994) Basic mechanisms of coding stimulus intensity in the olfactory sensory neuron. Neural Process. Lett., 1, 912.
Lemon, W.C. and Getz, W.M. (1997) Temporal resolution of general odor pulses by olfactory sensory neurons in American cockroaches. J. Exp. Biol., 200, 18091819.[Abstract]
Malaka, R., Ragg, T. and Hammer, M. (1995) Kinetic models of odor transduction implemented as artificial neural networkssimulations of complex response properties of honeybee olfactory neurons.Biol. Cybernet. , 73, 195207.[Web of Science][Medline]
Mascagni, M.V. and Sherman, A.S. (1998) Numerical methods for neuronal modeling. In Koch, C. and Segev, I. (eds), Methods in Neuronal Modeling: From Ions to Networks, 2nd edn. MIT Press, Cambridge, MA, pp. 569606.
Masson, C. and Mustaparta, H. (1990) Chemical information processing in the olfactory system of insects. Physiol. Rev., 70, 199245.
Moore, P.A. and Atema, J. (1991) Spatial information in the three-dimensional fine structure of an aquatic odor plume. Biol. Bull., 181, 408418.[Abstract]
Murlis, J., Elkington, J.S. and Cardè, R.T. (1992) Odor plumes and how insects use them. Annu. Rev. Entomol., 37, 50532.[Web of Science]
Reed, R.R. (1994) The molecular basis of sensitivity and specificity in olfaction. Cell Biol., 5, 3338.
Restrepo, D., Teeter, J.H. and Schild, D. (1996) Second messenger signaling in olfactory transduction. J. Neurobiol., 30, 3748.[Web of Science][Medline]
Rospars, J.P. (1988) Structure and development of the insect antenno-deutocerebral System. Int. J. Insect Morphol. Embryol., 17, 243294.
Rospars, J., Lánsky, P., Tuckwell, H.C. and Vermeulen, A. (1996) Coding of odor intensity in a steady-state deterministic model of an olfactory receptor neuron. J. Computat. Neurosci., 3, 5172.[Web of Science][Medline]
Singer, M.S. and Shepherd, G.M. (1994) Molecular modeling of ligand receptor interactions in the OR5 olfactory receptor. Neuroreport, 5, 12971300.[Web of Science][Medline]
Smith, B.H. and Getz, W.M. (1994) Nonpheromonal olfactory processing in insects. Annu. Rev. Entomol., 39, 351375.[Web of Science]
Smith, B.H. and Menzel, R. (1989) An analysis of variability in the feeding motor program of the honey bee; the role of learning in releasing a modal action pattern. Ethology, 82, 6981.
Smith, H. and Waltman, P. (1995) The Theory of the Chemostat. Cambridge University Press, Cambridge.
Trotier, D. (1994) Intensity coding in olfactory receptor cells. Cell Biol., 5, 4754.
Tuckwell, H.C., Rospars, J., Vermeulen, A. and Lánsky, P. (1996) Time-dependent solutions for a cable model of an olfactory receptor neuron. J. Theor. Biol., 181, 2531.[Web of Science][Medline]
Vermeulen, A. and Rospars, J. (1998) Dendritic integration in olfactory sensory neurons: a steady-state analysis of how the neuron structure and neuron environment influence the coding of odour intensity. J. Computat. Neurosci., 5, 243266.
Vermeulen, A., Rospars, J., Lánsky, P. and Tuckwell, H.C. (1996) Coding of stimulus intensity in an olfactory receptor neuron: role of neuron spatial extent and passive dendritic backpropagation of action potentials. Bull. Math. Biol., 58, 493512.[Web of Science][Medline]
Vickers, N.J. and Baker, T.C. (1994) Reiterative responses to single strands of odor promote sustained upwind flight and odor source
location by moths. Proc. Natl Acad. Sci. USA, 91, 57565760.
Accepted May 25, 1999
![]()
CiteULike
Connotea
Del.icio.us What's this?
This article has been cited by other articles:
![]() |
D. P. Dougherty, G. A. Wright, and A. C. Yew Computational model of the cAMP-mediated sensory response and calcium-dependent adaptation in vertebrate olfactory receptor neurons PNAS, July 26, 2005; 102(30): 10415 - 10420. [Abstract] [Full Text] [PDF] |
||||
![]() |
F. W. Grasso Invertebrate-Inspired Sensory-Motor Systems and Autonomous, Olfactory-Guided Exploration Biol. Bull., April 1, 2001; 200(2): 160 - 168. [Abstract] [Full Text] [PDF] |
||||
![]() |
W. M. Getz and P. Lansky Receptor Dissociation Constants and the Information Entropy of Membranes Coding Ligand Concentration Chem Senses, February 1, 2001; 26(2): 95 - 104. [Abstract] [Full Text] [PDF] |
||||
![]() |
J.-P. Rospars, V. Krivan, and P. Lansky Perireceptor and Receptor Events in Olfaction. Comparison of Concentration and Flux Detectors: a Modeling Study Chem Senses, June 1, 2000; 25(3): 293 - 311. [Abstract] [Full Text] [PDF] |
||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||





















