Modelling and sensitivity analysis of the reactions involving receptor, G-protein and effector in vertebrate olfactory receptor neurons

Abstract

A biochemical model of the receptor, G-protein and effector (RGE) interactions during transduction in the cilia of vertebrate olfactory receptor neurons (ORNs) was developed and calibrated to experimental recordings of cAMP levels and the receptor current (RC). The model describes the steps from odorant binding to activation of the effector enzyme which catalyzes the conversion of ATP to cAMP, and shows how odorant stimulation is amplified and delayed by the RGE transduction cascade. A time-dependent sensitivity analysis was performed on the model. The model output—the cAMP production rate—is particularly sensitive to a few, dominant parameters. During odorant stimulation it depends mainly on the initial density of G-proteins and the catalytic constant for cAMP production.

Appendices

Appendix A: Equations of the RGE model

$$ \frac{dR}{dt} = - k_1 \cdot R \cdot L + k_{- 1} \cdot \left( {RL} \right) $$

(A1)

$$ \frac{{d\left( {RL} \right)}}{dt} = k_1 \cdot R \cdot L - \left( {k_2 + k_{- 1} } \right) \cdot \left( {RL} \right) + k_{- 2} \cdot R_a $$

(A2)

$$ \frac{{dR_a }}{dt} = k_2 \cdot \left( {RL} \right) - k_{- 2} \cdot R_a $$

(A3)

$$ \frac{dG}{dt} = - e_{RG} \cdot R_a \cdot G + e_G \cdot G_r \cdot G_b $$

(A4)

$$ \frac{{dG_a }}{dt} = e_{RG} \cdot R_a \cdot G - e_{GE} \cdot G_a \cdot E - k_{Ga} \cdot G_a $$

(A5)

$$ \frac{{dG_b }}{dt} = e_{RG} \cdot R_a \cdot G - e_G \cdot G_r \cdot G_b $$

(A6)

$$ \frac{{dG_r }}{dt} = k_{Ga} \cdot G_a + k_{Ea} \cdot E_a - e_G \cdot G_r \cdot G_b $$

(A7)

$$ \frac{dE}{dt} = - e_{GE} \cdot G_a \cdot E + k_{Ea} \cdot E_a $$

(A8)

$$ \frac{{dE_a }}{dt} = e_{GE} \cdot G_a \cdot E - k_{Ea} \cdot E_a $$

(A9)

$$ y = k_c \cdot E_a $$

(A10)

In addition the following conservation laws apply:

$$ R_{tot} = R + RL + R_a $$

(A11)

$$ G_{tot} = G + G_a + G_r + E_a = G + G_b $$

(A12)

$$ E_{tot} = E + E_a $$

(A13)

