|
|
|
|
|
|
|
||
Address correspondence to R.D. Hamer, Smith-Kettlewell Eye Research Institute, 2318 Fillmore Street, San Francisco, CA 94115. Fax: (415) 345-8455; email: russ{at}ski.org
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONAPPENDIX A: APPENDIX B: APPENDIX C: REFERENCES |
|---|
Key Words: single-photon response reproducibility • stochastic model • multiple phosphorylation • phototransduction • rod responses
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONAPPENDIX A: APPENDIX B: APPENDIX C: REFERENCES |
|---|
). 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.
Evidence for Reproducibility
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 (CVampl) 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 (CVarea 
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.
A New Stochastic Model of Phototransduction Accounts for the Variability/Reproducibility of Single-photon Responses, as well as Other Key Electrophysiological Data
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, Ca2+-feedback, and Ca2+-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.
Theory for Mechanisms Underlying SPR 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.
|
MATERIALS AND METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONAPPENDIX A: APPENDIX B: APPENDIX C: REFERENCES |
|---|
and 
; Fig. 1)
, were modeled as discrete signals (time-dependent integers) due to the relatively small numbers of molecules involved. Simulated PDE* responses to single photoisomerizations were generated using a Monte-Carlo approach (see Simulations and Monte-Carlo Methods below).
|
|
, 
, and 
) that G·GDP, RK, and Arr bind competitively with R* (Pfister et al., 1983; Miller and Dratz, 1984; Buczylko et al., 1991; Pulvermuller et al., 1993; Krupnick et al., 1997; Gibson et al., 2000). R* is assumed to undergo multiple, sequential phosphorylation (Ohguro et al., 1993, 1994, 1995, 1996) at up to seven phosphorylation sites on the carboxy terminus of rhodopsin (Kühn and Wilden, 1982; Wilden and Kühn, 1982; Aton et al., 1984; Thompson and Findlay, 1984; Palczewski et al., 1991; Wilden, 1995; Eqs. 
). These reactions lead to a progressive reduction in R* activity, defined as the rate at which molecules of G* are activated per R*. Inactivation of R* is defined as the reduction in this activity due either to phosphorylation or arrestin binding, or any other mechanism that reduces the rate of G-protein activation by R* (e.g., feedback or local depletion of G-protein). R* inactivation is achieved by multiple phosphorylation according to the following equations, with the final, complete quench of R* activity occurring upon Arr binding to phosphorylated R*:
The notations "pre" and "post" in Eqs. distinguish RK-bound states of R* before and after phosphorylation.
Transducin is assumed to be activated by a conventional series of reactions (Eqs. 
; see Lamb and Pugh, 1992).
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.
 | (4) |
 | (5) |
The exponential rate, 
, was set to 0.6 from an exponential fit to the data in Fig. 2 A in Gibson et al. (2000).
Choice of Phosphorylation Dependence for R–RK Affinity
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 kRK1.
 | (6) |
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.
 | (7) |
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. 
) would cause, on average, an equal number of PDE* to be produced in each phosphorylation state, and it can be demonstrated analytically that this would minimize variability in the number of activated G-proteins produced per photoisomerization. In our implementation of the model, Arr is permitted to bind with increasing probability after the first phosphorylation (Fig. 4 D, inset). Consequently, derivation of the optimal phosphorylation-dependence for R–RK affinity is quite complicated. However, we have estimated the ideal affinity profile by numerical optimization techniques, and have verified that, with model parameters that provide a good account of the data, the affinity profile embodied in Eq. 6 yields nearly ideal minimization of variability.1
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.
Activation and inactivation of PDE (Fig. 1 B).
For the activation of PDE*, we assumed that activated transducin binds to the 
subunit of PDE (Eq. 
), and that the inhibition of cGMP hydrolysis by the PDE subunit is then relieved (Eq. 
). For the inactivation of PDE*, we did not include explicit equations for acceleration of transducin (and, hence, PDE) inactivation due to interactions with the proteins RGS9-1 and Gß5 (Skiba et al., 2001; Lishko et al., 2002; Martemyanov and Arshavsky, 2002; for review see He and Wensel, 2002). Instead, we assumed that these reactions were fast, and we incorporated their effect into a single step of PDE–inactivation (Eq. 
).2 We have examined the effects of including such reactions, and found the effects on the simulations to be negligible.
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.
Deterministic model of "back-end" reactions (Fig. 2)
.
The front-end steps contribute considerable amplification, and Leskov et al. (2000) report that a single amphibian rod R* molecule activates 150 PDE*/s.3 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*.
|
). These steps are: Ca2+-sensitive synthesis of cGMP by guanglyl cyclase and hydrolysis of cGMP by PDE* (Eq. 9); Ca2+ influx through cGMP-gated membrane cation channels and Ca2+ efflux via the Na+/Ca2+-K+ exchanger (Eq. 10); and sequestration and release of Ca2+ by intracellular buffers (Eq. 11).
 | (9) |
 | (10) |
 | (11) |
