Multiple Steps of Phosphorylation of Activated Rhodopsin Can Account for the Reproducibility of Vertebrate Rod Single-photon Responses

Reliable single-photon detection is a ubiquitous feature of vertebrate visual systems (Rodieck, 1998). Evolution has produced biochemical and biophysical machinery in photoreceptors that can support reliable single-photon detection despite the inherent high variability in all molecular reactions. An abiding problem in phototransduction has been to explain the observed reproducibility of vertebrate rod single-photon responses (SPRs), given that the amplitude and time course of the response are determined primarily by the lifetime and activity of a single activated rhodopsin molecule (R*). If R* were to inactivate in a single memoryless step with first-order kinetics, its lifetime would be highly variable, exhibiting an approximately exponential distribution with a coefficient of variation (CV = SD/mean) of unity. In the absence of other mechanisms, this highly variable R* lifetime would be reflected, to one degree or another, in the variability of the SPR amplitude and/or kinetics. The relative reproducibility of SPRs, despite the inherent variability of the underlying individual biochemical reactions, places strong constraints on any model of phototransduction.

Empirical measures have shown that the variability of vertebrate rod SPRs (assessed by any of several measures of amplitude or kinetics) is considerably lower than expected for a one-step R*-inactivation process. In this paper, we focus on four findings reported in the literature. The first finding is that the coefficient of variation of the amplitude of the SPRs (CV_ampl) is quite small (0.2–0.25: Baylor et al., 1979; Schneeweis and Schnapf, 1995; Rieke and Baylor, 1998a; Whitlock and Lamb, 1999; Field and Rieke, 2002). The second finding is that the ensemble variance of a set of dim-flash responses is very nearly proportional to square of the mean of the response over almost the entire response time course. This observation has been used to suggest that variability in the elementary response waveform must be small (Schnapf, 1983; Schneeweis and Schnapf, 1995; Rieke and Baylor, 1998a). That argument has been rejected by Whitlock and Lamb (1999), and we will elaborate on their insight that this finding alone does not imply that the SPR has low variability over its entire time course. The third finding is the low coefficient of variation of the area (CV_area ≈ 0.3) under the SPR (Field and Rieke, 2002). This measure has a natural pragmatic appeal, in that it captures variability in both the amplitude and time course of the SPR (Field and Rieke, 2002), and we will show furthermore that there is a theoretical basis for choosing variability in SPR area to estimate the number of inactivation steps of R* over measures of variability of either amplitude or kinetics alone. The fourth measure is based on the finding that the variance of the SPR is at least an order of magnitude smaller than the square of the mean response, until some time after the mean response reaches its peak (Rieke and Baylor, 1998a; Field and Rieke, 2002). Field and Rieke (2002) found that the relative times-to-peak of the SPR variance and SPR squared mean, as well as the width of the SPR variance waveform aided in discriminating between candidate models.

Based, in part, on the ideas and biochemical data of Gibson et al. (2000), we develop in this paper a kinetic model of phototransduction that includes a detailed stochastic simulation of the activation and inactivation of rhodopsin, G-protein (G), and phosphodiesterase (PDE). Simulations of the stochastic activation of PDE have been done previously by Felber et al. (1996) and Lamb (1994). Our modeling is distinguished from theirs in four main respects: (a) We model reaction kinetics in greater detail than the prior studies, including competitive, stochastic binding of R* with G-protein, arrestin (Arr), and rhodopsin kinase. (b) Rather than simulating the two-dimensional diffusional contact between molecules, we concentrate on the reactions of those molecules by assuming that the disc surface is well mixed. (c) Neither Lamb (1994) nor Felber et al. (1996) coupled the activated PDE to the downstream reactions generating photocurrent (including cGMP-hydrolysis, channel closure, Ca²⁺-feedback, and Ca²⁺-buffering). Our model does so, thus allowing direct comparisons between empirical electrophysiological responses and model responses under a wide range of conditions. (d) Lamb (1994) did not simulate inactivation reactions.

Our model embodies the following four hypotheses: (a) that R* undergoes a series of sequential phosphorylation steps, mediated by rhodopsin kinase (Kühn and Wilden, 1982; Wilden and Kühn, 1982; Aton et al., 1984; Thompson and Findlay, 1984; Palczewski et al., 1991; Wilden, 1995); (b) that competition for binding to R* occurs between three molecular species: inactive G-protein (G·GDP), rhodopsin kinase (RK), and Arr (Pfister et al., 1983; Miller and Dratz, 1984; Buczylko et al., 1991; Pulvermuller et al., 1993; Krupnick et al., 1997; Gibson et al., 2000); (c) that each sequential step of phosphorylation leads to a progressive decrease in affinity between R* and G·GDP, and concomitantly to a progressive increase in affinity between R* and Arr (Gibson et al., 2000); and (d) that the affinity between R* and RK also ratchets down with each step of phosphorylation (Buczylko et al., 1991; Gibson et al., 2000). We will refer to our full model based on these premises as the sequential phosphorylation model.

Using this model, we simulate ensembles of SPRs and dim-flash responses using Monte-Carlo methods. The model accounts for the four quantitative measures of SPR variability/reproducibility delineated above. In addition, despite the fact that in the model phosphorylation is the dominant process inactivating R*, the simulations reproduce key experimental results that Rieke and Baylor (1998a) interpreted previously as evidence that phosphorylation could not be the dominant deactivation mechanism.

Moreover, without any additional parameter adjustments, the model reproduces the salient qualitative features of SPRs obtained from four genetic knockout and transgenic studies: three studies in which mouse rods had been genetically modified with specific pathways of R* inactivation disabled (Arr−/−, Xu et al., 1997; RK−/−, Chen et al., 1999; CSM, Mendez et al., 2000), and one study in which the mechanism of feedback synthesis of guanylate cyclase had been disabled (GCAPs−/−; Burns et al., 2002). In general, we find that attempting to simulate the results of knockout and transgenic experiments can provide valuable insights into the kinds of mechanisms that must be present, or about candidate molecular schemes that can be ruled out even if, in principle, they could achieve the empirical SPR reproducibility in one or more of the four quantitative measures of reproducibility.

Finally, we present a theoretical analysis of mechanisms underlying SPR reproducibility. It reveals a heretofore unrecognized tradeoff between variability of SPR amplitude and variability of SPR kinetics that depends critically on R* lifetime in relation to the net kinetics of the downstream reactions in the cascade. We explain what the statistics of SPR area tell us about underlying molecular mechanisms. We show how the coefficient of variation of SPR area can be used, in the context of a certain class of models, to estimate a lower bound for the number of R* inactivation steps required to yield observed SPR reproducibility. Neither the CV of SPR amplitude nor the statistics of SPR duration variability alone can be used in this manner. An important conclusion of this theoretical analysis is that, in normal rods, across species, the kinetics of R* inactivation cannot be very different from the kinetics of the downstream reactions, including the inactivation of the activated transducin–phosphodiesterase complex, in the dim-flash regime.

To simulate the rod's response to dim flashes, we formulated an explicit model of the rod phototransduction cascade, including detailed stochastic implementation of the initial reactions. The time-varying quantities of reactants in these “front-end” reactions, defined as those from photon absorption and isomerization of R* through the activation and inactivation of PDE* (Eqs.

The stochastic front-end model is based on the representation of the underlying biochemistry illustrated in Fig. 1. We now present the assumptions on which we based this model, and we specify the individual reaction steps in terms of a set of equations. In these equations, n refers to the phosphorylation state of R*, and the values of the parameters are given in Table I.