The odorants (in the mucus or physiological solution) interact with the receptors (on the cilia). The mucus is a volume (3D) while the receptors reside on the ciliary surface (2D). In Eqs. (A1–A13) all species are membrane-bound, so that their amounts are given as membrane densities (molecules μm⁻²), except L and y which are in solution and given in μM. In Eq. (A1), the dimension of k ₁ ·L·R must be molecules μm² s⁻¹, indicating that k ₁ must have units s⁻¹ μM⁻¹. Similarly, if the cAMP production rate y (Eq. (A10) should have dimension (s⁻¹ μM), k _c must have units s⁻¹ μM molecules⁻¹ μm². The dimensions of the other rate constants (Table 2) are more straightforward to verify. The electronic supplementary material contains a Matlab-script for solving the dynamic equations of the RGE model (RGE.m).

Appendix B: The Sobol and Saltelli method for sensitivity analysis

The SA method developed by Sobol and Saltelli (Saltelli 2002; Saltelli 2004), summarized here, treats the model output y = f(x ₁ , x ₂ … x _k ) as a function of k input factors x _i which are the initial conditions and rate constants of the model. The output y may be a scalar or a time series. A time-dependent output (as in our case) can be analyzed by this method because sensitivity indices are computed at each time step separately.

The computation of the output sensitivity to the input factors involves two sample matrices

$$ M_1 = \left[ {\begin{array}{*{20}c} {x_1^{{(1)}} } \hfill & {x_2^{{(1)}} } \hfill & \ldots \hfill & {x_k^{{(1)}} } \hfill \\ {x_1^{{(2)}} } \hfill & {x_2^{{(2)}} } \hfill & \ldots \hfill & {x_k^{{(2)}} } \hfill \\ \ldots \hfill & \ldots \hfill & \ldots \hfill & \ldots \hfill \\ {x_1^{{(N)}} } \hfill & {x_2^{{(N)}} } \hfill & \ldots \hfill & {x_k^{{(N)}} } \hfill \end{array} } \right] $$

(B1)

and

$$ M_2 = \left[ {\begin{array}{*{20}c} {x_1^{{\left( {1\prime } \right)}} } & {x_2^{{\left( {1\prime } \right)}} } & {...} & {x_k^{{\left( {1\prime } \right)}} } \\ {x_1^{{\left( {2\prime } \right)}} } & {x_2^{{\left( {2\prime } \right)}} } & {...} & {x_k^{{\left( {2\prime } \right)}} } \\ {...} & {...} & {...} & {...} \\ {x_1^{{\left( {N\prime } \right)}} } & {x_2^{{\left( {N\prime } \right)}} } & {...} & {x_k^{{\left( {N\prime } \right)}} } \end{array} } \right], $$

(B2)

each containing a set of N Monte Carlo samples for the k input factors. From M ₁ and M ₂ , a third matrix N _j can be constructed by replacing the Monte Carlo set for input parameter j in M ₂ by the set from M ₁ .

$$ N_j = \left[ {\begin{array}{*{20}c} {x_1^{{\left( {1\prime } \right)}} } & {x_2^{{\left( {1\prime } \right)}} } & {...} & {x_j^{{(1)}} } & {...} & {x_k^{{\left( {1\prime } \right)}} } \\ {x_1^{{\left( {2\prime } \right)}} } & {x_2^{{\left( {2\prime } \right)}} } & {...} & {x_j^{{\left( {2\prime } \right)}} } & {...} & {x_k^{{\left( {2\prime } \right)}} } \\ {...} & {...} & {...} & {...} & {...} & {...} \\ {x_1^{{\left( {N\prime } \right)}} } & {x_2^{{\left( {N\prime } \right)}} } & {...} & {x_j^{{\left( {N } \right)}} } & {...} & {x_k^{{\left( {N\prime } \right)}} } \end{array} } \right] $$

(B3)

The first order and total order sensitivity indices can be computed according to:

$$ \hat{S}_j = \frac{{\hat{V}\left( {\hat{E}\left( {Y\left| {X_j } \right.} \right)} \right)}}{{\hat{V}(Y)}} = \frac{{\hat{U}_j - \hat{E}^2 (Y)}}{{\hat{V}(Y)}}, $$

(B4)

and

$$ \hat{S}_{Tj} = 1 - \frac{{\hat{U}_{- j} - \hat{E}^2 (Y)}}{{\hat{V}(Y)}}. $$

(B5)

U _j is obtained from the product of f-values computed from the sample matrix M ₁ and f-values computed from N _j where all factors except X _j have been re-sampled. Let f(M ^r ₁ ), f(M ^r ₂ ) and f(N ^r _j ) denote the output y obtained by using the r’th parameter set in M ₁ , M ₂ and N _j respectively. We then have:

$$ \begin{array}{*{20}c} {\hat{U}_j = \frac{1}{N - 1}\sum\limits_{r = 1}^N {f\left( {x_1^{{(r)}}, x_2^{{(r)}}, ...,x_k^{{(r)}} } \right)f\left( {x_1^{{\left( {r\prime } \right)}}, x_2^{{\left( {r\prime } \right)}}, ...,x_j^{{(r)}}, ...x_k^{{\left( {r\prime } \right)}} } \right)} } \\ { = \frac{1}{N - 1}\sum\limits_{r = 1}^N {f\left( {M_1^r } \right)f\left( {N_j^r } \right)} } \end{array} $$

(B6)

The rule is: Resample all parameters except the one (j) that should be estimated. The total effect of j is obtained by subtracting the effect of all other parameters from unity. Thus, U _−j measures the sensitivity to all parameters except j by re-sampling only j:

$$ \hat{U}_{- j} = \frac{1}{N - 1}\sum\limits_{r = 1}^N {f\left( {M_2^r } \right)f\left( {N_j^r } \right)} . $$

(B7)

The unconditional mean and variance (E(Y) and V(Y)) can be estimated in several ways, using either of the sets M ₁ or M ₂ . For an infinite number of Monte Carlo simulations, all definitions should be equivalent. However, since our simulations are time consuming, we limited the number of Monte Carlo simulations to 5,000 iterations, which is not very high for exploring a 13 dimensional parameter space. Due to the limited number of iterations, the SA-method is sensitive to how the sample sets (eqns. B1–B3) are combined in the expressions for the mean and variance. We found that some minor and statistically valid adjustments of the original method, lowers the probability of getting theoretically implausible results (e.g. negative sensitivity indices, or the sum of first order sensitivity indices exceeding unity). To make optimal choices, we consider two special cases: 1) when j is a non-influential variable, so that the result should not depend on whether j is re-sampled or not, and 2) when j is the only influential variable, so that the result should not be influenced by whether all other variables are re-sampled or not. These two cases are now studied for the computation of the first order sensitivity indices and the total order sensitivity indices respectively.