· GTP], g = [cGMP], c = [Ca2+], and cb 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; Kc, the concentration of Ca2+ at which cGMP is synthesized at half its maximal rate; m, the Hill coefficient for the action of Ca2+on cyclase; ßdark, the dark rate of cGMP hydrolysis; ßsub, the rate constant of cGMP hydrolysis per activated PDE subunit; fCa, the fraction of circulating current carried by Ca2+; F, the Faraday constant; vcyto, the effective cytoplasmic volume; Jdark, the dark circulating current; gdark, the steady-state concentration of cGMP in darkness; ncg, the Hill coefficient for channel opening by cGMP; Ca, the rate constant of calcium extrusion by the Na+/Ca2+-K+ exchanger; c0, the minimum intracellular Ca2+ concentration; k1 and k2, the association and dissociation rate constants of Ca2+ with intracellular buffers; and eT, the total calcium buffer concentration. The units and values of all parameters are given in Table I.
 | (12) |
We assumed that Ca2+ 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 Ca2+ 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.
Model Parameters
The values of the parameters for the sequential phosphorylation model (Eqs. 
) are given in Table I. Parameters were set at modern literature estimates when these were known. In general, all parameter values for which estimates were available in the literature were within a factor of two of the literature values. A dark internal calcium concentration of 500 nM was assumed (Lagnado et al., 1992; Gray-Keller and Detwiler, 1994; Sampath et al., 1998).
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. and and 
), based largely on the sequence delineated by Lamb and Pugh (1992). While empirical estimates were not available for each of these parameters, limits on some of these were estimated in Lamb and Pugh (1992) (see legend to Table I for details). The front-end rate constants and the rate constant of cGMP hydrolysis per PDE hydrolytic subunit, ßsub, were adjusted as needed to achieve a close fit (by eye) to the ensemble mean of the Whitlock and Lamb (1999) data, to generate a target number of PDE* activations per photoisomerization (200–300), as well as to achieve the appropriate qualitative behavior in the ATP/GTP manipulation experiments of Rieke and Baylor (1998a).
The back-end parameters Kc, ßdark, fCa, k1, k2, and eT 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.
Simulations and Monte-Carlo Methods
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 < kG1
kG1 + kRK1 + kA, RK binds if kG1kG1 + kRK1 + kA 
r < kG1 + kRK1/kG1 + kRK1 + kA, otherwise Arr binds. The lag-time in the free R* state is determined by generating a second random number from an exponential distribution with mean equal to 
kG1 + kRK1 + kA.
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. 
). Noise (both the continuous component of photoreceptor noise and instrument noise) was then added to the photocurrent response, as was done in Fig. 9 of Whitlock and Lamb (1999). Finally, a 2.5 Hz digital low-pass filter was applied to the responses (Whitlock and Lamb, 1999).
Simulation of Other Conditions
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.
Simulations of Transgenic and Genetic Knockout Data
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 kRK1[RK] was set to zero across all n in Eq. 1a. Arr-/- (Xu et al., 1997) was simulated by setting kA[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 kRK3 = 0 across all n.
Altered Levels of ATP and GTP
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 kRK3[ATP] and kG5[GTP] shown in Table I by 0.04 and 0.4, respectively, the same factors used in Rieke and Baylor (1998a) (see Eqs. 
and 
, respectively).
Early Saturation Model
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 kA(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 kP1[PDE] (Eq. ) by the factor (1 - PDE*/PDEmax), where PDE* is the running count of PDE-subunit activations and PDEmax is the total number of PDE subunits locally available. In simulations, we found that in order to achieve empirical levels of SPR variability (CVarea) using only this local saturation mechanism, the number of PDEs locally available to R* (PDEmax) had to be restricted to 300. In contrast, for the full sequential phosphorylation model, PDEmax is set to infinity, so that there is no local saturation (i.e., the response to a single R* will not cause any local depletion of PDE).
Calcium Clamping
Ca2+ clamp was simulated by replacing the time-varying calcium variable, c, with the constant, cdark (so that 
), which reduces the back-end deterministic equations (Eqs. ) to the single equation
 | (13) |
.
Data Analysis
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).
|
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.
Classification of Responses as Failures, SPRs, or Multiple-photon Responses (MPRs)
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).
Isolation of SPR Variance
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.
SPR Area
In the first new analysis, we were interested in CVarea 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. CVarea was calculated as the square root of the variance difference between SPRs and failure areas, divided by the mean SPR area.
Squared Mean SPR Versus SPR Variance
Another new analysis was a comparison of the time-varying ensemble variance (
SPR2) and squared mean (µSPR2), 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 CVampl, and a comparison of the variance of the full dim-flash ensemble (dim2) with the squared mean (µdim2).
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONAPPENDIX A: APPENDIX B: APPENDIX C: REFERENCES |
|---|
. 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.
Variability of SPR Amplitude
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 CVampl 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).
Dim-flash Variance Versus Dim-flash Squared Mean
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 dim2(red) and the µdim2 (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 dim2 and µdim2.
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 dim2 and µdim2 does imply nonstereotypic SPRs.
Variability of SPR Area
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.
SPR Variance Versus SPR Squared Mean
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), SPR2 is approximately an order of magnitude smaller than µSPR2 until after the peak of the squared-mean response. In addition, SPR2 peaks much later than µSPR2 (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.
Transduction Gain Manipulation by Alteration of Nucleotide Levels
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.
Genetic Knockout (KO) and Transgenic Manipulations
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).4 The GCAPs-/- responses rose with WT kinetics to peak at 4 times the WT amplitude at 300 ms.
The Sequential Phosphorylation Model Accounts for All the Data Examined
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).
|
Low SPR Variability Matches Empirical Values
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 
(blue) and the scaled dim2(red) in Fig. 4 C indicates that SPR and MPR variability is low enough to allow the ensemble variance to be dominated by variability in the number of photon absorptions.
The distribution of SPR areas is Gaussian like (Fig. 4 D, red overlay) with a low CV of 0.38 that nearly matches the CVarea 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 CVarea; 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. 
, DISCUSSION), and reflects the contribution of approximately seven of the eight possible inactivation steps. The reason that all eight possible steps did not contribute is that, in our model, the R–Arr affinity increases monotonically with the number of phosphorylations and there is a finite probability of Arr capping as early as the first phosphorylation. Thus, Arr capping occurs before the maximum number of phosphorylations. On average, R* was capped when 6.1 of the 7 possible phosphorylations had occurred (Fig. 4 D, inset), corresponding to a total of 7 shutoff steps (including Arr-capping), which is in agreement with the observed CVarea (
=0.38). The fact that CVarea achieves the theoretical limit for the mean number of phosphorylations at shutoff indicates that the relative rates of G* activation and R* phosphorylation across phosphorylation states were nearly optimal.
The sequential phosphorylation model also reproduces the empirical relationship between µSPR2and SPR2 (compare panel E in Figs. 4 and 3). The SPR2 waveform peaks at 1.6 times the time-to-peak of µSPR2, close to the peak shift when this analysis is applied to the Whitlock and Lamb data (1.5 times; Fig. 3 E).
Rieke and Baylor's Transduction Gain Manipulation Is Reproduced
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.
Qualitative Features of Transgenic and Genetic Knockout Data Are Reproduced
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.
Analytical Expression for R* Activity and the Time Course of R* Inactivation
Because G-protein competes with RK and Arr, and because of reversibility in some of these reactions (Eqs. 
and and ), G* activation rate (i.e., R* activity) is a complicated function of 10 front-end parameters. Thus, it is not feasible to identify a single parameter that controls the maximum G* activation rate or its dependence on n, the phosphorylation state of R*. However, we have solved the stochastic front-end equations to obtain an analytical expression for the mean steady-state activity of R* conditioned on the continued circulation of R* around the catalytic loop. This can be thought of as the rate of G* production at a late time (after the initial appearance of R* in a given phosphorylation state), averaged over all cases where neither capping nor further phosphorylation has yet occurred. The expression thus derived gives a theoretical rate of G* activation per R* (i.e., 
RG, in s-1) that is not physically observable, and therefore is not exactly the same as the chemical reaction rate of G* production. However, because R* typically circulates around the catalytic loop many times before being phosphorylated or capped, we can expect that the true maximum rate of G* production achieved will be close to this theoretical value. Although we do not show the full expression here, it can be approximated by the expression shown in Eq. 14.
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
 | (15) |