We assume (Eqs.

with the final, complete quench of R* activity occurring upon Arr binding to phosphorylated R*:

The notations “pre” and “post” in Eqs.

Transducin is assumed to be activated by a conventional series of reactions (Eqs.

where G_α·GTP is the activated form of transducin.

Based on the biochemical results of Gibson et al. (2000), the affinity of rhodopsin for G-protein is assumed to decrease exponentially with increasing phosphorylations (Eq. 4), while its affinity for Arr is assumed to increase linearly with n (Eq. 5). Arrestin is able to bind R*, with increasing probability, at any time following the first phosphorylation (n ≥ 1), to terminate whatever R* activity remains in that state.

The exponential rate, ω, was set to 0.6 from an exponential fit to the data in Fig. 2 A in Gibson et al. (2000).

Although the explicit dependence of R–RK binding affinity and rate on R* phosphorylation state has not been documented with the same detail as the R–G and R–Arr affinities (Gibson et al., 2000), there is biochemical support for the notion that R–RK affinity decreases systematically with the phosphorylation state of R*. For example, Buczylko et al. (1991) found that phosphorylated RK has significantly lower affinity for phosphorylated R* than for unphosphorylated R*. Moreover, the data and analysis of Gibson et al. (2000) shows that, given the opportunity for RK to interact with phosphorylated versus unphosphorylated rhodopsin, it preferentially interacts with unphosphorylated rhodopsin, consistent with the idea that R–RK affinity decreases systematically as R* becomes phosphorylated.

We set the R–RK affinity to have the same dependence on n as the affinity between R* and G-protein by varying k_RK1.

We treat the rate of incremental phosphorylation of the R–RK complex by ATP as n dependent, in that it reduces to zero when all available sites on the R* carboxy terminus are occupied.

In a model where RK and G bind R* competitively and Arr is only permitted to bind once R* is fully phosphorylated, this arrangement (Eqs.

We also examined the effect of having a different exponential phosphorylation dependence for the R–G and R–RK affinities. We found that (assuming the ω for R–G affinity in Eq. 4 was fixed at the empirical estimate of 0.6), if ω for R–RK affinity in Eq. 6 varied between ∼0.3 and ∼0.7, the predicted variability of PDE* production would only vary by ∼5 percent. Finally, we examined other profiles for R–RK affinity, such as linear decrease in affinity with n, or a flat profile (no dependence on phosphorylation state). These profiles led to a significant increase in the predicted variability in the number of PDE*s produced compared to the exponential profile used in the sequential phosphorylation model.

Because the phosphorylation dependence profile for R–RK affinity we adopted yields good statistical performance of the model (e.g., low SPR variability), and because some of the details of such phosphorylation dependence have yet to be fleshed out experimentally, our choice of phosphorylation dependence for R–RK affinity may be viewed as a testable biochemical prediction of our model.

For the activation of PDE*, we assumed that activated transducin binds to the γ subunit of PDE (Eq.

where PDE*·G_α · GTP represents the activated form of the transducin–PDE complex, referred to as PDE* throughout this paper, and PDE·G_α · GDP represents the inactivated form (PDE). In all simulations, unless otherwise stated, τ_PDE was 3 s. The mean stochastic lifetime of R* (with all effects of multiple phosphorylation and Arr-binding included) was ∼2.8 s. However, the effective first-order time constant (defined as the first moment of the mean R* activity function) approximating the mean R* inactivation waveform was ∼1.3 s. Thus, PDE–inactivation was (on average) the rate-limiting front-end reaction.

.

The front-end steps contribute considerable amplification, and Leskov et al. (2000) report that a single amphibian rod R* molecule activates ∼150 PDE*/s.³ The front-end parameters we used reproduce this PDE* activation rate up to the time of the first phosphorylation of R*. However, as phosphorylation of R* proceeds, the rate of PDE* production decreases (due to the phosphorylation-dependent decrease in R–G affinity, Eq. 4). The net result is that each R* leads to the production of, on average, ∼220 PDE*.

Since the number of activated PDE* molecules per R* is reasonably large, it is appropriate to treat the downstream “back-end” reactions as continuous signals, and to model them using a system of differential equations (Eqs.

where PDE* = [PDE*·G_α · GTP], g = [cGMP], c = [Ca²⁺], and c_b represents intracellular calcium buffers (all in μM), and the “dot” notation represents the time derivative. The parameters employed in these equations are: α_max, the maximal rate of cGMP synthesis; K_c, the concentration of Ca²⁺ at which cGMP is synthesized at half its maximal rate; m, the Hill coefficient for the action of Ca²⁺on cyclase; β_dark, the dark rate of cGMP hydrolysis; β_sub, the rate constant of cGMP hydrolysis per activated PDE subunit; f_Ca, the fraction of circulating current carried by Ca²⁺; F, the Faraday constant; v_cyto, the effective cytoplasmic volume; J_dark, the dark circulating current; g_dark, the steady-state concentration of cGMP in darkness; n_cg, the Hill coefficient for channel opening by cGMP; γ_Ca, the rate constant of calcium extrusion by the Na⁺/Ca²⁺-K⁺ exchanger; c₀, the minimum intracellular Ca²⁺ concentration; k₁ and k₂, the association and dissociation rate constants of Ca²⁺ with intracellular buffers; and e_T, the total calcium buffer concentration. The units and values of all parameters are given in Table I.

Photocurrent is given by

We assumed that Ca²⁺ feedback occurs at guanylate cyclase only. Having cyclase feedback as the dominant feedback mechanism is consistent with recent evidence from GCAPs−/− rod responses (Burns et al., 2002), suggesting that Ca²⁺ feedback via recoverin and RK onto R* phosphorylation rate does not significantly affect the dim-flash response, though it does affect responses at high intensities.

The values of the parameters for the sequential phosphorylation model (Eqs.

In order to capture the stochastic nature of the “front-end” of the cascade, in addition to R*'s interaction with Arr and RK, we implemented seven activation steps from photon absorption to PDE activation (Eqs.

The back-end parameters K_c, β_dark, f_Ca, k₁, k₂, and e_T were optimized so that the model would generate flash responses that matched a representative set of empirical flash responses obtained from toad rods (from Rieke and Baylor, 1998b) over a wide range of flash intensities. More details about the choices of parameter values can be found in the legend to Table I.

All aspects of the model and analyses were implemented using Matlab/Simulink (The Mathworks). The code is available upon request (contact first author).

Monte-Carlo simulations of 1,000 SPRs were run. Simulations were carried out using the Gillespie method (Gillespie, 1976, 1977; Felber et al., 1996). The probability that a molecular species takes a particular reaction pathway is given by the reaction rate for entering that pathway divided by the sum of rates of all available pathways. The dwell time or waiting time before taking one of the pathways is an exponentially distributed random variable with a mean given by the inverse sum of reaction rates for the possible pathways. For example, a free R* molecule's next event could be a binding with either G·GDP, RK, or Arr. To determine which of these happens, we generate a random number, r, uniformly distributed between 0 and 1. G·GDP binds if r < k_G1

To simulate experimental data we generated a series of PDE* responses to Poisson numbers of photoisomerizations. Failures are all-zero vectors, and multiples are generated by summing the PDE* responses from randomly selected single-photon responses. The PDE* response is transformed to photocurrent through the deterministic back-end of the model (Eqs.

In order to simulate alternate conditions (including effects of genetic modification of R* inactivation mechanisms), we made the following modifications to the sequential phosphorylation model and/or parameter values. Apart from the changes explicitly listed, all other parameters were unaltered. In each case, the mean number of photoisomerizations was 0.65, and we added noise as above.

The simulated genetically modified responses shown in Figs. 4 H, 7 F, and 8 F were generated using the means of 10 model SPRs in order to approximate the number of responses averaged in the corresponding experimental studies. These simulations represent the predictions for toad rods with the genetic manipulations used in four studies on mouse rods (Xu et al., 1997; Chen et al., 1999; Mendez et al., 2000; Burns et al., 2002).

In order to simulate rhodopsin kinase knockout (RK−/−) responses (Chen et al., 1999), the parameter k_RK1[RK] was set to zero across all n in Eq. 1a. Arr−/− (Xu et al., 1997) was simulated by setting k_A[Arr] = 0 across all n. In order to simulate the transgenic disabling of all phosphorylation sites on the carboxy terminus of rhodopsin (analogous to Mendez et al., 2000, “completely substituted mutant”, or CSM, in which serine and threonine residues were replaced by alanine), we set k_RK3 = 0 across all n.

In order to simulate the ATP and GTP manipulations used to generate the data in Fig. 14 of Rieke and Baylor (1998a), we scaled the parameters k_RK3[ATP] and k_G5[GTP] shown in Table I by 0.04 and 0.4, respectively, the same factors used in Rieke and Baylor (1998a) (see Eqs.

An early saturation model was implemented by restricting the amount of PDE locally available subsequent to a photoisomerization. R* shutoff was achieved in a single step by setting k_A(1) to infinity, so that after the first phosphorylation, Arr-capping was automatic and instantaneous. This ensured that the only mechanism available to reduce SPR variability was the local depletion of PDE.

The local depletion of PDE was achieved by scaling k_P1[PDE] (Eq.

Ca²⁺ clamp was simulated by replacing the time-varying calcium variable, c, with the constant, c_dark (so that

where

In order to establish an empirical baseline for SPR variability, four analyses of SPR variability of Whitlock and Lamb's (1999) original data were performed. Two of these analyses were also presented in Whitlock and Lamb (1999), and two were new analyses not presented in their paper. The results analyzed here were those presented in their Figs. 1–3, comprising responses of a toad rod to a set of 350 identical flashes that delivered on average ∼0.6 photoisomerizations per flash. The methods and procedures used in collecting Whitlock and Lamb's original electrophysiological data are described in detail in their paper (Lamb et al., 1986; Whitlock and Lamb, 1999).

Raw data records were first corrected for low-frequency drift, then objectively tested for being at baseline at the time of the delivery of the stimulus, and then zeroed by subtraction of a DC offset using the method we now describe.

The full data set consisted of 350 recording epochs of 10-s duration each, with the 20-ms (500 nm) flash presented 1 s into the recording epoch. Since the prestimulus portion of each recording comprises one tenth of the total epoch, we assume that at least 10% of each response should be at baseline levels, and extract the 10th percentile from the sorted current values for each record to form an estimate of the baseline drift over the experiment (∼1 h). This time-varying baseline drift waveform was low-pass filtered (−3 dB at 0.013 Hz) and subtracted from the raw responses.

The means of each of the 350, 1-s prestimulus intervals from the drift-corrected data were computed, and their distribution was fit with a Gaussian probability density function. Using the parameters of the Gaussian fit, responses whose pre-stimulus means were more than three standard deviations from the overall mean were excluded from further analysis (“bad-zero” records). 36 responses were excluded on this basis. The mean of the remaining 314 prestimulus means (a scalar DC offset for the entire experiment) was calculated and subtracted from each record, ensuring that all responses used for analysis start from zero, on average.

One additional response containing a large negative artifact was manually excluded, leaving 313 dim-flash response records for analysis.

SPRs were identified from analysis of the histogram of response amplitudes. Amplitude was defined as the scaling factor providing the best (least squares) fit of the normalized mean response to each individual response. Fits were carried out over the time interval between t = 0 (time of the stimulus) and the time to peak (1.91 s) of the ensemble mean response (Field and Rieke, 2002). A sum of Gaussians model (Del Castillo and Katz, 1954; Baylor et al., 1979) was then fit to the histogram, yielding estimates of the mean number of photon absorptions, mean SPR amplitude, SPR amplitude variance, and additive background noise variance, for the underlying distribution. These four parameters were used to determine analytically a range of amplitudes where the probability that the response resulted from a single photon absorption was ≥50%. Responses were classified as either failures, SPRs, or MPRs on the basis of whether their amplitude fell below, within, or above this range. In simulations, this approach turned out to provide high sensitivity and positive predictive value (>95% for each).

To separate variability due to background noise from stochastic variability in the actual underlying response to photon absorption, we subtract the variance of responses classified as failures from the variance of those classified as SPRs. In any subsequent discussion of SPR variability, this correction is implied. This technique is used in the CV calculations for SPR amplitude and SPR area, as well as in the time-varying SPR variance plot described below.

In the first new analysis, we were interested in CV_area for the SPR (Field and Rieke, 2002). For each failure and SPR (identified by the method just described) we integrated the response from time zero until the end of the recording epoch at t = 9 s. CV_area was calculated as the square root of the variance difference between SPRs and failure areas, divided by the mean SPR area.

Another new analysis was a comparison of the time-varying ensemble variance (σ_SPR²) and squared mean (μ_SPR²), for identified SPRs only. This analysis was presented in Rieke and Baylor (1998a), and featured prominently in Field and Rieke (2002).

The other two variability measures we used had also been used in Whitlock and Lamb (1999): a histogram of dim-flash response amplitudes and calculation of CV_ampl, and a comparison of the variance of the full dim-flash ensemble (σ_dim²) with the squared mean (μ_dim²).

Fig. 3 shows the empirical data with which the model results are compared. Eight data panels are shown. Fig. 3, A–E, shows analyses of original data from Whitlock and Lamb (1999). Fig. 3 A shows the waveforms of 101 responses classified as SPRs (red), 51 responses classified as multiple photon absorptions (MPRs, green), and 161 waveforms classified as failures (gray). The solid blue curve is the ensemble mean of the SPRs. With the help of the color-coding, it is evident that the dim-flash responses are fairly well grouped, with the SPRs and MPRs being distinguishable from the failures, indicating the regularity of the elementary response.

In our analyses, we use four empirical measures of SPR variability/reproducibility to evaluate each model.

A classical method of quantifying the variability of the SPR amplitude is to plot a histogram of the dim-flash amplitudes (Fig. 3 B) and fit a statistical (sum-of-Gaussians) model (black curve) to the histogram of dim-flash response amplitudes. As described in materials and methods, we used the model fit to define amplitude boundaries with which to classify responses statistically as most likely resulting from zero, one, or multiple photon absorptions. The subset of the histogram resulting from SPRs identified by this method is shown as a red overlay. The CV of SPR amplitudes for the Whitlock and Lamb data estimated in this manner is 0.15, which is somewhat lower than some previously published values (∼0.20; Baylor et al., 1979; Schneeweis and Schnapf, 1995; Rieke and Baylor, 1998a; Whitlock and Lamb, 1999; Field and Rieke, 2002). The relatively low CV_ampl is due to the definition of amplitude we used and our method of classification of responses. Larger values would be obtained using our classification scheme if amplitude were defined as in Whitlock and Lamb (1999), or if we had estimated CV from the parameters of the sum-of-Gaussians fit. The SPRs identified from the amplitude histogram were used in the analysis of SPR area (Fig. 3 D) and SPR variance over time (Fig. 3 E).

A second measure of variability is based on comparing the light-evoked ensemble variance increase and the square of the ensemble mean (Schnapf, 1983; Rieke and Baylor, 1998a). Fig. 3 C shows this classical comparison between the σ_dim²(red) and the μ_dim² (blue) for the Whitlock and Lamb data. The variance increase was calculated by subtraction of the variance of the failures from the variance of all responses. As in Whitlock and Lamb (1999), our reanalysis of their data shows virtually the same relationship between σ_dim² and μ_dim².

This measure of SPR variability has some limitations. Whitlock and Lamb (1999) pointed out that, although a stereotyped elementary response will necessarily lead to a close match between the variance and squared mean responses, a close match between these does not necessarily imply a high degree of SPR reproducibility. This is because, for dim flashes, the variance is dominated by Poisson noise stemming from the quantal nature of light. Thus, the result in Fig. 3 C only shows that the SPR variability over time was not of sufficient magnitude to dominate the Poisson variability of the light stimulus.

In appendix A we extend Whitlock and Lamb's insight about the limitations of this analysis, and provide a quantitative interpretation of the relationship between the squared mean response and variance of the response. Despite these caveats, this analysis can help to evaluate models since failure of proportionality between σ_dim² and μ_dim² does imply nonstereotypic SPRs.

A third gauge of SPR variability is the variability of the response area for SPRs. The motivation for the area analysis comes from the fact that, as Field and Rieke (2002) pointed out, the area under the response waveform is expected to be a good gauge of overall response variability since it includes the effects of response variability throughout the entire response waveform. Apart from this pragmatic consideration, there is a theoretical basis for using the variability in response area (over other measures) in the analysis of mechanisms of SPR reproducibility. In discussion, we show why it is the CV of SPR area, and not the CV of either SPR amplitude or duration, that tracks the variability in integrated R* activity, and hence, can provide an estimate of the number of underlying R* inactivation steps.

A histogram of dim-flash response areas for the Whitlock and Lamb data is shown in Fig. 3 D. This analysis was not presented in the original Whitlock and Lamb (1999) paper. The red overlay in the histogram shows the subset of areas from the responses that had been classified as SPRs in Fig. 3 B. The CV for area for the SPRs derived in this fashion was 0.36, similar to values recently reported for mammalian rods (∼0.3, Field and Rieke, 2002). Using our sequential phosphorylation model, we found that after 1000 random additions of simulated noise to a single Monte-Carlo run of 350 trials (dim flashes), 95% of the estimates for CV of SPR area fell in the interval 0.30 to 0.44.

The fourth measure of SPR reproducibility comes from an analysis of the time-dependent residual variability of the SPRs (Fig. 3 E). This analysis was introduced by Rieke and Baylor (1998a)(see their Fig. 5) and was featured prominently in a recent paper by Field and Rieke (2002). The panel depicts the time course of the noise-corrected SPR variance (red curve; SPR variance minus failure variance) and the square of the mean SPR (blue curve). The SPRs were categorized in the same manner as for the amplitude histogram in Fig. 3 B. As noted in Rieke and Baylor (1998a) and Field and Rieke (2002), σ_SPR² is approximately an order of magnitude smaller than μ_SPR² until after the peak of the squared-mean response. In addition, σ_SPR² peaks much later than μ_SPR² (1.5 times), and is broader. Field and Rieke (2002) emphasized that these features provide constraints on models of SPR reproducibility, aiding in the discrimination between models.

Aside from the variability/reproducibility measures described above, other data in the literature provide additional constraints on any candidate model.

Fig. 3, F and G, reproduces data from Fig. 14 of Rieke and Baylor (1998a) showing the effects of lowering transduction gain in the presence of normal (500 μM) or low (20 μM) ATP levels in dialyzed toad rod outer segments. With normal ATP (Fig. 3 F), lowering GTP by a factor of 2.5 caused the dim-flash response amplitude to decrease with no effect on the response kinetics. The kinetics were slowed by the GTP manipulation, however, if phosphorylation was substantially slowed by reducing ATP (Fig. 3 G). Rieke and Baylor (1998a) interpreted these results to imply that neither phosphorylation nor arrestin-binding controlled the majority of rhodopsin's cumulative activity. We will show that this pattern of results under GTP and ATP manipulation does not rule out phosphorylation dominating R* inactivation, contrary to Rieke and Baylor's reasoning.

Implementation of a detailed stochastic model offers an opportunity to simulate genetic KO or transgenic substitution experiments targeting mechanisms expected to affect activation, inactivation or feedback in the phototransduction cascade. In the present study, we evaluate the ability of the model to predict the results from three recent experiments that genetically manipulated R* shutoff mechanisms, plus one experiment in which feedback synthesis of cGMP was disrupted. The results of these studies are reproduced in Fig. 3 H. The panel shows data obtained from mouse rods that had six major phosphorylation sites on the rhodopsin C-terminus disabled by substitution of alanine for the WT serine and threonine residues normally occurring at these sites (green: complete substitution, or CSM, rods; Mendez et al., 2000), that had Arr knocked out (red: Arr−/−; Xu et al., 1997), or RK−/− (blue: Chen et al., 1999), or GCAPs−/− (orange: Burns et al., 2002). The WT responses from each of these studies are shown in Fig. 3 H as thin curves. The WT responses were scaled to the same relative peak amplitude (1.0), but the relationship to the corresponding genetically manipulated responses in each case was not altered.

Both the RK−/− and CSM responses rose at the same rate as the WT until ∼100 ms, and continued to rise until they reached a peak ∼2 times the WT peak amplitude. The SPRs in both CSM and RK−/− rods were step-like, shutting off abruptly at highly variable times. Histograms of the duration of SPRs from CSM and RK−/− rods were approximately exponential, with time constants of 5.1 and 3.3 s, respectively. The Arr−/− responses reached approximately the same peak amplitude as the WT responses, then exhibited a partial recovery with nearly normal kinetics. When viewed on a long time scale, the mean Arr−/− responses manifested the initial, relatively rapid recovery phase, and then settled into a long, slow recovery phase (with a mean recovery time constant = 51 s; results not shown here).⁴ The GCAPs−/− responses rose with ∼WT kinetics to peak at ∼4 times the WT amplitude at ∼300 ms.

The sequential phosphorylation model dramatically reduces SPR variability relative to the expected behavior of a single-step inactivation model, bringing it to empirical values. The model yields the correct qualitative behavior in all tests, including reproduction of the Rieke and Baylor (1998a) transduction gain manipulation experiments (Fig. 4, F and G), as well as the response features from four genetic knockout and transgenic studies: three in rods in which R* inactivation mechanisms had been genetically disrupted (Xu et al., 1997; Chen et al., 1999; Mendez et al., 2000), and one in which feedback synthesis of guanylate cyclase has been disabled (Burns et al., 2002).

The model generates responses that look closely similar to the real data (compare Fig. 4 A with Fig. 3 A). The similarity between the simulated and real responses may be evaluated more quantitatively by the four variability tests, and by simulation of other experimental manipulations.

The relatively low response variability can be seen in raw model responses (Fig. 4 A), where the simulated response failures, SPRs and MPRs, can be distinguished. The distribution of amplitudes (Fig. 4 B) reproduces the behavior obtained empirically (Fig. 3 B). The subset of amplitudes classified as SPRs by the statistical method used in analysis of the Whitlock and Lamb data is shown as a red overlay. The solid blue curve in panel B depicts the histogram of the amplitudes of the actual SPRs (which, unlike in the electrophysiological data, can be identified perfectly in the simulated responses). This illustrates that our method of amplitude measurement and response classification identifies SPRs in the presence of (simulated) recording and photoreceptor noise with high statistical accuracy (sensitivity and positive predictive value both exceed 95%).

The CV for SPR amplitudes identified statistically is 0.16, nearly identical to the value we derived from the Whitlock and Lamb data (0.15; Fig. 3 B). The close match between the

The distribution of SPR areas is Gaussian like (Fig. 4 D, red overlay) with a low CV of 0.38 that nearly matches the CV_area in the Whitlock and Lamb toad rod data (0.36; Fig. 3 D). This value is somewhat larger than the only other value reported in the literature (∼0.30 for mammalian CV_area; Field and Rieke, 2002).

The CV for SPR area produced by the model (0.38) is close to the theoretical limit for an eight-step model (0.35; see Eq.

The sequential phosphorylation model also reproduces the empirical relationship between μ_SPR²and σ_SPR² (compare panel E in Figs. 4 and 3). The σ_SPR² waveform peaks at 1.6 times the time-to-peak of μ_SPR², close to the peak shift when this analysis is applied to the Whitlock and Lamb data (1.5 times; Fig. 3 E).

The sequential phosphorylation model captures the transduction gain manipulation data from Fig. 14 of Rieke and Baylor (1998a). Fig. 4 F shows the results of a simulation under the control ATP condition. The decrease in GTP by a factor of 2.5 affects transduction gain (shown in inset) without significantly altering the kinetics of the response (shown by the larger normalized curves). However, when ATP is lowered by a factor of 2.5 (Fig. 4 G), the same GTP manipulation slowed the kinetics of the response in addition to decreasing the gain (inset). This finding is significant and will be discussed further in the discussion.

The sequential phosphorylation model reproduces the salient qualitative and quantitative features of the Arr−/−, CSM, RK−/− and GCAPs−/− genetic manipulations (Fig. 4 H). Allowing for the difference in timing between mammalian and amphibian rods, the model does well at capturing these features. These simulated responses represent predicted responses if the same genetic manipulations were performed in toad rods as were done in the mouse rods.

Because G-protein competes with RK and Arr, and because of reversibility in some of these reactions (Eqs.

Here, R* activity depends on 6 rate constants, plus the phosphorylation state of R*, n. Substituting for the rate constants in Table 1, and combining Eq. 4 with Eq. 14 we obtain

which gives a theoretical rate of activation of 146 G*/s at n = 0, consistent with recent estimates in the literature (e.g., Leskov et al., 2000). From these equations, we can predict the R* activity for each phosphorylation state, and, in particular, if Arr-quench occurs between n = 6 and n = 7 phosphorylations, the R* activity will have decayed to ∼35% or 23% of its initial level, respectively. Thus, theory shows that decay of R* activity is dominated by phosphorylation in our model.

Using our model, and simulations of genetic manipulations, we can illustrate the relative contributions of phosphorylation and Arr-binding to the inactivation process. To do this, we recorded the mean G-activation rate from 1,000 simulated SPRs under three conditions, the results of which are all illustrated in Fig. 5: (a) WT (i.e., using our sequential phosphorylation model; blue curve), where both phosphorylation and Arr-capping contribute to R* inactivation. (b) CSM (green curve), depicting the case where both phosphorylation and Arr-quench are disabled. In the absence of phosphorylation, R* activity goes to a fixed, steady-state level close to the theoretical initial maximal R* activity of ∼146 G*/s (derived from Eq. 15, for n = 0). (c) Arr−/− (red curve) depicting the case where final Arr-quench is disabled, but phosphorylation is not. Thus, this curve depicts the mean reduction in R* activity due to phosphorylation per se.

The three R* activity curves in Fig. 5 demonstrate qualitatively that phosphorylation dominates R* inactivation. In order to quantify the fractional decrease of R* activity due to phosphorylation alone, we calculated the decrement in R* activity at the time of Arr-quench for each of the 1,000 Monte-Carlo trials, and then averaged these decrements. The decrease in R* activity due to phosphorylation measured in this way was 66% of the unphosphorylated activity, with Arr-quench accounting for the remaining 34% of R* activity reduction. This shows, by numerical simulation, that phosphorylation dominates R* inactivation in the sequential phosphorylation model.

These results are expected from general theoretical considerations if optimal variability reduction is to be achieved: when Arr quench of R* activity is only one final step out of eight potential inactivation steps, the multiple steps of phosphorylation preceding Arr-binding must dominate R* inactivation. Otherwise, the stochastic nature of the final Arr-binding would contribute a disproportionate amount of variability to the overall inactivation process.

We have presented analytic and computational methods with general application for testing theories of the reproducibility of the SPR of retinal rods. The computational methods are based on detailed stochastic modeling of the underlying biochemical kinetics and utilize Monte-Carlo simulations. Our emphasis was on testing the sequential phosphorylation model inspired by the biochemical experiments of Gibson et al. (2000) showing that the affinity of phosphorylated rhodopsin for G-protein declines exponentially with phosphorylation number, while the affinity of phosphorylated rhodopsin for arrestin grows linearly with phosphorylation number. In addition, as discussed below, we have simulated some experiments and analyses of Whitlock and Lamb (1999) in order to evaluate their evidence that calcium feedback plays a central role in the reproducibility of rod responses.

Previous authors have suggested that the processes of phosphorylation and arrestin binding may constitute multiple steps of R* inactivation that contribute to reducing SPR variability (e.g., Whitlock and Lamb, 1999; Gibson et al., 2000; Mendez et al., 2000; Field and Rieke, 2002). Using a full stochastic, biochemical model, our simulations and analyses demonstrate quantitatively for the first time that multiple phosphorylation of R* (plus Arr-binding) can account for the SPR reproducibility observed in vertebrate rods. Sequential, phosphorylation-dependent ratcheting down of R* activity, and ratcheting up of inactivation rate, can change the distribution of R* lifetimes from approximately exponential (with CV ≈ 1) to a much less variable distribution, substantially reducing the variability of SPR amplitude, kinetics, and area below the levels that would otherwise occur if R* were shut off in a single, memoryless step (see appendix B; compare Figs. 4 and 8).

We found that the sequential phosphorylation model accounts for the four measures of reproducibility and generates responses that exhibit the correct features in almost all details. The model also accounts for the signature qualitative features of four genetic experiments in mouse rods (see Fig. 4 H). The match to these data was achieved without any additional parameter adjustments, other than the simulation of the genetic manipulation. The model also reproduces the results of one of Rieke and Baylor's (1998a) transduction gain experiments that was thought to rule out phosphorylation (and/or Arr-binding) as the dominant mechanism in the deactivation of R*.

Rieke and Baylor (1998a) reasoned that the low GTP concentration (40% of control) in their transduction gain experiments would increase the amount of time a G-protein spent bound to R*. Because inactive G-protein, RK, and Arr are thought to compete in a mutually exclusive manner for R*, this should reduce the availability of R* for phosphorylation and arrestin capping. Consequently, if phosphorylation and Arr binding were responsible for a large portion of R* inactivation, one would expect both a reduction in amplitude and a slower response recovery in low GTP. They interpreted the result, that transduction gain was reduced and that response kinetics were unchanged, to mean that phosphorylation and arrestin binding do not make a major contribution to the inactivation kinetics of R*. They further tested this hypothesis by lowering ATP to slow phosphorylation so that it would be expected to control a significant fraction of R*'s cumulative activity. Under this condition, lowering GTP by the same amount did significantly slow the response (Fig. 3 G).

Our simulations reveal that Rieke and Baylor's nucleotide manipulation experiments do not, in fact, rule out R*'s inactivation being dominated by phosphorylation. In the sequential phosphorylation model, R* is inactivated entirely by phosphorylation followed by arrestin capping, with phosphorylation accounting for ∼66% of the total R* activity reduction, and the final Arr-quench accounting for the remaining 34% (Fig. 5). Yet, as discussed in the following section, when reversibility of some of the early reactions in the cascade is taken into account, we find that the model can reproduce the transduction gain results of Rieke and Baylor.

Rieke and Baylor's reasoning would be correct if the reactions governing the interaction between R* and G-protein were unidirectional. The reasons that reversibility in these early reactions allows the sequential phosphorylation model to reproduce the Rieke and Baylor (1998a) nucleotide manipulation results may be understood in the following way.

Reversibility of the interactions between R* and G-protein creates two possible outcomes for a G-protein molecule that has formed a complex with R*. It can dissociate from R* either in its active form, after replacement of the GDP by GTP (the catalytic route), or while still in the inactive GDP-bound form (see Eq.

Lowering the intracellular GTP concentration decreases the likelihood of G-proteins taking the catalytic route, so that a G·GDP molecule requires a greater number of R* encounters, on average, before it is activated. This, of course, reduces the rate of G·GTP formation, i.e. the gain of the response. In addition to decreasing the probability of catalysis, decreasing the GTP concentration will also increase the average time between the formation and dissociation of an R–G·GDP complex (whether by the catalytic or non-catalytic pathway). The increase in lifetime of R–G complexes will perforce slow phosphorylation and Arr-binding, due to the competitive nature of the interaction with R*.

Thus, Rieke and Baylor (1998a) were correct to conclude that, in addition to lowering the gain, the increase in lifetime of the R–G complexes would lead to a slowing of inactivation kinetics, if phosphorylation and arrestin controlled R* shutoff. However, we found that if the rate constants in Eq.

When the concentration of ATP is low, the slowing of kinetics elicited by the lowered GTP concentration will no longer be negligible. The lowered ATP concentration reduces the rate of phosphorylation of R* by RK, so that a given R* molecule will bind a substantially greater number of G-proteins in each of its phosphorylation states. Since each and every G-bound state of R* represents a timeout from the inactivation reactions, summing across a greater number of such timeouts magnifies the absolute difference in inactivation rate between the control and low GTP conditions.

When the R–G reactions are assumed to be unidirectional, the R–G complexes are unable to take the noncatalytic route, so that the GTP manipulation leading to a gain reduction will necessarily cause a comparable increase in the lifetime of the complex, and this in turn will slow the rate of phosphorylation and Arr capping, and hence the overall response kinetics. Thus, when reversibility is not included in these reactions, the only way of accounting for Rieke and Baylor's GTP/ATP manipulation data would be if R* inactivation occurred primarily by mechanisms that were not competitive with G-protein (e.g., Rieke and Baylor's putative multiple R* transitions).

For the sake of completeness we implemented such a modification to our stochastic model, and we found (as expected) that it failed to account for the results presented by Rieke and Baylor (1998a) in their Fig. 14 (unpublished data).

Whitlock and Lamb (1999) analyzed variability in the rising phase of SPRs by fitting their discrete model to each response and extracting from each fit Lamb and Pugh's (1992) amplification constant, A. They found that the CV of A was 0.13, and hypothesized that this variability in apparent transduction gain could be due to variability in the packing density of proteins on the 1,000 or so disc surfaces upon which the photoisomerizations occur across stimulus trials (see Calvert et al., 2001).

We have applied this same analysis to the SPRs generated by the sequential phosphorylation model (by fitting Whitlock and Lamb's discrete model to the simulated SPRs in the absence of recording or photoreceptor noise), and find that, as in the physiological data, the simulated SPRs do not have stereotypical rising phases—the CV for A is 0.19. Since the model does not have explicit implementation of packing density and diffusion kinetics, nor a mechanism simulating longitudinal variations in transduction gain (Schnapf, 1983), and since recording and photoreceptor noise were excluded from this analysis, this result demonstrates that a substantial proportion of the variability in photocurrent activation in the dim-flash regime may be due to the underlying stochastic variability in G-protein and PDE activation per se. The balance of contribution of these different mechanisms remains to be evaluated.

Biochemical studies in vitro indicate that seven or more phosphates per rhodopsin are incorporated (Kühn and Wilden, 1982; Wilden and Kühn, 1982; Aton et al., 1984; Thompson and Findlay, 1984; Palczewski et al., 1991; Wilden, 1995). However, there are well-known technical difficulties in determining sufficiently quickly the number and identity of the rhodopsin sites phosphorylated under dim-flash conditions in vivo. For one thing, dephosphorylation events subsequent to physiologic quench can readily cause the number of phosphorylations to be underestimated (Ohguro et al., 1995, 1996; Hurley et al., 1998; Kennedy et al., 2001; Maeda et al., 2003). Thus, in vivo studies have not been able to demonstrate the incorporation of more than one (Ohguro et al., 1995, 1996) or, at most, three or four phosphates per rhodopsin (Kennedy et al., 2001). Moreover, these studies were conducted at intensities orders of magnitude above the single-photon level, where rhodopsin kinase may well be saturated, and therefore they may not accurately reflect the incorporation of phosphates under dim-flash conditions.

A recent study by Mendez et al. (2000) has provided evidence that, under physiological conditions, all the available phosphorylation sites are required for the normal kinetics of deactivation of rhodopsin. They used transgenic techniques to substitute various serine and threonine residues in the carboxy-terminal region with alanines, and they found that the rate of recovery of the dim-flash response increased systematically as the number of phosphorylation sites available was increased, implying that all the native phosphorylation sites were needed to support normal inactivation kinetics. While this in vivo study did not actually show that each of these sites underwent phosphorylation, it did show that the presence of all the native sites is important for normal response kinetics. Accordingly, the Mendez et al. (2000) results are consistent with our assumption that all seven sites are available for phosphorylation, and furthermore we are not aware of any convincing evidence against the notion that, under single-photon conditions, as many as seven phosphorylations do indeed occur.

Mendez et al. (2000) concluded that three phosphorylation sites are necessary for normal SPR reproducibility, and we agree with this, but they further concluded that more than three sites do not further improve SPR reproducibility, even though the rate of deactivation may be decreased. However, from theoretical considerations we can show that, in the absence of contributions from other mechanisms, the number of R* inactivation steps needed to achieve the observed SPR reproducibility must be greater than three. For example, our analyses of the sequential phosphorylation scheme show that if multiple, sequential phosphorylation (with Arr capping) is the mechanism that reduces SPR variability, then three phosphorylation sites cannot be sufficient to account for the observed variability in SPR area. The lowest CV for SPR area that three phosphorylations plus Arr-capping could support is 0.50, considerably higher than the empirical values (0.30, mammalian, Field and Rieke, 2002; 0.36, amphibian, Fig. 3 D) or the value from simulations with our full model (0.38).

We will now show that there is not a straightforward relationship between the statistics of the random lifetime of a single activated rhodopsin molecule and the statistics of either dim-flash response amplitude or response kinetics alone. In particular, there is a theoretical tradeoff between variability in SPR amplitude and variability of SPR duration that depends on the kinetics of R* inactivation relative to the kinetics of reactions downstream to R*. Moreover, this tradeoff could be informative in determining what inactivation reactions might be rate-limiting in the recovery of rod responses in the dim-flash regime, since the relative CVs of SPR amplitude, duration, and area may indicate the extent to which R* is or is not rate-limiting in SPR recovery.

The variability of the photocurrent response in any biochemical kinetic scheme (in which non-linearities subsequent to PDE* production do not play a significant role) is determined by variation in the number of PDE* molecules produced, the time at which they are produced, and the variability in the lifetimes of individual PDE* molecules. In the sections below, we present a theoretical basis for the superiority of CV of SPR area as a measure of variability, and show why variability in SPR amplitude or kinetics alone are, in principle, less informative gauges of the variability.

Let N be a stochastic variable representing the cumulative R* activity during the SPR, i.e., the number of transducin molecules activated during R*'s lifetime, and hence the number of activated PDE* subunits. If SPR reproducibility were to derive purely from the regularization of R* activity by a sequence of n + 1 inactivation steps (e.g. n phosphorylations plus arrestin), then under optimal conditions, the lowest attainable CV_N (i.e., σ_N/μ_N) would be given by

Below we show that CV_area of the SPR provides a good approximation of the CV_N, and hence can be used to estimate a lower bound on the number of inactivation steps of R*.

Consider the net PDE* response to a flash to be the superposition of N discrete PDE* responses, each with unit amplitude, random onset time, and random (exponentially distributed) lifetime. The CV of the area under the PDE* waveform (A_PDE) will be determined entirely by the statistics of N. Assuming a linear dim-flash response subsequent to PDE* production, CV_area for the SPR will equal CV(A_PDE), which we show in appendix C (Eq. C5) can be written as

The first term under the square root can be thought of as the component of CV_area due to variability in the number of PDE* activations per R*. The second term under the square root represents variability, for a given mean number of PDE*, due to stochastic variation in individual PDE* lifetimes. In practice, we have empirical estimates for CV_area and μ_N, whereas the variability in N has not been directly measured, but can be estimated by a rearrangement of Eq. 17.

Field and Rieke (2002), and our present analyses, find CV_area in the range of 0.26–0.35. Since the number of PDE* produced by a single activated R* during its lifetime is expected to be on the order of hundreds (e.g., Yee and Liebman, 1978; Heck and Hofmann, 2001), we expect the term μ_N⁻¹ to be at least an order of magnitude smaller than CV_area². Consequently, CV_area itself provides a close approximation to the variability in integrated R* activity, and therefore the minimum required number of inactivation steps (or, in this case, n phosphorylations) can be estimated from

To illustrate how variability in R* activity could manifest as variability in either the amplitude or duration of the response, we now consider two opposing, limiting-case scenarios using, for the sake of clarity, the following simple model: R* activity is assumed to be a rectangular pulse of fixed height and variable duration, and the SPR depends linearly and deterministically on this pulse. In the first scenario, the kinetics of the downstream cascade are assumed to be very fast compared to the typical R* inactivation rates. Here, the response is able to track the R* activity function with little lag, and the photocurrent response itself will also approximate a rectangular pulse. In the limiting-case, the SPR is just a scaled version of the R* activity, so there is no variation in SPR amplitude, and SPR duration will be equal to, and thus have the same variability as, the lifetime of R*.

In the second scenario, R* inactivation is very fast compared to downstream kinetics, so that with respect to the timescale of the response, R* activity is effectively an impulse. Here, the amplitude of the SPR will scale linearly with the random lifetime of R*. In particular, the coefficient of variation of the SPR amplitude will be the same as that for the random R* lifetime.

For example, if R* were shut off abruptly, following a series of x nonactivity-changing transitions with equally distributed waiting times (and there were no other mechanisms to reduce SPR variability), then the CV of R*'s lifetime would equal

In order to illustrate this tradeoff within the context of a full stochastic biochemical model, we ran Monte-Carlo simulations using seven different values of τ_PDE, ranging from 0.75 s (R* recovery highly rate-limiting) to 12 s (PDE* recovery highly rate-limiting). The resulting CVs for SPR amplitude, duration, and area are shown in Fig. 6 A. The behavior of the full biochemical model clearly manifests the tradeoff behavior predicted from the limiting case scenarios: CV_ampl (open squares) is lowest when R* inactivation is rate-limiting, and monotonically increases as the downstream reactions become rate-limiting; conversely, CV_dur (open triangles) is high and decreases monotonically as the rate-limitation switches from R* to downstream recovery kinetics. However, the corresponding CV_area (open circles) does not depend on what rate-limits recovery (CV_area is constant over all τ_PDE values).

The CVs of SPR amplitude and SPR duration are not adequate measures of SPR variability, and cannot be used to estimate the number of R* inactivation steps

The analysis above (see Eqs. 17 and 18) shows that for a given mean and variance of integrated R* activity, there will be a unique value of CV_area. However, as illustrated by the two tradeoff scenarios, neither the CV of SPR amplitude nor SPR duration are predetermined by the statistics of G* production per se; depending on the (average) speed of R* inactivation in relation to the kinetics of the downstream reactions, they may take on any value from near-zero to CV_area (which may be as high as ∼1 under some conditions, e.g., RK−/−).

Thus, without a priori knowledge of the relevant reaction rate constants, it is inappropriate to use either SPR amplitude or duration variability to estimate the minimum number of R* inactivation steps. This is illustrated in Fig. 6 B, where we have plotted the inferred number of R* inactivation steps corresponding to the CVs in Fig. 6 A. Only the number of steps inferred from the CV_area measures (open circles) match the actual number of R* inactivation steps used (X, dashed line).

In general, CV_ampl and CV_dur measures will tend to overestimate the number of R* inactivation steps, and may do so severely. For example, a typical empirical CV_ampl of the single-photon response is ∼0.2 (Baylor et al., 1979; Schnapf, 1983; Baylor et al., 1984; Schneeweis and Schnapf, 1995; Rieke and Baylor, 1998a; Whitlock and Lamb, 1999; Field and Rieke, 2002), which would correspond to an estimated minimum of 25 inactivation steps. The present analyses provide another striking example of the fallacy of this approach. The value of CV_ampl obtained from the sequential phosphorylation model was 0.16, which would correspond to an estimate of 39 inactivation steps, far greater than the number of steps actually used in the model. The value of CV_area, however, was 0.38, commensurate with the ∼7 inactivation steps that were used (at Arr-binding, an average of 6.1 phosphorylations occurred, plus Arr-binding = 7.1 steps, corresponding to a CV_area of 0.38; Fig. 4, B and D).

In practice, measures of CV of SPR amplitude and SPR duration have been reported to be of similar magnitudes, in the ranges 0.2–0.25 and 0.2–0.4, respectively (Baylor et al., 1979; Schnapf, 1983; Baylor et al., 1984; Schneeweis and Schnapf, 1995; Rieke and Baylor, 1998a; Whitlock and Lamb, 1999; Field and Rieke, 2002). These findings suggest that the rates of R* inactivation and downstream kinetics are not highly disparate in normal vertebrate rods.

Four basic classes of models have been introduced to account for SPR reproducibility, and they are not mutually exclusive: late saturation, early (local) saturation, feedback, and multistep inactivation of R*. As we now discuss, none of the first three classes alone appear to be viable as the primary mechanism underlying SPR reproducibility.

Single-photon responses might be rendered insensitive to variability in R* lifetime due to localized depletion of either cGMP or the number of available open cGMP-gated membrane channels. Some data from Field and Rieke (2002) argue against both of these possibilities in mammalian rods. They found that responses from local stimulation of the outer segment (1–2 μm, corresponding to ∼40 discs) were indistinguishable from responses obtained using diffuse illumination of the entire outer segment, and that, in both conditions, responses summed approximately linearly up to 3 R*. Photocurrent waveforms did not indicate any interaction between responses, even for two or three photo-isomerizations occurring in the same region of the outer segment. Since the PDE*s elicited by multiple photon absorptions within the local stimulation area were arguably competing for the same overall pool of cGMP (and cGMP-gated channels), these experiments indicated that activation of a single R* in normal (diffuse-light) conditions does not tax the limits of available channels or cGMP.

It is possible that the situation is different in rods of some lower vertebrates since Lamb et al. (1981) found that, even at very low intensities (≤3 R*/flash), restricted stimuli generated smaller responses than diffuse stimuli eliciting an equal number of photoisomerizations in toad rods. However, Rieke and Baylor (1998a) provided evidence against late saturation in amphibian rods by showing that manipulations that increased the amplitude of the SPR (e.g., lowering ATP to slow phosphorylation, or clamping Ca²⁺) did not generate responses with any evidence of saturation. Responses showed no clipping or other nonlinear behavior. In addition, when they lowered transduction gain by decreasing GTP (and, in some cases, also increasing GDP), the responses had the same shape as control responses when scaled to match in amplitude. If late saturation had played a significant role in shaping the responses in normal gain conditions, then the reduction in gain should have resulted in differently shaped responses (i.e., less distorted by saturation).

In principle, responses could also be rendered insensitive to variability in integrated R* activity due to depletion of an early transduction intermediate (PDE or G-protein). This appears to be ruled out as a main mechanism by empirical data from at least four studies (Xu et al., 1997; Chen et al., 1999; Mendez et al., 2000; Rieke and Baylor, 1998a), by the theoretical analyses in Field and Rieke (2002), and by our own analyses and simulations, as explained below.

First, a model in which SPR reproducibility is achieved solely by R* running out of G-protein (or G-proteins running out of PDE) in a local neighborhood on the disc surface would fail to account for the data from genetic experiments in which R* phosphorylation is disrupted. Assuming that in Arr−/−, RK−/−, and CSM rods rod structure is normal, and other proteins are unaffected by the genetic intervention, one would expect the same local saturation to occur in these rods as in WT rods. Thus, if saturation were responsible for the bulk of time-dependent effective R* activity reduction, one would expect the SPRs from Arr−/− (Xu et al., 1997), RK−/− (Chen et al., 1999), and CSM (Mendez et al., 2000) rods to exhibit recovery with kinetics similar to WT, contrary to what is observed (Fig. 3 H).

Second, the transduction gain manipulation experiments of Rieke and Baylor (1998a) argue against depletion of either G or PDE as the main variability-reducing mechanism for the same reasons as cited in the section above.

Third, in their theoretical analyses, Field and Rieke (2002) showed that local depletion of PDE predicts an incorrect relationship between the time to peak of the SPR variance and SPR mean (σ_SPR² peaks too early), and predicts a waveform for SPR variance that is much less broad than the observed waveform.

Finally, we have implemented a stochastic biochemical model in which early saturation (local depletion of PDE) was the sole mechanism that reduced SPR variability (Fig. 7 ; see materialsandmethods for model details). We systematically reduced the total pool of locally available PDE until the CV_area (0.32; Fig. 7 D) and CV_ampl (0.18; Fig. 7 B) matched empirical levels, as well as the CVs obtained from the full sequential phosphorylation model (Fig. 4). In addition, the variance of the ensemble dim-flash responses was close to the square of the ensemble mean (Fig. 7 C), and the individual simulated responses (Fig. 7 A) looked similar to the real data (Fig. 3 A).

In these respects, the performance of the model is not readily distinguishable from the performance of the full sequential phosphorylation model (compare with Fig. 4, A–D). However, as predicted, the early saturation model fails completely to reproduce the correct Arr−/−, RK−/−, and CSM response features; these all recover fully with nearly WT kinetics (Fig. 7 F), in striking contrast to the data (Fig. 3 H). In addition, the predicted SPR variance peaks too early and is too narrow (compare Fig. 7 E with Fig. 3 E), confirming the analysis of Field and Rieke (2002).

The early saturation model also produces an inordinate number of small-amplitude, small-area SPRs that overlap substantially the distribution of response failures. This is not evident if the raw model responses are analyzed using the same statistical methods used to analyze real data (red overlays in Fig. 7, B and D). However, when we plot the distribution of veridical SPR amplitudes and areas (solid blue curves, Fig. 7, B and D), the substantial overlap with the corresponding distributions of response failures (gray distributions centered at zero on the abscissas of Fig. 7, B and D) can be seen. The variability of these veridical distributions is much higher than the CVs derived when the responses are analyzed as empirical data (CV_ampl = 0.45 vs. 0.18; CV_area = 0.57 vs. 0.32). More detailed discussion of this issue is presented in appendix B, where we show the results of implementing a model with a 1-step inactivation of R* and no other variability reducing mechanisms.

Feedback regulation of R* lifetime or activity, whether mediated by intracellular Ca²⁺ or some other messenger, seems to be ruled out as the main mechanism of SPR reproducibility by several lines of evidence. The evidence provided by Rieke and Baylor's (1998a) Ca²⁺-clamp and high-gain/low-gain experiments argue against any kind of feedback messenger (including Ca²⁺) being the sole mechanism for SPR reproducibility. Field and Rieke (2002) examined the predictions of both Ca²⁺-mediated and non-Ca²⁺–mediated feedback models. The non-Ca²⁺–mediated model predicted SPR variance waveforms that were significantly more narrow than their data or than the predictions from a multi-step model. Moreover, the SPR variance waveforms peaked too soon relative to the SPR mean waveforms.

In addition to the experiments and analyses by Rieke and colleagues just cited, several other lines of evidence argue against Ca²⁺-mediated feedback as the main mechanism reducing R* lifetime variability. Briefly, some of the major arguments against a Ca²⁺-feedback regulation of SPR variability are:

(a) Burns et al. (2002) found that SPRs from GCAPs−/− mouse rods were nearly unaffected by application of the strong Ca²⁺ buffer, BAPTA. The logic here is that, in the absence of GCAPs, the main sites left for Ca²⁺ feedback are recoverin-RK and Ca²⁺ effects at the cyclic nucleotide gated channel (Hsu and Molday, 1993, 1994; Gordon et al., 1995; Sagoo and Lagnado, 1996; Haynes and Stotz, 1997; Rebrik and Korenbrot, 1998; Rebrik et al., 2000). Since they found that BAPTA had no effect on the SPR amplitude or kinetics when cyclase feedback was taken out of the picture, the implication was that any Ca²⁺ feedback onto recoverin-RK in the GCAPs−/− rods was not strong enough to affect the responses in the single-photon regime. That such feedback does exist and is functional in GCAPs−/− rods was suggested by their finding that GCAPs−/− responses showed evidence of Ca²⁺-mediated effects at high intensities.

(b) At least two theoretical analyses argue against a Ca²⁺-feedback model. First, Field and Rieke (2002) simulated the effects of early Ca²⁺-feedback and found that the predicted σ_SPR² waveform peaked too late relative to the μ_SPR² in comparison to the empirical data or the predictions of their multi-step model. Second, we have simulated Ca²⁺-feedback onto R* lifetime via recoverin-RK with a deterministic version of the sequential phosphorylation model. We find that SPR amplitudes and kinetics are nearly unaffected by the Ca²⁺ changes induced in the single-photon regime, even when local Ca²⁺ changes are simulated by assuming that all the change in internal Ca²⁺ concentration is restricted to only 10% of the OS volume (Lamb et al., 1981; Whitlock and Lamb, 1999; unpublished data).

(c) The sequential phosphorylation model reproduces some key results from Whitlock and Lamb (1999) that were interpreted as evidence for Ca²⁺-mediated regulation of R* lifetime variability (via recoverin-RK). They fit a single-photon, discrete version of the Nikonov et al. (1998) model to their SPR waveforms (Eqs. 2 and 3 in Whitlock and Lamb, 1999) collected from control rods as well as from rods infused with BAPTA. The variability in their model was governed by a parameter, t_life, the stochastic lifetime of (all-or-none) R* activity. Whitlock and Lamb (1999) found that both CV_ampl and the mean of the distribution of t_life values increased under BAPTA, suggesting (within the context of their model) that the slowing of Ca²⁺ dynamics had increased R* lifetime variability.

There is an alternative interpretation of the results of the Whitlock and Lamb's (1999) BAPTA experiments that derives from the tradeoff between amplitude and duration variability discussed in the previous section. The effect of BAPTA (or full Ca²⁺-clamp) would be to slow the net kinetics of the downstream reactions relative to R* inactivation by slowing (or stopping) the feedback synthesis of cGMP. Our tradeoff analysis would predict that, in addition to increasing mean t_life, BAPTA (or Ca²⁺-clamp) would increase CV_ampl, but would decrease the CV of t_life (a measure SPR duration). This is, in fact, what Whitlock and Lamb (1999) observed. We confirmed this prediction in additional simulations with the sequential phosphorylation model (which has no Ca²⁺-feedback onto R* lifetime or activity). When full Ca²⁺-clamp was applied to the model, both CV_ampl and mean t_life increased, but the CV of t_life decreased (unpublished data).

The above interpretation provides a mechanistic explanation for Field and Rieke's (2002) caution that “…the inferred increase in rhodopsin lifetime [by Whitlock and Lamb] in BAPTA-loaded rods could have…been caused by a slowed Ca²⁺ feedback to cGMP synthesis (p. 743).”

(d) Some empirical evidence that Ca²⁺ is not affecting R* variability comes from Field and Rieke (2002), who found that the CV of SPR area in BAPTA-loaded rods was indistinguishable from that measured under control conditions. We have shown in the present paper that the CV of SPR area provides a gauge to underlying variability of integrated R* activity (Eqs. 17 and 18). If calcium feedback were contributing significantly to the reduction of variability of R* shutoff, then slowing down this feedback signal with BAPTA should have made single-photon responses more variable, and increased the CV_area, but it did not. The sequential phosphorylation model also reproduces this result (the CV of SPR area is unchanged under Ca²⁺-clamp; unpublished data).

Models based on multiple deactivation steps seem to the only viable remaining class of model.

Rieke and Baylor (1998a) applied an array of experimental tests which, along with some modeling, they argued provided strong evidence against any form of saturation or feedback as the sole mechanism conferring SPR reproducibility. This left multistep rhodopsin shutoff as the only remaining viable mechanism. Although Rieke and Baylor (1998a) acknowledged that the known processes of R* inactivation—phosphorylation and Arr-capping—would contribute to recovery and hence to the overall response variability, their experiments in which nucleotide levels were manipulated led them to conclude that these two processes could not dominate R* inactivation, a conclusion that our results directly dispute (see Fig. 4, F and G, and conclusions about the sequential phosphorylation model above). Their rejection of phosphorylation as the dominant inactivation mechanism (in conjunction with their arguments against saturation or feedback) is what led Rieke and Baylor (1998a) to propose that 10–20 phosphorylation-independent intrinsic conformational changes of R* occur to alter its catalytic activity, even though evidence for these steps in the known biochemistry and biophysics of rhodopsin is lacking. Because they did not have an explicit stochastic biochemical model, they could not test their hypotheses regarding the role of phosphorylation and Arr-binding, or the role of putative phosphorylation-independent R* states.

For the sake of completeness, we implemented a model that incorporated Rieke and Baylor's (1998a) putative conformational R* transitions within the context of a stochastic, biochemical model structure. In order for this model to capture dim-flash response behavior, as well as the Rieke and Baylor (1998a) gain manipulation data as well as the knockout and transgenic data, we had to make a number of assumptions. First, we assumed that the putative rhodopsin transitions required a phosphorylation step, and would only occur after phosphorylation was complete. If phosphorylation and the transitions proceeded concurrently, the model could not reproduce the qualitative features of the RK−/− or CSM data; simulation of these responses would exhibit recovery not observed in the data. We implemented the model using only one phosphorylation step prior to the R* transitions, in order to capture the scenario suggested by Rieke and Baylor (1998b).

In order to be able to capture the Arr−/− data, we assumed that, in WT rods, after phosphorylation and the R* transitions were complete, there still remained ∼25% of R*'s activity. In this case, when Arr−/− was simulated, the model could exhibit the partial recovery to ∼50% baseline seen in the data. If, on the other hand, we assumed that in WT rods, each of the 10–20 total inactivation steps (phosphorylation, R* transitions, and Arr included) ratcheted down R* activity by an equal amount, the simulated Arr−/− responses would recover nearly to baseline, in contrast to the data.

Under these constraints, we verified that a Rieke/Baylor-like multistep model can reproduce the observed properties of empirical SPRs, including the four variability measures, the transduction gain experiments of Rieke and Baylor (1998a), and the transgenic and KO data examined here (unpublished data). Despite this success, the model is not grounded in well-established biochemistry; there is no experimental evidence for a long sequence of postphosphorylation molecular transitions of rhodopsin. In contrast, the sequential phosphorylation model has substantial experimental (and now, theoretical) support.

The sequential phosphorylation model is consistent with a broad array of electrophysiological and biochemical data (compare Figs. 3 and 4), and reproduces all the empirical results examined, including the nucleotide-manipulation experiment of Rieke and Baylor (1998a) that led them to reject such a scheme. There is thus no compelling need to invoke a long series of non-phosphorylation dependent state changes of R* in order to account for SPR reproducibility.

One of the most difficult and enduring problems in the history of research on photoreceptor physiology has been to understand the single-photon response, and in particular its relatively high degree of reproducibility (Pugh, 1999; Pugh and Lamb, 2000) despite the high variability inherent in all molecular reactions. Such reproducibility in the single-photon regime imposes strong constraints on any model of phototransduction.

We introduce a stochastic biochemical model, the sequential phosphorylation scheme, that is able to account for all aspects of observed SPR behavior examined here (Fig. 4), including four quantitative measures of reproducibility, and the transduction gain manipulations of Rieke and Baylor (1998a). The latter results and associated analyses call into question Rieke and Baylor's rejection of phosphorylation and arrestin-binding as the dominant mechanisms of R* inactivation, and suggest the testable biochemical hypothesis that one or more of the early reactions between R* and G-protein has a significant reverse component.

The paper also utilizes the analysis and simulation of knockout and transgenic experiments as an invaluable tool for testing and rejecting candidate models. The sequential phosphorylation model was able to reproduce salient response features from rods either with various R* inactivation mechanisms genetically disabled (RK−/−, CSM, Arr−/−), or with the mechanism of feedback synthesis of cGMP disrupted (GCAPs−/−). The simulations of the knockout and transgenic data illustrate vividly how some models may be able to reproduce the empirical levels of SPR reproducibility (CV_area), but may fail completely to reproduce fundamental qualitative features of a suite of extant (or simulated) knockout data, and can thus be rejected (e.g., see Fig. 7).

The sequential phosphorylation model has the virtue that it achieves all this while being based on a broad range of biochemical evidence. We believe that the model provides a solid foundation for further development and testing of models. It is clear that the phototransduction cascade is complicated enough so that stochastic modeling techniques of the sort used in this paper will play a central role in the evaluation of future hypotheses and data.

Whitlock and Lamb (1999) have pointed out that, although a stereotyped elementary response will necessarily lead to a close match between the variance and squared mean responses, a close match between these does not necessarily imply a high degree of SPR reproducibility. This is because, for dim flashes, the variance is dominated by Poisson noise stemming from the quantal nature of light.

This idea can be stated in a quantitative manner as follows. The variance of the dim-flash response is related to the mean of the dim-flash response, the variance of the single-photon response, the mean of the single-photon response, and the mean number of photoisomerizations (λ), by

Note that the right side of the equation represents the scaling between μ_dim² and σ_dim², and the term added to 1 in parentheses (⁠

Baylor et al. (1979), Schnapf (1983), and Rieke and Baylor (1998a) assumed that the proportionality between σ_dim² and μ_dim² implied that the SPR had a stereotypic waveform, corresponding to an assumption that σ_SPR²=0 in the equation above (and, hence CV_SPR²

Whitlock and Lamb's (1999) insight may also be recast in terms of our equation: although the CV_SPR²

Our analysis extends the theoretical understanding of the problem by showing that, in principle, the CV_SPR(t) could be of any magnitude as long as it was not time-varying, and the proportionality between σ_dim² and μ_dim² would still hold.

In practice, data in the literature show that CV_SPR(t) is time varying, but is small until late in the response (Field and Rieke, 2002; also, see Figs. 3 E and 4 E). Thus, the ensemble variance and ensemble squared mean tend to be closely matched, at least up to the peak of the squared ensemble mean response (Schnapf, 1983; Schneeweis and Schnapf, 1995; Xu et al., 1997; Rieke and Baylor, 1998b; Whitlock and Lamb, 1999; Mendez et al., 2000).

For more than two decades it has been understood that, without other mechanisms to regularize responses, the observed variability/reproducibility of SPRs was inconsistent with single-step inactivation of a single R* molecule (Baylor et al., 1979) given the stochastic nature of biochemical reactions. It is useful to show the consequences of single-step R* inactivation within the context of a detailed stochastic phototransduction model. Some of these consequences can be predicted a priori, but others cannot.

Fig. 8 shows the behavior of a stochastic model in which R* was forced to inactivate in a single step. This was achieved by setting k_A(1) to infinity, so that after the first phosphorylation, Arr-capping was automatic and instantaneous. The result was that R* had fixed activity with approximately exponentially⁵ distributed lifetimes. In all other respects the model was identical to our implementation of the sequential phosphorylation model, including the simulation of recording and photoreceptor noise.

As expected, the responses exhibit a high degree of variability that is evident by inspection of the waveforms in panel A. The responses contain many long-duration responses never seen in normal rods. The variability of response waveform manifests as a severe mismatch between the σ_dim² and the μ_dim² (Fig. 8 C), and an inordinately high SPR variance (Fig. 8 E).

We know that the R* lifetimes associated with these responses must have an approximately exponential distribution, with a CV of ∼1. In this case, where R* activity is an all-or-none pulse, we have shown that R* lifetime is proportional to the number of PDE*s activated, and hence to SPR area as well (appendix C). Consequently, the variability of the distribution of R* lifetimes must be reflected fully in the analogous distribution of SPR areas. We can thus predict in advance that, in the absence of noise, the CV of SPR area will be ∼1. This is confirmed by simulations without added noise, and where the SPRs were identified perfectly based on a priori knowledge from the model. The model generates an approximately exponential distribution of R* lifetimes with CV ≈ 1, which leads to the expected approximately exponential distribution of SPR area with CV of ∼1 (not depicted).

However, when recording and photoreceptor noise and failures are simulated, and the data are analyzed as was done for the Whitlock and Lamb (1999) data, the resulting distribution of SPR areas is not exponential, and CV_area is less than one (0.77; red overlay, Fig. 8 D). This result is not due to the added noise affecting the intrinsic variability of the SPRs. Rather, it is due to the empirical method of classification of responses (see materials and methods) when we treat the model responses as real data in which one does not know in advance which responses are SPRs, MPRs, or failures. It turns out that this method leads to a biased classification of responses when one or more of the true underlying distributions of response amplitude categories is markedly asymmetric. The effect of this bias may be seen by comparing the amplitude histograms for the two cases, one in which the empirical classification method was used (red overlay in panel B), and one in which responses were classified perfectly (solid blue curve in B). The empirical classification identified SPRs incorrectly, so that the histogram of SPR amplitudes appeared Gaussian-like, with a low CV of 0.24. The CV of response areas calculated from this group of SPRs was also biased low (0.77; compare blue curve with red histogram in panel D). The group of actual SPR amplitudes (blue curve in B) is highly asymmetric and non-Gaussian, and when CV_area is calculated from these responses, it is close to the expected value of ∼1 (0.97; blue curve, panel D). The rounded onset of the area histogram is a result of the addition of noise to the model; without noise, the histogram is approximately exponential (not depicted).

We checked to see how much this factor might have biased the calculation of CV_area for the responses from the sequential phosphorylation model. In this case, the underlying distribution of response amplitudes was relatively symmetrical. Hence, the derived CV_ampl and CV_area did not change very much when perfect response identification was used instead of our empirical classification method (the CV_area increased modestly from 0.38 to 0.42). To the extent that the sequential phosphorylation model is a good representation of the underlying biology, this result reassures us that the empirical method did not significantly bias the estimate of CV_area from the Whitlock and Lamb data (Fig. 3 D).

From our understanding of the role of phosphorylation and Arr-binding in R* inactivation, we can predict the results of simulating the three genetic experiments (CSM, RK−/−, Arr−/−). In particular, since the single shutoff step is achieved by a single phosphorylation/Arr-binding event, we expect the simulated Arr−/− responses to parallel the simulated RK−/− and CSM responses. This, again, is confirmed by our modeling. Fig. 8 F shows that the simulated Arr−/− responses (red) peak at about twice the WT amplitude, and do not exhibit the partial recovery seen in the data (Fig. 3 H). In addition, although the simulated RK−/− and CSM responses are similar to the behavior seen in the data, a model in which R* terminates in a single step would not be able to reproduce any of the intermediate cases examined in the Mendez et al. (2000) experiments (where only one or two or three phosphorylation sites had been transgenically disabled).

Other behaviors cannot be predicted in advance. For example, due to the tradeoff that we now know exists between amplitude and duration variability, we cannot know a priori how the variability in R* lifetime will parse out between these two domains: we can only make a clear prediction about the variability of SPR area in the absence of noise, since this measure is not dependent on the relative “front-end” and “back-end” inactivation kinetics. In addition, the quantitative details of how response variability would manifest over time in the ensemble and SPR variances was not known in advance.

We assume that activated transducin molecules activate PDE* subunits with one-to-one stoichiometry, and define N as the random number of transducin or PDE activations per SPR, i.e., the integrated R* activity. We also assume that the active lifetimes of G-PDE* complexes are exponentially distributed with a mean of τ_PDE. The area of a PDE* response (A_PDE) is determined by summing the N random PDE* lifetimes elicited by the R*. For a given value of N, the expected value and variance of A_PDE are given by

When variability in N is taken into account, the overall expected value and variance of A_PDE are

Combining Eqs. C3 and C4, the CV of A_PDE can be derived

where CV_N and μ_N are the coefficient of variation and mean of N, respectively.

We thank Terry Hegarty, M.S., for her assistance in preparation of the manuscript.

Supported by NEI Grant EY115-13 and Smith-Kettlewell Eye Research Institute Fund 5809-0100 (R.D. Hamer), and an ARC Federation Fellowship to T.D. Lamb.

Olaf S. Andersen served as editor.