The first order indices

If j is a non-influential variable, f(N ^r _j ) should be identical to f(M ^r ₂ ) in Eq. (B6), and S _j should be zero. If j is the only important parameter, f(N ^r _j ) should be identical to f(M ^r ₁ ), and S _j should be unity. These two conditions can be written:

$$ \frac{1}{N - 1}\sum\limits_{r = 1}^N {f\left( {M_1^r } \right)f\left( {M_2^r } \right)} - \hat{E}^2 (Y) = 0, $$

(B8)

and

$$ \frac{1}{N - 1}\sum\limits_{r = 1}^N {f\left( {M_1^r } \right)f\left( {M_1^r } \right)} - \hat{E}^2 (Y) = \hat{V}(Y). $$

(B9)

The former case suggests that we use the valid estimate

$$ \hat{E}^2 (Y) = \frac{1}{N}\sum\limits_{r = 1}^N {f\left( {M_1^r } \right)f\left( {M_2^r } \right)} $$

(B10)

for the squared unconditional mean (which satisfies the condition to a factor N/(N  −  1)). The second case suggest the usual definition of variance, so that

$$ \hat{V}(Y) = \frac{1}{N - 1}\sum\limits_{r = 1}^N {f^2 \left( {M_1^r } \right)} - \hat{E}^2 (Y), $$

(B11)

with Ê ² (Y) defined as in (B10). Note that a truly non-influential variable in this case would be overestimated due to factor N instead of N  −  1 in the definition of the mean.

The total effect indices

If j is a non-influential variable, f(N ^r _j ) should be identical to f(M ^r ₂ ) in Eq. (B7), and S _−j should be zero. If j is the only important variable, f(N ^r _j ) should be identical to f(M ^r ₁ ), and S _−j should be unity. These two conditions:

$$ \frac{1}{N - 1}\sum\limits_{r = 1}^N {f\left( {M_2^r } \right)f\left( {M_2^r } \right)} - \hat{E}^2 (Y) = \hat{V}(Y), $$

(B12)

and

$$ \frac{1}{N - 1}\sum\limits_{r = 1}^N {f\left( {M_1^r } \right)f\left( {M_2^r } \right)} - \hat{E}^2 (Y) = 0 $$

(B13)

give the corresponding definitions:

$$ \hat{E}^2 (Y) = \frac{1}{N}\sum\limits_{r = 1}^N {f\left( {M_1^r } \right)f\left( {M_2^r } \right)}, $$

(B14)

and

$$ \hat{V}(Y) = \frac{1}{N - 1}\sum\limits_{r = 1}^N {f^2 \left( {M_2^r } \right) - \hat{E}^2 (Y)} $$

(B15)

for the mean and variance. The only difference is that the first summation in the variance in this case is based on M ₂ , as opposed to M ₁ in the first order indices. Note that a truly non-influential variable in this case would result in the exact estimate S _−j  = 0.

In the original work (Saltelli 2002, 2004), the more common definition

$$ \hat{E}^2 (Y) = \left[ {\frac{1}{N}\sum\limits_{r = 1}^N {f\left( {M_2^r } \right)} } \right]^2 $$

(B16)

was used in the computation of the total effect indices. These various definitions of the expected values and variance are statistically plausible and should give identical result in the case of an infinite number of iterations. Due to the limited data set and independent computations of the indices, the stochastic possibility of obtaining negative sensitivity indices was still realized at some time points, but then only for very small sensitivity indices (see Table 4). Still, the reported conclusions on important vs. unimportant input parameters were found repeatedly in several independent simulations.

Another way of making the simulations more robust, proposed in Saltelli (2004), is to subtract the mean from Y, and replacing Y with the variable Y  −  Ê(Y). That is, replace the variable f(M ₁ ) with

