Computational Model of the Insect Pheromone Transduction Cascade

A biophysical model of receptor potential generation in the male moth olfactory receptor neuron is presented. It takes into account all pre-effector processes—the translocation of pheromone molecules from air to sensillum lymph, their deactivation and interaction with the receptors, and the G-protein and effector enzyme activation—and focuses on the main post-effector processes. These processes involve the production and degradation of second messengers (IP3 and DAG), the opening and closing of a series of ionic channels (IP3-gated Ca2+ channel, DAG-gated cationic channel, Ca2+-gated Cl− channel, and Ca2+- and voltage-gated K+ channel), and Ca2+ extrusion mechanisms. The whole network is regulated by modulators (protein kinase C and Ca2+-calmodulin) that exert feedback inhibition on the effector and channels. The evolution in time of these linked chemical species and currents and the resulting membrane potentials in response to single pulse stimulation of various intensities were simulated. The unknown parameter values were fitted by comparison to the amplitude and temporal characteristics (rising and falling times) of the experimentally measured receptor potential at various pheromone doses. The model obtained captures the main features of the dose–response curves: the wide dynamic range of six decades with the same amplitudes as the experimental data, the short rising time, and the long falling time. It also reproduces the second messenger kinetics. It suggests that the two main types of depolarizing ionic channels play different roles at low and high pheromone concentrations; the DAG-gated cationic channel plays the major role for depolarization at low concentrations, and the Ca2+-gated Cl− channel plays the major role for depolarization at middle and high concentrations. Several testable predictions are proposed, and future developments are discussed.


