Spatial Distribution of Calcium-Gated Chloride Channels in Olfactory Cilia

Background In vertebrate olfactory receptor neurons, sensory cilia transduce odor stimuli into changes in neuronal membrane potential. The voltage changes are primarily caused by the sequential openings of two types of channel: a cyclic-nucleotide-gated (CNG) cationic channel and a calcium-gated chloride channel. In frog, the cilia are 25 to 200 µm in length, so the spatial distributions of the channels may be an important determinant of odor sensitivity. Principal Findings To determine the spatial distribution of the chloride channels, we recorded from single cilia as calcium was allowed to diffuse down the length of the cilium and activate the channels. A computational model of this experiment allowed an estimate of the spatial distribution of the chloride channels. On average, the channels were concentrated in a narrow band centered at a distance of 29% of the ciliary length, measured from the base of the cilium. This matches the location of the CNG channels determined previously. This non-uniform distribution of transduction proteins is consistent with similar findings in other cilia. Conclusions On average, the two types of olfactory transduction channel are concentrated in the same region of the cilium. This may contribute to the efficient detection of weak stimuli.


Introduction
We model the behavior of cytoplasmic Ca 2+ , buffer and membrane potential in an olfactory cilium during experiments that involve the Ca 2+ -gated Cl − ion channels. The primary goal of this work is to elicit information on the distribution of the Ca 2+ -gated Cl − ion channels.
Our main goal in this supplement is to restate the mathematical model which was introduced in Badamdorj [46] and Badamdorj et al. [20] as well as the rapid buffer simplification.

Mathematical Model
Our model consists of equations for the membrane potential v = v(x, t) and calcium concentration c = c(x, t) where 0 < x < L and L is the length of the cilium. The point where x = 0 is the open (proximal) end of the cilium, and x = L corresponds to the closed (distal) end. We assume that time t is in the range of several seconds.
The membrane potential satisfies 1 r a where r a is ciliary intracellular axial resistance and J is the transmembrane current flow through the Ca 2+ -gated Cl − channels. We have assumed the capacitance term is negligible as well as the background conductance (leak current). The Cl − current is given by where g Cl is the single Ca 2+ -gated Cl − channel conductance and ρ = ρ(x) is the distribution of the Ca 2+ -gated Cl − channels along the length of the cilium. The Hill function represents Ca 2+ molecules binding and activating the Ca 2+ -gated Cl − channels [47]. Here n is the Hill constant and K 1/2 is the half-maximal concentration. To complete the description of the membrane potential problem we append the boundary conditions below v(0, t) = v Bulk and ∂v ∂x (L, t) = 0.
Here v Bulk is the voltage at which the membrane potential is clamped at the open end. The behavior of the cytoplasmic calcium and buffer complex b(x, t) = [CaB] can be modeled by the following initial/boundary value problems [48] ∂c ∂t

S1
and where k + and k − are the rate constants for uptake and release of Ca 2+ from the buffer. Here D Ca and D B are the diffusivities for the Ca 2+ and buffer, respectively. We neglect any Ca 2+ pumps, since no ATP is supplied. B T is the total concentration of buffer (complex plus free buffer). We neglect calmodulin, which has multiple Ca 2+ binding sites associated with the CNG channels [49]. The model is not very sensitive to the number of Ca 2+ binding sites (see the Results section of the main paper).
Assuming that the rapid buffer approximation (RBA) is appropriate (buffer reaction is faster than the other processes involving calcium) we set K B = k − /k + then combine (4) and (5) to obtain [48,50] This model for the evolution of calcium is complete after adding the boundary conditions and initial condition c(x, 0) = 0.
We have represented the fact that there is a free calcium concentration c Bulk outside of the pipette (at the open end of the cilium) and initially the concentration inside is so small that we take it to be zero. We are especially interested in the behavior of the current, which is experimentally measurable. Our belief that the RBA is appropriate arises from a nondimensionalization of (4) using K 1/2 for c and B T for b. This reveals that the buffer reaction term is much larger than the other terms that model the diffusion and binding of Ca 2+ . (In fact, the factor is around 10 6 ). So the reaction term can be solved independently to determine b in terms of c and then substituted into the new equation obtained by adding (4) to (5) leading to (6).
Continuing to follow Keener & Sneyd [48], we can reformulate our model by defining where φ −1 is a one-to-one function which we can invert by solving an intermediate quadratic equation Again we append boundary conditions and initial condition w(x, 0) = 0.
Remark: In situations where 2c/K B << 1, it can be shown that equation (6) is given approximately by and . If we further assume that the impact of the binding is small and the cilium is long so we can neglect the boundary condition at the tip (resulting in a half-line problem), then a closed form solution for c can be determined. This reduction is reasonable in the case when the buffer is HEDTA since, then, K B = 14.8 µM.