$$ \tilde{f}\left( {M_1 } \right) = f\left( {M_1 } \right) - \frac{1}{N}\sum\limits_{r = 1}^N {f\left( {M_1^r } \right)}, $$

(B17)

in Eqs. (B6–B15) above, and do the corresponding operations for f(M ₂ ) and f(N _j ). This was used in all our simulations.

Appendix C: Sensitivity analysis and parameter estimation in the RGE model

The SA considers parameter variations within an interval around the estimated parameter set. In our study, these variations were sampled from a normal distribution, using the estimated parameters in Table 2 as the mean and 15% of mean as standard deviation (std). We wanted a sampling interval that could reflect some realistic parameter variations (e.g. due to measurement errors), and not only small perturbations. At the same time, we wanted the sampling interval to be narrow enough to reflect the system sensitivity for the specific estimated parameter set. For instance, the RGE(D05) model, where G ₀  < E ₀ , was always most sensitive to G ₀ , while the RGE(TK05) model, where G ₀  > E ₀ , was most sensitive to E ₀ when the odour concentration was very high. If the sampling interval was too broad, such differences in the scenario would not be reflected in the SA, so that an interval std = 0.15·mean. is a good compromise. This will also make the probability of sampling negative parameter values close to zero.

The sensitivity of the model output to a given parameter i has two related implications. First, if the model is used to predict an output for which no empirical data exist, an accurate empirical measurement of parameter i is crucial for a good prediction. Second, if the model is calibrated to existing data (as in our case), it means that the model can estimate parameter i with high accuracy, since any small deviance from its correct value would give a strong discrepancy between the model output and the empirical data. Our model should thus be able to estimate G ₀ and k _c with high accuracy. As can be seen in Fig. 8, the parameters k ₋₁ and k _Ea are mainly influential after the stimulus has been turned off. Since the TK05 experiments were conducted with continuous stimulus, we used the D05 data to estimate the generic parameter k _Ea , and then used the same value in both data sets. For k ₋₁ , we did not demand similar values for the two data sets, and k ₋₁ was estimated for the TK05 set during stimulation. As can be seen in Fig. 8, k ₋₁ is mainly influential for the post-stimulus cAMP, for which there is no data, indicating that the post-stimulus cAMP signal for the TK05 set should be regarded with some caution.

Conversely, the least important parameters in our RGE model (i.e. k _Ga e _G , e _GE and E ₀ ) can not be estimated with any accuracy. This means that the model is robust to perturbations of these parameters. Furthermore, comparison of stimulations with low and high odorant concentrations (Fig. 10) indicates that the relative importance of the various parameters on the model response depends on the odorant concentration. Parameter k _c is best estimated at high odorant stimulation both for the TK05 and D05 ORN (Fig. 10(b, d)). G ₀ is best estimated at high odorant concentration in the D05 ORN (Fig. 10(b)), and at intermediate concentration for the TK05 ORN (Fig. 9(b)), while E ₀ is only influential at high odorant stimulus for the TK05 ORN (Fig. 10(d)). The remaining parameters are best estimated at low concentration both for the D05 and TK05 ORN, because in this case the output is more evenly sensitive to several of the input parameters (Fig. 10(a, c)). The sensitivity to R ₀ is highest at low concentrations, indicating that a high density of receptor molecules may only be important for detecting weak signals. However, the sensitivity to R ₀ was quite low also in this case, and at intermediate stimulus, the sensitivity to R ₀ was lower than for the other proteins even when it was the sparsest species (Section 3.3). Our estimate of R ₀ (in the range 700–1,700 molecules µm⁻²) is relatively low in comparison with R ₀  = 6,000 molecules.µm⁻² estimated in moths (Kaissling 2001).

As pointed out earlier, a complex model may have several parameter combinations that may give a good fit to a single output channel. Since the parameter estimation procedure does not guarantee a global minimum, we cannot exclude the possibility of another set of parameters that might give a good fit to the cAMP data. As shown in Section 3.1, the estimated parameters represent a local fit that depends on the initial conditions of the estimation algorithm. An accurate estimate of a parameter is thus relative to this local optimum. This point is most critical concerning the protein densities, for which we are not aware of any empirically estimated values in vertebrate ORNs. If future research should indicate that the real values for G ₀ and E ₀ diverge strongly from ours, the model parameters should be updated, and another local optimum should be searched for.