Introduction
Olfactory receptor neurons (ORNs) are essential for the recognition of odor molecules. In vertebrates this recognition is performed by several hundreds olfactory receptor proteins (ORs) borne by the ORN plasma membrane, each ORN expressing a single type of receptor [1]. In insects a smaller number of ORs have been identified [2][3][4]. In male moths, ORNs housed in antennal sensilla trichodea ( Figure 1) can detect female-released sexual pheromone with exquisite sensitivity, specificity and efficiency [5]. These ORNs have been the subject of intensive studies during the last fifty years using molecular, radiochemical, pharmacological, electrophysiological, calcium imaging, behavioral and modeling techniques (reviewed in [6][7][8]). The latter contribution has been significant and ORNs have experienced a rich history of modeling, since reports that a male moth can find a pheromone releasing female from several miles away [9,10] and that a single pheromone molecule is sufficient to elicit an action potential in the moth sensory neurons [11]. The system has been modeled at the level of behavior [12,13], at the level of antenna as biomechanical filter for odor molecules [14][15][16], at the level of electrical circuits that give rise to action potentials (e.g. [17,18]), and at the level of biochemical processes that lead to neuronal activation [19][20][21][22]. The most detailed model yet published is that of Kaissling [23] which attempted to account for the production of the ''receptor potential'' through the interactions of a process generally referred to as ''perireceptor events''. In fact this process consists of a biochemical network of the carrier proteins (pheromone binding proteins, PBPs), ORs and odor degrading enzymes [8,24] which occupy a common space surrounding the outer dendritic receptive membrane of ORNs.
A reasonably complete picture of the transduction processes emerges from these studies, although some of the molecular and ionic channel mechanisms underlying the transduction process still remain elusive. After perireceptor processes, the pheromone bound ORs are believed to interact with a G-protein which in turn activates the effector enzyme phospholipase C-b (PLCb) [25]. This enzyme catalyzes the production of second messenger molecules, inositol 1,4,5-triphosphate (IP 3 ) and diacylglycerol (DAG), which trigger the opening of a cascade of various ionic channels. The resulting ionic currents generate the receptor potential (RP) which passively propagates to the ORN soma and axon where it generates action potentials. Recently, this classical metabotropic mechanism has been challenged in insect ORNs and a direct coupling of the OR to a cationic channel has been proposed in parallel or in replacement [26][27][28]. These new developments are important from molecular, physiological and evolutionary points of view.
The full description of such a complex signaling network, involving both feedforward and feedback processes, is a daunting task. Modeling can contribute to this description by integrating various effects and displaying quantitatively what results from the interplay of all molecular actors. The knowledge accumulated on the pheromonal ORN is sufficient to start building a model of its transduction cascade, and to test whether it can effectively link together some of the known facts and suggest new experiments. Thus, the first aim of our investigation was to develop a qualitative model of the pheromone transduction cascade integrating the known molecular and ionic mechanisms. The second aim was to translate these mechanisms, wherever possible, into a set of differential equations and to determine the quantitative values of their parameters. We made a systematic search of known values and determined the unknown values by fitting the model output to the properties of the experimentally measured RP. These properties were systematically determined in response to ''square'' pulses of pheromone of constant duration at several intensities [23,29,30]. They offer the most precise data on the transduction cascade available so far. These responses are characterized by a rapid rising phase, a slow falling phase, especially at high concentrations, and an extremely wide dynamic range of about 6 decades from threshold to saturation.
In insects, most modeling efforts have been dedicated to the perireceptor and receptor processes in moth pheromone sensilla [23,[31][32][33][34][35]. Although interactions of ORs, G-proteins and effectors have been recently studied [36], no model has been proposed yet for post-effector processes in insects. The model we present here focuses on these processes in male moth pheromone ORNs and takes advantage of the modeling studies available on olfactory transduction in vertebrates [37][38][39][40][41][42].
Beyond fitting adequately the experimental dose-response curves we addressed the following related questions. What are the functional roles of the various currents? In particular, what could be the respective roles of the direct (ionotropic) and indirect (metabotropic) gating mechanisms of the initial cationic current? What are the mechanisms behind the characteristics of the concentration-response curves (broad dynamic range, short rising time and long falling time)? What are the processes that contribute most to the amplifying function of the cascade?

Results
In the first three subsections a formal model of pheromone transduction is presented. In the next three subsections the model is fitted to experimental data and its properties are studied.

Qualitative Model of Pheromone Transduction
Based on experimental results obtained in moth ORNs, complemented when necessary with data coming from other animal species and some reasonable assumptions, we developed a global qualitative model of pheromone transduction. A schematic diagram of the model is shown in Figures 2 and 3. This model is summarized in this section. Some of the experimental results and the main assumptions (denoted A to F) on which it rests, are briefly mentioned and listed in Table 1. Complementary justifications, references and comments are provided in the Discussion section.
Perireception (Figure 2, steps 1 and 2). After adsorption on the cuticle the pheromone molecules enter the hair lumen through micropores in the sensillum wall. Within the aqueous sensillum lymph that fills the lumen, they bind to PBPs which carry them to the ORs borne by the ORN outer dendritic membrane. They are also degraded by enzymes [21,43]. These processes can be fundamentally viewed as two competing effects, one which is the entrance of molecules from the outside, corresponding to an uptake measured in micromole of pheromone per liter per second,

Author Summary
All sensory neurons transduce their natural stimulus, whether a molecule, a photon, or a mechanical force, in an electrical current flowing through their sensory membrane via similar molecular and ionic mechanisms. Olfactory receptor neurons (ORNs), whose stimuli are volatile molecules, are no exception, including one of the best known: the exquisitely sensitive ORNs of male moths that detect the sexual pheromone released by conspecific females. We provide a detailed computational model of the intracellular molecular mechanisms at work in this ORN type. We describe qualitatively and quantitatively how the initial event, the interaction of pheromone molecules with specialized receptors at the ORN surface, is amplified through a sequence of linked biochemical and electrical events into a whole cell response, the receptor potential. We detail the respective roles of the upward activating reactions involving a cascade of ionic channels permeable to cations, chloride and potassium, their control by feedback inactivating mechanisms, and the central regulatory role of calcium. This computational model contributes to an integrated understanding of this signalling pathway, provides testable hypotheses, and suggests new experimental approaches. Figure 1. Moth pheromone-sensitive sensillum trichodeum in tip-recording conditions. The sensillum is a small organ typically composed of 2 ORNs and 3 auxiliary cells (thecogen Th, trichogen Tr and Tormogen To), housed within a porous cuticular hair. The tight junctions between cells separate the ORN extracellular environment in two parts with different ionic compositions, the sensillar lymph bathing the outer dendritic segment (sensory) and the hemolymph bathing the inner dendrite and soma. In experimental conditions the pheromone is delivered close to the hair. The ORN electrical response is recorded extracellularly with an electrode slipped on the cut hair tip. Figures 2  and 3 give detailed views of the ORN membrane processes at the molecular level. Figure 6 gives an overview of the global electrical organization of the sensillum. Modified from [8]. doi:10.1371/journal.pcbi.1000321.g001 and the other which is the degradation or hypothetic deactivation [23] of pheromone molecules. When the system is stimulated by a square wave of pheromone, all pheromone molecules are not immediately removed so their concentration grows until there is an exact balance between uptake and removal. When stimulation ends, uptake returns to zero but removal continues until all pheromone molecules are removed and their concentration quickly falls to zero. This system, called flux detector by Kaissling [31], would not work without removal because the pheromone molecules are trapped inside the perireceptor space.
Reception and the two coupling mechanisms (Figures 2 and 3, steps 3 to 6 and 9b). The pheromone molecule binds to an OR (step 3) then activates it (step 4), which presumably corresponds to a conformational change of the OR. In our model, the pheromone-activated OR (R*) can follow two possible pathways (assumption A). In the first pathway, it binds a Gprotein to give an activated G-protein (G*, step 5) which itself combines with an effector enzyme, PLCb [25], to produce an activated effector (denoted E*, step 6). The G-protein is involved in a loop which returns it to its initial state and the cycle can start again. The three proteins, R, G and E, can encounter one another and interact because they diffuse in the membrane. Moreover, each activated OR can activate several G-proteins when it diffuses and so contributes to signal amplification. The name of ''random walk amplifier'' [44] was given to this concept. In the second pathway (step 9b, top left of Figure 3), R* directly binds a cationic channel [26][27][28].
Second-messenger production (Figures 2 and 3, step   7). The effector enzyme catalyzes the cleavage of phosphatidylinositol biphosphate (PIP 2 ) producing IP 3 and DAG. The respective roles of IP 3 and DAG have not yet been completely clarified [45][46][47][48]. Generation of IP 3 induced by pheromones was found to be species-and tissue-specific; it occurs only in male antennae [49,50]. The involvement of this enzyme in insect ORN responses was demonstrated by the fast and transient production of IP 3 after incubation of moth antennal homogenates with pheromone compounds [49,51,52] as well as with non-pheromonal odorants in locust and cockroach [25,50]. Its implication has also been demonstrated by a genetic approach in Drosophila [53]. Upon application of pheromone, the concentration of IP 3 increases very rapidly reaching a maximum after about 50 ms, declines quickly to a lower plateau, then declines further with a slower time course to the basal level within a few hundred ms [25]. The production of IP 3 is GTP-dependent [25,49].
Opening of cationic channels (Figure 3, steps 9 and 9b). First, DAG activates a non-specific cationic channel (step 9). These DAG-gated cationic channels were observed in vivo from outer dendritic segments in A. polyphemus [46] and in cultured ORNs of Spodoptera littoralis [48]. Also, the perfusion of sensilla trichodea with DAG increases the firing activity of ORNs in A. polyphemus and B. mori [47,55]. Second, ORNs express an unusual member of the insect OR family, known as OR83b in D. melanogaster and also found in several moth species [56]. The coexpression of OR83b with conventional ORs is necessary to get odor-evoked responses both in vivo and in vitro [26,57]. Both proteins interact with one another to form a heteromeric receptor complex [58]. OR83b, alone [27] or heteromerized with the OR [28], was recently identified as a cationic channel. So, two cationic channels are apparently involved, one which can be directly activated by R* (step 9b) and the other by DAG (step 9).
Opening and closing of Ca 2+ -dependent Cl 2 channels ( Figure 3, step 13). An increase in intracellular Ca 2+ activates Cl 2 currents in moth ORNs [59][60][61]. We hypothesized that, as in vertebrates [62], the Cl 2 current is depolarizing in insect ORNs (assumption C). No experimental evidence of the indirect inhibition of this Cl 2 current by Ca 2+ (for example via activated protein kinase C, PKC*) was found in S. littoralis [61]. For this reason we did not include any feedback regulation of this current in the basic model. However, since this is the only current without feedback, we examined a variant where it is inhibited by PKC*, in agreement with experimental data in Xenopus oocytes [63]. . Ca 2+ also binds to a complex of protein kinase C (PKC) and DAG. The resulting activated complex PKC* (steps 14 and 15) can phosphorylate PLC (step 16) which down-regulates its activity. In antennal homogenates from A. polyphemus, pheromone stimulation induces a 6-fold increase in PKC activity [55]. First, PKC, possibly activated by DAG and intracellular Ca 2+ [64,65], appears to be involved in the termination of the pheromone-dependent rise of IP 3 since PKC inhibitors prolonged the pheromone-induced transient IP 3 rise [66,67]. Second (step 16), in many systems, a PKC-dependent feedback regulation of PLCb has been observed [68]. In S. littoralis it has been shown that antennal PLCb has PKC binding sites (Chouquet et al., in preparation). Other effects of PKC* are described below. Ca 2+ extrusion ( Figure 3, step 17). Ca 2+ must be extruded from the ORN after stimulation. In frog [69] and squid [70] ORNs, as well as in other cell types [71], the presence of a Na + / Ca 2+ exchanger (NCX) has been demonstrated. Another extrusion mechanism found in other cell types involves a PMCA (plasma membrane ATPase pump). In insect ORNs the mechanisms of Ca 2+ extrusion are not known, which led us to compare the NCX and PMCA mechanisms.
Opening of Ca 2+ -dependent K + channels ( Figure 3, step 18). Intracellular Ca 2+ combined with depolarization activates K + channels. The largest current in ORNs of the moths Mamestra brassicae and S. littoralis is a voltage-gated and Ca 2+ -activated current [72,73]. This is a fast activating and sustained current with an outward rectification; K + flows out resulting in membrane repolarization. The conductance of the K + channels is 66 pS in M. sexta [74] and 180 pS in L. migratoria [75]. The location of these channels is unknown. In the model we assumed they are on the inner dendritic segment and soma (assumption D) because their repolarizing role is incompatible with the K + concentrations on both sides of the outer dendritic membrane (see below, paragraph ''Equilibrium and resting potentials'').

Biophysical Model
The qualitative description above, although indispensable, is not sufficient to gain a proper understanding of pheromone transduction. We must now turn to a formal description of the various steps involved. Note that abbreviations in roman (e.g. G*, E*, IP 3 etc.) denote chemical species, whereas the corresponding symbols in italics (e.g. G*, E*, IP 3 etc.) denote concentrations.
Pre-effector steps. A formal description of the perireceptor and receptor stage (steps 1 and 2, [23,34]) and the RGE stage (steps 3 to 6, [36]) were given previously and will not be repeated here. Briefly, it gives for a square pulse of pheromone of any duration and intensity, expressed in concentration L air (molarity in air) or, better, in uptake U (mole per liter per second; U = k i L air , with k i = 10 4 s 21 in the experimental conditions considered here), the concentration of the activated effector E* (and other intermediate species, including R* and G* for example) as a function of time. This system involves 13 chemical species and 12 reactions. It is described by a set of 13 ordinary differential equations and 4 conservation equations involving 17 parameters (4 initial protein concentrations, 10 reaction rate constants and 3 reaction rate constants limited by diffusion) which are given as equations (12)-(28) in the Methods section. Although very simplified, this model gives the same time-course of activated receptors as the more realistic model [23] and, likely, as the latest development of this model (Kaissling, manuscript in preparation). The rest of this section is devoted to a formal description of the post-effector network of reactions involving diffusible modulators as well as ionic channels from which the evolution of the membrane potential can be derived.
Diffusible modulators. The post-effector biochemical reactions involve five modulators IP 3 , DAG, Ca 2+ , CaCaM and PKC* (the latter results from the association of PKC, DAG and Ca 2+ ). These reactions are depicted schematically in Figure 3 and represented in standard biochemical notation in Figure 4.
The cleavage of PIP 2 by activated effector enzyme E * producing IP 3 and DAG is inhibited by PKC*. This is the only feedbackregulated reaction of the RGE stage in the model. The rate of production v of IP 3 and DAG was modeled by the following equation: where the variables are E*(t), the concentration of activated effector enzyme at time t, and PKC*(t), the concentration of activated PKC at time t. In the absence of PKC*, the reaction rate is maximal, v = s M E * , where s M is the maximal (uninhibited) production rate. In the presence of PKC*, the other constant parameters are K is , the concentration of PKC* needed for half-maximal inhibition, i.e. v = s M E * /2, and n is , the Hill coefficient of the inhibitory PKC*-E* reaction. The most important single modulator is Ca 2+ which acts as a second messenger to open Cl 2 and K + channels and acts as an inhibitor of PLC, IP 3 -and DAG-gated channels through the CaM and PKC pathways. We considered all three initial reactants, PIP 2 , CaM and PKC, as external species, i.e. available in unlimited quantity. All reactions were modeled as standard bidirectional reactions, with a forward production and a backward degradation. Their expression as a set of first-order differential equations is straightforward, see equations (29)- (34) in Methods.
Ionic channels. The ionic currents can be classified according to their gating mechanisms (molecule, ion and/or voltage) and ion permeability. The same formal description was applied to all of them, except for the OR83b cationic channel. All channels have an agonist Y, which triggers their activation, and some of them have an antagonist Z, which mediates their feedback inhibition. For example Y is DAG and Z is CaCaM for the cationic DAG-gated channel. When the concentration Y of the agonist increases the conductance of the channels G j for ionic current j (which can be Ca 2+ , cations, Cl 2 and K 2 ) increases according to a sigmoid Hill function ( Figure 5): where G Mj is the maximum ionic conductance of the channels, K j is the concentration of Y producing their half-maximal conductance, and n j is the Hill coefficient of the agonist-channel interaction. The antagonist moves this curve to the right, i.e. decreases its sensitivity by acting on K j (competitive inhibition, assumption E). This action involves another Hill function: In the absence of antagonist (Z = 0), the EC 50 is minimum (K j = K mj ). In the presence of Z, the EC 50 increases (K j .K mj ) and the concentration-response curve G j (Y) is shifted to the right (see example in Figure 5A). Knowing the channel conductance G j , the equilibrium potential of the permeating ion E j and the membrane potential V, the corresponding electrical current I j is given by Ohm's law: where DV is the potential difference between both sides of the membrane. This is in agreement with the linearity of the experimentally measured I-V curves of the unspecific cationic [76] and Ca 2+ -gated Cl 2 [62] channels. These equations where applied to describe the three currents having an agonist and an antagonist, i.e. Ca 2+ current I Ca (agonist IP 3 , antagonist CaCaM), cationic current I cat (agonist DAG, antagonist CaCaM), Cl 2 current I Cl (agonist Ca 2+ ). For the Ca 2+ extrusion current I x , in the absence of experimental data in insect ORNs, we considered both the PCMA and NCX hypotheses. The Ca 2+ current I x driven by the PCMA does not depend on the membrane potential. It is given by where E x is the maximal pump capacity and G x is the conductance of the pump. On the contrary, the net current I x driven by the NCX depends on the membrane potential according to eq. (4) [77]. In both the PMCA and the NCX we assumed that the variable conductance G x is given by eq. (2), with Ca 2+ as agonist Y and no antagonist Z, so that K x is a constant (assumption F).
Modified versions of these equations were used for the K + and leak currents. For the K + current I K (agonist Ca 2+ , no antagonist), which is also voltage-dependent, we used a modified version of the non-inactivating Ca 2+ where A K is a constant and the variables are the membrane potential (V) and the Ca 2+ concentration (Ca). Finally, the conductance G ld of the leak current I ld is also a constant given by the inverse of the membrane specific resistance at rest. The complete set of functions describing all 6 currents is given as equations (35)- (49) in Methods.
The OR83b cationic channel was not introduced in the present quantitative model because no formal description of its gating and regulating mechanisms is presently available. To our knowledge, no similar channel has been described in other neurons, which prevented extrapolation from known examples. The consequences of this approximation are examined in the Discussion section.
Receptor potential. If all channels were located in a patch of outer dendritic membrane and if this patch could be considered in isolation, the dynamics of the membrane potential V, defined as the difference of potential between inside and outside (taken as zero), would be given by: where C md is the capacitance of the membrane, I dep is the depolarizing current and I rep is the repolarizing current However, this simple model is not applicable for two reasons. First, the difference of potential experimentally recorded is between the sensillar lymph, bathing the outer dendrite, and the hemolymph (reference electrode), bathing the inner dendrite and soma ( Figure 6). These two media are separated by auxiliary cells which introduce a supplementary potential -the transepithelial potential. Second, the K + channel is located in the inner dendritic segment, bathed by the hemolymph, which constitutes another compartment. Therefore, a three-compartment model distinguishing outer dendrite, inner dendrite and soma, and auxiliary cell, is needed for an adequate description of the system. Besides the potentials, leak and ionic currents described above, this introduces three new potentials (outside the outer dendrite V ed , inside the inner dendrite and soma V is , and outside the auxiliary cell V ea ) and four new currents (from outer dendrite to soma I i , leak at soma I ls , through auxiliary cell I a and along sensillar lymph I e ). The functions giving these four currents are given in equations (50)-(53) and the set of differential equations linking potentials to currents is given in equations (54)-(57) which generalize equation (6) (see Methods). Potential V ed given by eq. (55) is the most important in practice because it corresponds to the difference of potential between the recording electrode, in contact with the sensillar lymph, and the reference electrode, in contact with the hemolymph. Potentials were obtained by numerical integration of equations (54)- (57). Finally the RP was calculated as the difference of potentials between the two sides of the outer dendritic membrane DV = V id 2V ed during stimulation and at rest However, as shown in Figure 6, the experimentally known potential is V ed , not DV. So, we computed its difference during stimulation and at rest, the so-called sensillar potential SP (SP is nearly proportional to RP and often called ''receptor potential'' in the literature).

Initial Values of the Variables and Values of Parameters
The post-effector model described above includes 65 values: the initial values of the concentrations of modulators and the membrane potential (10 values given in Table 2) and 55 parameter values. These parameters fall in 3 different categories: (1) the rate constants of the reactions involving the modulators (10 values); (2) the characteristics of the 5 main currents with 4 to 7 parameters per current (28 values); (3) the 17 parameters describing the dendritic morphology (surface and volume), the electrical properties of the dendritic membrane and the conversion factors from currents to ionic fluxes.
However, from a practical point of view, the most important distinction is between parameters with a known value (27) and those which were unknown (38). We considered as known, and used without change, any parameter value determined in a moth ORN, especially in Antheraea polyphemus. If a parameter value was known in a non-moth species, especially in the frog ( [40] and references therein), or as a ubiquitous component presents in any cell, we used it as a starting value. The fixed parameters based on experimental determinations or calculated from experimental data are given in Table 3. The 38 fitted parameters are given in Table 4 for modulators and Table 5 for ionic currents. Their final values were obtained as explained below and in the Methods section by comparison with experimentally known response characteristics.
Conversion factors. Some molecules (receptor, G protein, effector enzyme and DAG) are membrane bound. Their density, in molecules/mm 2 , can be expressed in intracellular concentration, in mM, using the following conversion factor where N A is Avogadro's number; S cut and V cut are the lateral area and volume respectively of the outer dendrite after cutting the hair tip (see Table 3). A similar formula was used for the conversion of extracellular concentrations (with same value of S cut but V cut replaced with the volume of the sensillum, see Table 2 in [36]). Ca 2+ appears in biochemical equations as a modulator and in equations of electrical currents as a permeable ion. The latter equations describe in electrical units the movement of Ca 2+ ions through the IP 3 -gated channels, the DAG-gated channels and the Ca 2+ extrusion exchangers. The relationship between current I (in pA, i.e. pC/s) and the chemical flux J (in mM/s) is J = fI, where the conversion factor f (expressed in mM/pC) is given by f = 1/zFV cut where z is the charge of the Ca 2+ ion, F the Faraday's constant (96484610 6 pC?mmole 21 ) and V cut the volume of the external dendrite (see Table 3). Different conversion factors, f Ca , f cat , and f x were applied to the three currents. Factor f Ca converts the inward IP 3 -gated Ca 2+ current into a Ca 2+ flux and takes also into account the buffering capacity of the intracellular medium. This In this and other tables (and in text) a dual notation is used: roman (e.g. PKC) for a species and italic (e.g. PKC) for its concentration. 2 For all other variables (L 0 , LN 0 , RL 0 , R * 0 , G * 0 , G b0 , G r0 , E * 0 , IP 30 , DAG 0 , Ca 0 , CaCaM 0 , PKCDAG 0 , PKC * 0 ) the initial values were taken as zero. doi:10.1371/journal.pcbi.1000321.t002 is necessary because a large proportion of the free Ca 2+ entering the cell (95 to 99% in adrenergic neurons; [79]) is rapidly bound to various molecules. Factor f cat converts the DAG-activated cationic current into a flux, takes into account the buffering effect and the fact that only a fraction of the cationic current is carried by Ca 2+ ions. Factor f x depends on the detailed mechanism of extrusion. For the PCMA we took f x = f. The NCX removes one Ca 2+ ion for 3 Na + ions. This produces a net transfer charge of one positive charge which contributes to the depolarization. Therefore the conversion factor for the Ca 2+ flux through the NCX is f x = 2f.
Equilibrium and resting potentials. The intracellular Ca 2+ concentration is constrained to be smaller than 0.01 mM because Ca 2+ -dependent channels, which are closed at rest, start to be activated at this concentration [61]. The concentration of the other permeating ions is not precisely known in the outer dendrite. Na + concentration was estimated at <1 mM and K + concentration at <150 mM in a moth ORN [74]. For Cl 2 no estimate was found in insects (in vertebrates its concentrations is <14 mM [80]).
The extracellular concentrations of these ions in the sensillum lymph bathing the outer dendrite are different from their  concentrations in the hemolymph. This difference creates the transepithelial potential mentioned previously. For Na + it is <25 mM and for K + it is <200 mM [17,81]. For Ca 2+ it is in the range 1 mM [17] to 6 mM [73]. For Cl 2 it is unknown. Reversal potentials were calculated from these concentrations: E K is close to zero, E Na is about 83 mV and E Ca is about 140 mV. The reversal potentials of Cl 2 and of the Ca 2+ extrusion mechanism being unknown, E Cl and E x were fitted (see below). As mentioned before, it follows from these values that the K + channels, experimentally known to be repolarizing, cannot be located in the outer dendrite, where the equilibrium potential of K + ions (<0 mV) is too high for such a role. However, the extracellular concentrations of K + in the hemolymph bathing the inner dendrite, soma and axon is 3.1 mM [82] or 20 mM [83], corresponding to an equilibrium potential E K <2100 or 250 mV, compatible with its expected role. This is the only ionic concentration in hemolymph used in the present work. To avoid extraneous complications in modeling repolarisation, we took E K equal to the resting potential.
The resting potential measured in vitro is <262 mV [72,74]. In vivo, the contribution of the transepithelial potential must also be taken into account. It is estimated at <35 mV [81]. As a result, the difference of potential at rest between the intracellular compartment of the outer dendrite and the sensillar lymph is <297 mV.

Comparison of Simulations with Experimental Measurements
Given the initial concentrations (Table 2) and the fixed parameter values (Table 3), computer simulations of the model were carried out. We searched for values of the unknown parameters, listed in Tables 4 and 5, yielding responses in accordance with experimental observations. Three sources of information were used. First, parameter values must remain in their physiological range. Second, the kinetic features of the second messengers and ionic currents must reproduce qualitatively the experimental observations. Third, the time evolution of the SP at various pheromone concentrations must agree quantitatively with the in vivo measurements of the SP performed in A. polyphemus [23,30].
Following a 2-s square pheromone pulse, the SP grows to a maximum then returns progressively to zero ( Figure 7H). This simple kinetics can be summarized with three numbers (Figure 8D), its height, its rising time, measured by the time it takes to reach half maximum, and its falling time, measured by the time to fall    [23,30] from the end of the stimulation to half-maximum. These three quantities depend on the dose of pheromone delivered to the system, measured either as a concentration in air in micromoles per liter or better as an uptake in mM per second. Zack [30] systematically determined the amplitude, rising time and falling time of the SP at various uptakes, from threshold to saturation. On this range, the amplitude increases, the rising time decreases 10 times and the falling time increases 10 times. These three doseresponse functions (Figure 8, dotted lines) were our main criteria for the fine tuning of the parameter values because they are the only ORN responses experimentally measured with precision in vivo.
We modified the unknown constants in the model to fit these experimental curves by using the same stimulation conditions as used by Zack [30]. We found a set of parameter values in agreement with known facts and giving good fits (see Methods), provided the mechanism of Ca 2+ extrusion depends on membrane potential. With the pump mechanism (PMCA), the falling time could be fitted only on a restricted range; for example if correct at low uptakes, it was much too large at higher uptakes. Therefore, all following results are given for the potential-dependent extrusion only (NCX). The values of the fitted parameters are given in Tables 4 and 5. The corresponding simulated dose-response curves are illustrated in Figure 8 (solid red lines) for the SP and in Figure 9 for the other variables (chemical species and currents).
Dose-response characteristic functions of the SP. Figure 8 shows that the simulated SP reproduces adequately the experimentally measured SP. First, it has the same wide pheromone stimulation range, from 10 24.75 to 10 1.5 mM/s ( Figure 8A solid red line). Second, it has short rising times ( Figure 8B

Main Properties of the Model
With the parameter values at hand the main properties of the model can be described. Some aspects require specific attention: Figure 9. Comparison of dose-response curves for height of the chemical species and currents. Height (left column; see definition in Figure 8D) and relative height (right column) for chemical species (top row) and for depolarizing and repolarizing currents (bottom row). (A) Ca 2+ (solid blue line) is the most abundant species (concentration divided 10 fold to be shown on the same scale as other species). (B) Responsiveness of all chemical species is much smaller than that of SP (curves shifted to the right of the SP curve shown as a solid red line), larger than that of effector enzyme E * (dashed red line) at low doses and smaller than E* at high doses. (C) Cl 2 (solid green), K + (dashed red) and cationic (dashed blue) currents are the most intense currents. (D) Responsiveness of the Cl 2 and cationic currents is higher than that of the effector enzyme (cf. (B)) and the IP 3 -gated Ca 2+ current (solid blue). In particular, the cationic current curve is close to that of SP (solid red curve, same as in B) and K + (confounded with SP) at all doses, while the curve of the Cl 2 current (solid green) is on the right of the SP curve at low doses (smaller responsiveness) and close to it at high doses. doi:10.1371/journal.pcbi.1000321.g009 the kinetics of currents and chemical species, their relative importance, their relative responsiveness and the explanation of the dose-response curves.
Kinetics of currents and chemical species. The kinetics of activated receptor, effector enzyme, second messengers, main depolarizing and repolarizing currents, and receptor potential are shown in Figure 7. The concentration of IP 3 (and of DAG; they are nearly identical due to the cleavage of PIP 2 ) increases very rapidly and then declines quickly to a lower plateau in response to middle and high uptakes ( Figure 7C). Very different types of kinetics are found. Some are phasic, others are tonic. For example, DAG production ( Figure 7C) and DAG-activated cationic current I cat ( Figure 7E) are phasic, whereas Ca 2+ increase ( Figure 7D) and Ca 2+ -activated Cl 2 current I Cl ( Figure 7F) are tonic. The phasic kinetics of I cat results from a quick and strong inhibition by Ca 2+ via CaCaM ( Figure 5A) whereas the tonic kinetics of I Cl results from the absence of inhibition. Interestingly, when inhibition of I Cl by Ca 2+ via PKC* was added to the model, a good fit with experimental SP data was obtained only for a very weak inhibition ( Figure 5B).

Relative importance of currents and chemical
species. The concentrations of the diffusible species and the intensities of the various currents are very different from one another. Their relative importance can be estimated on plots of height versus uptake ( Figure 9A and 9C). Calcium is the most abundant species at all uptakes ( Figure 9A). The DAG-activated cationic current I cat is the major depolarizing current at low uptakes, and the Ca 2+ -activated Cl 2 current I Cl is the major depolarizing current at medium and high uptakes. Their activation at different uptakes results ultimately from the values of their EC 50 s (see legend of Figure 5). The IP 3 -activated Ca 2+ current I Ca is a minor current at all uptakes: its maximum conductance is small (0.13 nS, Table 5), i.e. 10 to 60 times smaller than those of I cat and I Cl .
Relative responsiveness of the chemical species and currents are apparent on plots of relative amplitude versus uptake. They present significant differences. All species have a relatively low responsiveness, similar to that of the activated effector PLC, much lower than that of the SP (not shown). The EC 50 s of all species (ca. 0.5 mM/s) are much lower than that of the SP (less than 0.1 mM/s) ( Figure 9B). This difference in responsiveness between species and SP can be observed at all uptakes, although it decreases at higher uptakes. The ratios of the diffusible species concentrations with respect to E* are .1 at low uptakes and ,1 at high uptakes, which indicates that the dominant effect is amplification at low uptakes and inhibition of the secondmessenger production at high uptakes. This dual effect is well illustrated by the curve of the IP 3 -gated calcium current ( Figure 9D, solid blue line) which is close to the PLC curve ( Figure 9B, dashed red) and presents the same type of responsiveness as the diffusible species. On the contrary, the cationic and Cl 2 currents have a similar responsiveness to that of SP. Figure 9D shows that I cat is the most sensitive current. These results mean that the DAG-gated cationic current plays a major role in depolarization at low uptakes and a minor role at high uptakes.
Relative contribution of currents to the SP. To analyze more precisely the relative importance of the major depolarizing (cationic and Cl 2 ) and repolarizing (K + ) currents in the generation of the SP we selected four typical uptakes at regular intervals from low to high. The kinetic of the absolute values of these currents were compared in picoampere ( Figure 10) and after normalization with respect to their maxima ( Figure 11). Figure 10 shows that the cationic current is the most important in both amplitude and duration at low uptakes, the Cl 2 current takes over the dominant role at medium and high uptakes. The curve of the K + current is close to that of the cationic current at low uptakes ( Figure 10A) and close to the curve of Cl 2 current at high uptakes ( Figure 10C and 10D). However, the cationic current rises faster than both the Cl 2 current and the K + current (insets of Figures 10 and 11) at all uptakes. This means that the rapid rise (short half-rise time) of the SP should be attributed to the initial depolarization induced by the cationic current.
Phase space behavior of the modulators presents noteworthy properties. In the model, DAG and Ca 2+ are the two major modulators as they gate the two main depolarizing currents and activate the feedback inhibitors PKC* and CaCaM. In addition, Ca 2+ enters mainly through DAG-gated channel (the ratio of f Ca I Ca /f ct I cat very rarely exceeds 5% and its mean across all times and doses is 1.3%). The phase portraits in the DAG-Ca 2+ plane ( Figure 12) and the E * -SP plane ( Figure 13) show how the relations between the concentrations of DAG and Ca 2+ , and between E * and SP evolve in time at different uptakes. Let's consider first the DAG-Ca 2+ relationship. At low uptakes ( Figure 12A and 12B Figure 12C) the activation and inactivation parts start to separate, indicating a more complicated relationship. Finally ( Figure 12D and 12E) there is an almost complete separation, the phase portrait taking a characteristic L-shape: high concentrations of DAG (up to 1.4 mM) are associated with small concentrations of Ca 2+ (less than 20 mM) whereas high concentrations of Ca 2+ (up to 90 mM) are associated with small concentrations of DAG (less than 0.2 mM). The trajectory in the phase plane shows also that the rising speed is higher than the falling speed, particularly at high uptakes, as indicated by the short times to reach the maxima of DAG and Ca 2+ concentrations (these times are given in Figure 12). The same description holds true for the phase portrait in the E * -SP plane (Figure 13). At high uptakes it takes a characteristic upside-down L-shape.
Amplification mechanisms. A major property of the cascade is to amplify a relatively weak input signal into a strong output. This overall property can be partly quantified by the ratio SP r /E * r (subscript ''r'' stands for ''relative'') of the normalized SP (output) to the normalized concentration of activated effector enzyme (input). The normalization is necessary because the two quantities are not expressed in the same units. With this definition the total amplification of the cascade depends on the pheromone uptake: it is large at low uptake then progressively decreases at higher uptakes. A noteworthy consequence of the large amplification factor at low uptakes is the leftward extension of the dynamic range. Table 6 gives the contributions of each step at uptakes 10 24 , 10 23 and 10 22 mM/s respectively. DAG-gated and Ca 2+ -gated channels amplify the signal in different ways and they dominate the depolarization at different uptakes. At low uptakes, DAG-gated channels amplify the signal with a short rising time and the high amplification factor obtained (5300 at 10 24 mM/s) makes extremely weak signals detectable. At higher uptakes, IP 3and DAG-gated currents become transient and mainly work to provide the Ca 2+ entry and quick initial depolarization, and the Ca 2+ -gated Cl 2 current takes over as the dominant depolarization current.
Transduction delay. In order to estimate the relative contribution of extra-and intra-cellular reactions in SP rising and falling times, we stimulated directly the modeled cascade by 2-s square pulses of R* from 0.1 to 10 3 molecule/mm 2 , instead of a square pulse of pheromone, so removing the time taken by perireception and reception processes. We found that the rising time of E* decreases from 41 to 8 ms and that of SP decreases from 32 to 18 ms. (not shown). The falling time of E* increases from 42 ms to 102 ms and that of SP increases monotonically from 0.15 s to a maximum of 3.05 s. Since the experimentally measured values of the rising time decrease from 400 to 47 ms and that of the falling time increase from 1.5 s to 17 s on the same range, this means that most of the rising (92% to 66%) and falling times (90% to 82%) result from the extracellular reactions. Moreover, as far as the intracellular reactions are concerned, the contribution of the pre-effector steps to the falling time (28% to 3%) is much smaller than that of the post-effector steps, mostly because of the slow return of the intracellular Ca 2+ concentration at its resting level, particularly at high uptakes.

Sensitivity Analysis of Model Parameters
The sensitivity of the system to the parameters controlling each biochemical and electrical step was analyzed as explained in the Methods section. The effects produced on the SP responses by a change in the value of a single parameter at a time were examined. The main results of this analysis can be summarized as follows (see Text S1). (i) The sensitivity of SP to the parameters depends on the characteristic (height, rising or falling time) and on the dose. (ii) Each parameter has its greatest influence on one of the characteristics (Tables 4 and 5, rightmost column). The most influenced characteristic is usually the falling time (63%). (iv) The 12 most influential parameters are K is , G Mx , K mcat , K k , K icat , K pd , n cat , n Cl , G MCa , n x , G Mcat and K iCl . They are thus the best determined parameters. The 7 least sensitive parameters are s M , k cc2 , n iCa , i MCa , A k , K mCl and i Mcat . Tables in Text S1 list the parameters which most and least influence each characteristic.
Finally, we determined the importance of the feedback controls on second-messenger production, main ionic channels and calcium extrusion, by removing them one at a time (see Text S1). We compared the action of PKC* on PLC for different types of activation and found that it is inhibitory when activated by Ca 2+ only, but not when activated by DAG only.

Discussion
In this work, we propose a detailed model of the biochemical and electrical processes generating the receptor potential in the . DAG-gated current I cat (dashed blue) is the main depolarizing current at low dose (A). Ca 2+ -gated current I Cl (solid green) takes over the major role at high doses (C and D). The kinetic response of the repolarizing current I K (dashed red) is close to that of I cat at low dose (A) and close to that of I Cl at high doses (C and D). As shown in the insets, I K closely follows I cat at the beginning of the rising phase. doi:10.1371/journal.pcbi.1000321.g010 moth pheromone ORN. The role of the various currents and the molecular mechanisms of transduction are discussed in the first two subsections and the ORN response characteristics and the integrated cellular functions in the next two subsections. In the present state of knowledge, uncertainties remain on several of these processes. Therefore, the model studied is clearly incomplete and its features are not all equally well established. We discuss these uncertainties, qualitatively in the subsections ''Cationic currents'' and ''Calcium, chloride and potassium currents'' and quantitatively in the subsection ''Validity of parameter values''. However, the model helps to interpret ORN properties and suggests new experiments, as discussed in the last two subsections.

Cationic Currents
In the model the DAG-gated cationic current is the first ionic current in the transduction cascade (we attribute a different function to the parallel IP 3 -gated current, see below). A widespread assumption (our assumption A) is that this channel is indirectly coupled to the ORs via metabotropic pathway involving G-proteins and second messengers. This is based on the fact that like vertebrate ORs, insect ORs belong to the G-protein coupled receptor (GPCR) superfamily characterized by the canonical 7 transmembrane topology of its members [84]. However, recent bioinformatics and experimental investigations have revealed that the membrane topology of at least some Drosophila ORs differs from other GPCRs with an intracellular N-terminus and an extracellular C-terminus [58,85]. The structural distinction between insect and mammalian ORs, with different membrane topology, put into question the coupling of insect ORs to G proteins [86,87]. Recently, two studies indicated that, in expression systems, odorants can activate insect ORs and generate sensory currents independently of known G protein-coupled second messenger pathways, through a so-called ionotropic pathway [28], which may involve the cationic channel OR83b [27]. Proteins ortholog to OR83b are also expressed in moth pheromonal ORNs [56], which indicates that the ionotropic pathway may also be present in this neuron type.
However, several experimental observations made both in vivo and in vitro provide strong support for the involvement of G q proteins and PLC in the moth pheromone transduction cascade. First, G q proteins are present in moth antennae [88,89] and these proteins were localized in the outer dendrite of ORNs [54]. Second, the G protein activator, NaF, activates the firing activity of ORNs both in vivo [89] and in vitro (Lucas, unpublished results), demonstrating that G q proteins are functional in moth pheromone-responding ORNs. Third, Xenopus oocytes transfected with B. mori conventional pheromone receptors and without OR83b do not respond to pheromone stimuli unless they are co-transfected with G q proteins [90]. Fourth, pheromone stimulation activates the PLC activity as indicated by IP 3 production [25,49,50,52].
Two other arguments, obtained from non-pheromonal ORNs, strengthen the observations above in pheromonal ORNs. First, in addition to the fast and transient ionotropic pathway, at low concentration odorants also activate G proteins and the production of second messengers in expression systems [27]. This metabotropic current develops after a longer latency and with a slower kinetics but is more sensitive to odorants than the ionotropic pathway. Second, strong genetic evidence supports a role for a Gq-and PLCb-mediated signaling cascade during olfactory transduction in Drosophila [53].
The results at hand suggest the coexistence of two signaling pathways in both pheromonal and non-pheromonal ORNs, one ionotropic, the other metabotropic. Remarkably, both pathways gate a cationic channel, although it is not known whether the cationic channels of both pathways are the same (OR83b) or not, and if different whether their conductances are the same. The ionotropic signal, not amplified, is rapid and transient, and the metabotropic, amplified, is sensitive and prolonged, the former being more visible at high odorant concentration and the latter at low concentration. So, the relative intensities and timings of the two currents might be significantly different. If these properties are confirmed in future studies, they would suggest that the intensity of the cationic current found in the present model should not be entirely attributed to the DAG-gated channel because part of it comes from the receptor-gated channel. An advantage of the ionotropic pathway, especially at high uptake when the number of activated receptor is the largest and the ionotropic current is expected to become significant, is that its energetic cost for the neuron is lower. According to this view one might expect that the G-protein pathway is more inhibited at high uptake than found in the present model. The suppressed current would be compensated by the energy-saving direct coupling mechanism and the global effect would be the same. At the present time, the relative contribution of the two channels cannot be specified, if only because the reaction rate of R* with OR83b is not known. A second advantage of the OR83b cationic current is its time course. The receptor-gated current is expected to appear first and therefore to trigger an action potential faster than the DAG-gated current. The ionotropic pathway could contribute to explain the high speed of response of ORNs, especially at high pheromone concentration, and the very fast behavioral response of moths, down to 150 ms, both to contact and loss of pheromone filaments during their oriented flight to calling females [91]. The latency of the initial response was not considered in the present study and deserves more attention in the future.

Calcium, Chloride, and Potassium Currents
It remains uncertain whether the role of the Cl 2 current is depolarizing or repolarizing, because the intracellular and sensillar concentrations of Cl 2 are not known. We assumed that this current is depolarizing (assumption C) for three reasons. First, in vertebrates, the Cl 2 equilibrium potential is more positive than the resting potential [92] and a similar Ca 2+ -activated Cl 2 current amplifies ORN depolarization ( [93][94][95][96], see review [97]). Second, insect PBPs present in the mM range provide organic anions to the sensillum lymph [8] which decreases the sensillar Cl 2 concentration due to the principle of charge neutrality. Third, pheromone responses (SP and action potential firing) recorded in vivo were significantly higher when the sensillar Cl 2 concentration was lowered from 215 to 18 mM (Lucas, unpublished results).
One of the main uncertainties on channels (and other proteins) concerns their exact spatial location, which entails uncertainty on the mechanisms of Ca 2+ entry. Because the present model depends only on time, any diffusion or translocation is expressed in temporal (not spatial) terms and reflected in the reaction constants. This simplification has generally no incidence, except for two channel types. The first one relates to IP 3 -gated channels which have been located in the outer dendritic membrane, based on electrophysiology [45] and immunocytochemistry [54]. This location implies an inward flow of Ca 2+ from the sensillum lymph. However, besides or in replacement of this dendritic membrane channel, channels located in the endoplasmic reticulum may be considered. In this case Ca 2+ would come from intracellular stores. However this hypothesis is weakened by the fact that no intracellular Ca 2+ stores have been found in the outer dendrite.
In fact the present model is compatible with both possibilities and cannot discriminate them. A third possibility is discussed in the next paragraph.
The second channel for which the spatial localization is important is the repolarizing K + channel. A priori it can be located either in the outer dendrite or in the inner dendrite and soma region (as we did, our assumption D). The main objection against its outer dendritic location is that the equilibrium potential there is close to 0 mV, so that K + cannot have a repolarizing function. Against the inner dendritic location one can mention their modulation by Ca 2+ [72,74,75] because Ca 2+ would have to diffuse from the outer to the inner membrane to trigger them which is unlikely due to the poor diffusive ability of this ion. A possible solution to this problem is that IP 3 -gated channels are in  the inner dendrite. The IP 3 synthesized in the outer dendrite would have to diffuse to the inner dendrite and soma, which is compatible with the diffusive property of this second messenger, and, there, it would gate the Ca 2+ channel, possibly from Ca 2+ stores. Ca 2+ and voltage would then trigger the repolarizing K + current with a delay due to IP 3 diffusion.
This distinguishes the functions of these two different channel types in the model. IP 3 -gated channels let Ca 2+ ions flow into the cytoplasm, whereas DAG-gated channels also depolarize the ORN at low uptakes (see below). In fact it can be shown in the model that the IP 3 pathway is unnecessary for the depolarization.
On the basis of available experimental data we have included a feedback regulation of IP 3 channels (via CaCaM, our assumption B), but not of Cl 2 channels (via PKC * ). However, in the latter case, we studied the hypothesis of a feedback inhibition via PKC* [63]. We found that both regulations are not essential ( Figure 5B). For reducing the computation time we ignored the known feedback regulations on the pre-effector steps [66,98,99]. They are likely important for repetitive pheromone pulses that occur in natural conditions. In single-pulse condition, as in this study, this simplification presents no inconvenience provided the inhibitory regulation of, say, the receptors is not much faster than that of the effector enzymes.
All actors involved in the model have been shown to exist in the ORN membrane. The only exception is the Ca 2+ extrusion pathway which has not yet been described in moth ORNs, although it plays an essential role in the falling phase of SP. Because of the poor intracellular diffusion of Ca 2+ ions, it is likely that these mechanisms are located in the outer dendrite. We have shown that only a voltage-dependent extrusion is compatible with the experimental data, which strengthens the hypothesis of a Na + -Ca 2+ exchange (NCX), without ruling out the presence of an ATPase pump (PMCA). We found that the NCX pump needs no negative feedback control and that its reversal potential is E x <217 mV. Knowing the reversal potential of Ca 2+ (E Ca <140 mV) and the relationship between E x , E Ca and E Na , see eq. (46) in Methods, this value of E x implies E Na <88 mV, which is a reasonable value. Further experimental and theoretical investigations are needed to clarify the Ca 2+ extrusion mechanisms.
For simplicity we have not included the Ca 2+ -gated Ca 2+permeable cationic channels described in M. sexta ORNs [100] because with Ca 2+ entering the ORN and at the same time gating the channel, the control of the falling phase of the SP is made extremely difficult. Moreover, Ca 2+ -activated cationic channel could not be found in S. littoralis ORNs [61]. Perhaps qualitative differences exist in ion channel expression in ORNs across moth species providing different voltage-and time-dependent downregulation mechanisms.
Finally, when stimulated repetitively or for a longer time the ORN adapts and its response characteristics are different from those analyzed here [30]. Adaptation is beyond the scope of the present work as it may involve reactions other than those built in the present model.

Validity of Parameter Values
With the parameter values given in Tables 4 and 5 the model obtained accounts for several experimentally known properties of pheromone transduction in moths: The predicted SP reaches a maximum value and follows a time course, both in its rising and falling phases, that quantitatively agrees with the measured characteristics (depolarization, rising time and falling time) of the SP as a function of pheromone uptake (Figure 8).
The model also accounts quantitatively for the transient course of IP 3 production, as described in stop-flow experiments [25], with a very rapid increase followed by a quick decline at middle and high uptakes ( Figure 7C). It follows that the kinetics of DAG and the currents gated by both IP 3 and DAG must also be transient in this range of uptakes, which is the case in the model ( Figure 7E). The transient course of second messengers at middle and high uptakes results from a strong feedback inhibition. At low uptake the inhibition is much weaker and thus the course is not transient.
On the contrary, the intracellular Ca 2+ concentration ( Figure 7D) and the Ca 2+ -gated Cl 2 current ( Figure 7F) become sustained in the range of uptakes where the second messengers and their gated currents become transient: they increase at a slow rate then decrease gradually. These features are in agreement with experimental findings. Moreover, this Cl 2 current (I Cl ) is the major component of the depolarizing currents at middle and high uptakes, which qualitatively agrees with experimental data in the frog [101].
The phase portraits E*-SP ( Figure 13) show that at most uptakes except the highest, the maximum of E* and SP are reached at the same time and the trajectories of the rising and falling phases are close. This is equivalent to the fact, described by [23] that the concentration R 1 * giving a certain value of SP during the rising phase and the concentration R 2 * giving the same value of SP during the falling phase are equal, except at high uptake.
The experimental variability of SP is the highest at high uptake [30]. This can be interpreted, in the framework of our model, because we have shown that relatively small changes in some parameter values in different ORNs can lead to relatively large SP changes in this range of uptakes (Table 5).
However, several limitations affect the determination of the parameter values. Besides qualitative limitation regarding the completeness of the model, two other kinds of limitations must be taken into account. First, the model is based on ordinary differential equations depending only on time which entails two limitations: space is neglected, as discussed above, and very small concentrations are not adequately described because all chemical species must be in sufficiently large number to be considered as continuous variables. The total number of activated receptors per ORN is 30 when the uptake is 10 22.5 mM/s [36]. The bottleneck of the whole cascade being the receptors, all uptakes greater than this can be considered as adequately described by ordinary differential equations. Presently, for uptakes less than 10 22.5 mM/s, only mean values are obtained. A complete description in this range will require a stochastic approach, at least at the receptor level.
Second, supposing the model qualitatively valid for the stimulus used, the problem of the precision of the parameter estimation arises. The sensitivity analysis we performed permitted a classification in two categories (see tables in Text S1): the parameters whose modification changes significantly the SP response (right columns), which were therefore estimated with good precision, and those which do not influence much the SP (left columns), which are less well estimated. For example the equilibrium potentials need not be precisely known. This analysis one parameter at a time gives only a partial view because some parameters are linked. For example, the concentration E* and the maximum synthetic rate of s M of activated effector appear as a product, see equations (29) and (30) in Methods, so they cannot be known independently. The value given for s M is valid under the assumption that E 0 * = 0.136 mM [36]. However, it must be realized that the problem of parameter estimation is extremely constrained because the range of acceptable values of most parameters is restricted and because of the many feedback reactions. As a consequence the important parameters are not the same at all uptakes, especially for falling time. For this reason, finding a solution that works at all uptakes proved very difficult and suggests that a significantly different set of parameters in agreement with the experimental data available will not be easy to find.

Interpretation of the Global ORN Properties
The present model helps to interpret global properties of the ORN, especially its performance at efficiently encoding the stimulus. The response characteristics (amplitude, rising and falling times) of the SP in the pheromonal ORN present three remarkable properties: a wide dynamic range, a short rising time which decreases with pheromone concentration, and a long falling time which increases with concentration. The proposed model explains all these features. They can be analyzed from two different points of view: the relative contribution of the pre-and post-effector steps and the mechanisms by which the post-effector steps contribute to the observed features. (1) First of all, the posteffector cascade contributes considerably to the large dynamic range which extends over about 6 decades from 10 24.75 to 10 1.5 mM/s. Indeed, the dynamic range of the effector response is 3.25 decades only. The post-effector extension from 3.25 to 6 decades results from the collaboration of the two main currents, the DAG-gated cationic current and the Ca 2+ -gated Cl 2 current, which present a large difference in responsiveness. The cationic current has an (uninhibited) efficient concentration EC 50 of 0.01 mM of DAG, which is reached at a pheromone uptake of 10 24.25 mM/s ( Figure 5A), whereas the Cl 2 current has an EC 50 of 81 mM of Ca 2+ , which would be reached at an uptake of 50 mM/s ( Figure 5B), a value not actually reachable because the perireception system saturates at <30 mM/s. This means that the cationic current is most active at low uptakes whereas the Cl 2 current is most active at high uptakes. They complement one another and contribute, by the separation of their EC 50 s, to widen the dynamic range of the ORN. (2) The fast rising time of the SP results from the cationic current because at all uptakes it is faster than the Cl 2 current (compare the dashed and solid curves of insets in Figure 10). (3) The long falling phase of the SP is also explained by the two main currents. It reflects primarily the time course of the cationic current at low uptakes and the time course of the Cl 2 current (which closely follows the time variation of intracellular Ca 2+ ) at high uptakes. The long persistence of Ca 2+ at high uptakes suggests that Ca 2+ extrusion is not fast enough. However, as shown by the direct stimulation of receptors in the model, the bottleneck which limits the speed of rise and fall of SP is in the extracellular processes. Although involved and with many steps, the intracellular processes, as modeled here, are very fast, which confirms Kaissling's analyses [23].
An important function of the cascade is to transform a weak initial signal (pheromone binding to OR) into a strong local signal (RP and its corresponding SP). As already shown in a previous work [36] the amplification provided by the pre-effector stage, as quantified by the ratio E * r /R * r , is relatively small at any uptake (always less than 7.5). Consequently most of the amplification is provided by the post-effector stage, as quantified by SP r /E * r ( Table 6). The main amplification is in the electrical stages (I cat and I Cl ), since the amplifications of all chemical stages are relatively small. Cationic and Cl 2 channels amplify the signal in different ways and they dominate the depolarization in response to different uptakes. Cationic channels amplify the signal at low uptakes with a short rising time. Cl 2 channels amplify the signal by a larger maximal conductance and a longer duration of depolarization. The high amplification factor at low uptakes involves several mechanisms acting inversely on the activation and inactivation processes. Moreover these mechanisms are not the same at different uptakes. At extremely low stimulation uptakes, inactivation processes have weaker effects. This can be partly explained by the value of Hill coefficients which is #1 for n Ca , n cat , n x , and n icat but .1 for n is and n iCa . This means that the inhibition of the production of second messenger (IP 3 and DAG) and of the IP 3 -gated conductance develops on a relatively narrow range of concentration of their modulator, i.e. changes relatively abruptly in a threshold-like manner. On the contrary, the activation of the IP 3 -and DAG-dependent conductances and Ca 2+ exchangers develops on a wider range of concentrations, i.e. more smoothly.

Perspectives
The field of insect olfactory transduction has generated an ever increasing amount of data but dispersed and fragmentary. The complication and sometimes confusion that result from these circumstances justify an attempt to unite the parts in a comprehensive view and formal synthesis of how things might work. However, because it is clearly incomplete and involves several assumptions, the model presented here should not be considered uncritically as a faithful description of reality. Like any model, it is intended primarily as a starting point for addressing new questions, designing new experiments, and offering a tentative framework for their interpretation.
This model calls for three kinds of experiments. First, experiments on specific components. For example, are the contributions of the ionotropic and metabotropic pathways to olfactory transduction in agreement with our tentative proposal? Is the membrane repolarized by a K + channel depending on an IP 3gated Ca 2+ channel, both located in the inner dendrite as suggested here? What is the exact mechanism of Ca 2+ extrusion? Second, can the model be extended to account for adapted or periodically stimulated neurons? Can it be induced to oscillate like in mammals [41,42,102]. Finally, can qualitative differences in all these mechanisms be found between different ORN types in the same and different biological species? All these experiments have the potential to confirm or invalidate the views presented here and so to yield significant progress in our integrated view of the ORN functioning.

Model Equations
Differential equations for pre-effector events. These equations describe the uptake, perireception, reception and early amplification in the moth pheromone receptor neuron. The species and reactions are defined in Figure 2. The values of rate constants are given in [36]. The initial values of the species are given in Table 2.
Leak current I ld at outer dendrite K + current I K at inner dendrite and soma Leak current I ls at inner dendrite and soma Longitudinal currents from outer dendrite to soma, in auxiliary cell and in sensillar lymph Differential equations for potentials. Potentials inside (V id ) and outside (V ed ) the outer dendrite, inside the inner dendrite and soma (V is ), and outside the auxiliary cell (V ea ), see Figure 6.
Numerical integration. The system of differential equations given above was integrated with the Matlab ode45 solver (The MathWorks, Natick, USA).

Parameter Estimation
The unknown parameters of the model were estimated by utilizing various search methods based on the following criteria. First, we imposed that each parameter be in a physiologically acceptable range of values compatible with the properties of our qualitative model of transduction: the order of magnitude of Hill coefficients is one; the IP 3 -gated channel is permeable to Ca 2+ only while the DAG-gated channel is permeable to Ca 2+ and other cations, so that f Ca .f cat ; the reversal potential of the Cl 2 channel E Cl must be more positive than 297 mV to be depolarizing. Second, we considered a parameter set as acceptable if the predicted kinetics of the sensillar potential were close to the experimentally measured kinetics at all uptakes. For checking this condition, we minimized a cost function based on the three response characteristics, height (H i ), rising time (t rise,i ) and falling time (t fall,i ) at a series of uptakes i for which these characteristics were determined experimentally. The differences, DH i = |H i 2Ĥ i |, between the values Ĥ i predicted by the model for a given set h of parameter values and the experimental values H i , were determined at every uptake i. The differences Dt rise,i and Dt fall,i were determined in the same way. Because the three characteristics vary on different scales, the differences were weighted and summed to produce a single cost function where n = 26 is the number of uptakes. Third, a solution was finally accepted only if it was in qualitative accordance with other available experimental facts: the transient feature of the kinetics of the second messengers and IP 3 -gated currents, the sustained property of the Cl 2 currents and K + currents, and the condition that intracellular Ca 2+ concentration must not exceed 200 mM. Two search methods were utilized in sequence to find the parameter values. First, for a global exploration of the parameter space, we relied on a trial-and-error method. We compared a few thousands parameter sets, drawn from sets h obeying the first criterion above, at 3 uptakes (low, medium, high). Most sets led to unacceptable cost functions E(h) and were rejected. The best sets were further selected on the third criterion then tested at more uptakes. Eight presumptive solutions tested at all uptakes with E(h) in the range 2.74-5.38 were found fulfilling the three criteria. Second, the best presumptive solution for which E(h) = 2.74 was locally optimized utilizing the Matlab unconstrained minimizer fminsearch based on the Nelder-Mead simplex (direct search) method. The algorithm converged on the set of estimated parameters h 0 shown in Tables 4 and 5. With these parameter values E(h 0 ) = 2.61.

Sensitivity Analysis
The sensitivity of a model response M to a single parameter b i can be expressed as a sensitivity function S b i ð Þ~LM=Lb i . This partial derivative was estimated as the central finite difference (59) using both the forward and backward differences This equation is only valid for an infinitesimal variation (perturbation) of b i (Db i ?0). Practically, Db i was implemented as the product Db i = jb i , where b i is the nominal parameter value as estimated above and j is a perturbation factor. We took j = 0.01, large enough to avoid numerical inaccuracies and small enough to prevent the nonlinearity of the model to play a role in the sensitivity calculations. In our case, the model responses M u are the height, half-rise time and half-fall time at different uptakes log 10 U (denoted here as subscript u), which have different units and take values of different orders of magnitudes. In order to compare their sensitivities S u (b i ) we normalized them by the model response M u For each fitted parameter, the normalized sensitivities S ru of the 3 characteristics were calculated at 26 values of U from threshold to saturation. The characteristic and the uptake giving the largest absolute value S r (b i ) of the S ru (b i ) were recorded (see Tables 4  and 5). It is conceivable that in the optimal parameter set h 0 , the low sensitivity of specific parameters is a result of the local optimization procedure. We checked that this was not the case for each lowsensitivity parameter k by calculating the cost function E(h9) where h9 = h 0 , except for k whose value was taken 10% smaller (and 10% larger) than its optimal value. We verified that in all cases E(h9)<E(h 0 ).

Supporting Information
Text S1 Functional significance of nonlinear mechanisms and sensitivity analysis of model parameters. Found at: doi:10.1371/journal.pcbi.1000321.s001 (0.08 MB DOC)