Inverse Solver
Previous studies [9,19,20] suggest that tall, narrow Gaussian channel distributions provide current evolution profiles that fit the experimental current data well (see also Fig. 2 of Results). We thus search for x 0 , ρ 0 , and δ w so that the current I(t) resulting from in our model is a good match to the measured current data for a given specific cilium. Here ρ 0 is the peak height, x 0 is the location of that peak and δ w is a width. Note the following approximate equality that holds if the Gaussian is truly narrow with respect to the length of the cilium L and not centered close to the ends where T Cl is the total number of channels. So the matching involves choosing x 0 , ρ 0 , and δ w so that is minimized where (t 1 , I Data 1 ), . . . , (t NData , I Data NData ) are the time and current data points arising from the experiments (e.g. those plotted in Fig. S4 below). Our search procedure uses asymptotic formulas from Badamdorj et al. [20] in concert with a dichotomous search algorithm to provide a good guess for x 0 followed by the MATLAB fminsearch code that uses the Nelder-Mead algorithm to find the three points x 0 , ρ 0 , and δ w . We describe the dichotomous search part of the overall algorithm more precisely. The reduced problem studied in Badamdorj et al. [20] indicated that if t 1/2 is the time the current profile has reached half its peak height I Cl (see Fig. S4) and D Avg = 1 2 (D Ca + D B ) then So we search for a minimizing argument x 0 of where we used (12) and the second formula in (13) for ρ 0 in the sum of squares error formula S.

Dichotomous Search Algorithm:
Choose [a D , b D ], ǫ D , and i Max . Set i = 0 and δ w = 1.
The dichotomous search code assumes the function F M is unimodal. It is a bracketing approach like the bisection method for rootfinding. Successive subintervals are determined by bisection and calculation on the half subinterval in which the minimum resides. A distinguishability coefficient ǫ D is used for this (0 < ǫ D << 1). A bracketing interval [a D , b D ] ⊂ [0, L] is initially found from the guess for x 0 in (13). We have generally chosen δ w = 1 in this search. Of course in the MATLAB search we allow all three parameters (x 0 , ρ 0 , δ w ) to vary.
Usually, only a moderate number of iterations are done, since there is the second Nelder-Mead search as well. This fminsearch or Nelder-Mead search starts with a guess using the x 0 from the dichotomous search, a value for ρ 0 from (12) and (13), and δ w = 1.

Example
Here we display the graphs from a sample run of our code which searches for an optimum (x 0 , ρ 0 , δ w ) set of parameters. The results in Fig. S5 are for a 25-µm cilium. The inverse solver was used and produced an answer with a relative least squares error where E 2 = 0.024 with Fig. S5A is a plot of the channel distribution function ρ with x 0 = 14.4 µm, δ w = 0.917, and ρ 0 = 1489 channels. Here the channel density was 96.8 channels/µm and total number of channels T Cl = 2420. In Fig. S5B are shown the current data resulting from a forward solution using ρ (gray) as well as the experimental current data (black). In Fig. S5C there is a plot of F M (x) vs x revealing that it is unimodal in this example. Notes: • Where uncertainties are shown they represent SEM. Uncertainties given for n and K 1/2 are unpublished values from the data in Kleene & Gesteland [21].
• For the conversion factor above, α = 1/ (N A A) where A is the cross-sectional area of the cilium and N A is Avogadro's number.
• The reaction rates are apparent dissociation constants between Ca 2+ and the buffer, i.e. the inverses of the apparent association constants determined as described in the Methods section of the main paper.