Quantitative modeling of the molecular steps underlying shut-off of rhodopsin activity in rod phototransduction

Sequential phosphorylation of equivalent residues, followed by arrestin binding

As in previous studies, we assume that all the phosphorylation sites, whether Ser or Thr, are equivalent to each other. In the general case illustrated in Figure 1A, the R* molecule passes sequentially through a series of states in which one extra phosphate is added at each step (green arrows). Because the sites are assumed equivalent, it is only necessary to monitor the number of phosphates attached, rather than keep track of which sites they have been attached to. Arrestin can bind to any of these states (blue arrows), thereby terminating R* activity.

The active states, R*·P(n), are shown in red, and the inactivated states, R*·P(n)·Arr, are shown in black. The integer n indicates the number of phosphates attached and runs from zero up to the number of residues that can be phosphorylated, that is, n = 0 ... N, where typically N = 6, as shown explicitly in Figure 1B. The initial activate state, R*·P(0), represents non-phosphorylated metarhodopsin II, which is formed rapidly after photoisomerization. Each state R*·P(n) may undergo one of two reactions: it may be phosphorylated to become R*·P(n+1), or it may bind arrestin to become R*·P(n)·Arr. We denote the rate constant at which R*·P(n) undergoes phosphorylation as ν(n), in green, and the rate constant at which it undergoes arrestin binding as μ(n), in blue. The mean rate at which the molecule exits state R*·P(n) is ν(n)+μ(n), while the probability that the exit reaction corresponds to phosphorylation rather than to arrestin binding is ν(n)/{ν(n)+μ(n)}. Because each of the arrestin-bound states R*·P(n)·Arr has zero ability to activate transducin, these states can all be ignored, and there is no need to indicate any reactions involving them.

Note that we have not yet needed to specify the catalytic activity of R*·P(n). That activity only needs to be known when we solve the downstream deterministic equations.

Analysis of arrestin binding times

We now examine the time to arrestin binding, t_Arr, and its probability distribution, p(t_Arr). For the arrestin binding reaction, we follow recent studies by assuming that some number, M, of phosphates need to have been attached before arrestin can bind and that subsequent states bind arrestin with a fixed rate constant. This assumption is based on the biochemical experiments of Vishnivetskiy et al. [16]. The rate constant μ(n) of arrestin binding to state R*·P(n) is

μ ( n ) = 0 ,  for  n = 0   ...   M - 1 ; μ ( n ) = μ ,  for  n = M   ...   N .            ( 0.1 )

For example, with M = 3, as shown in Figure 1B, the states with zero, one, or two phosphates cannot bind arrestin, whereas each of the subsequent states with three or more phosphates binds arrestin with the fixed rate constant μ. At this stage, we leave the values for the phosphorylation rate constants completely general and simply denote them as ν(0), ν(1) ... ν(N−1). Note that ν(N) = 0 because the fully phosphorylated state cannot be further phosphorylated. The parameters N and M and the required rate constants are described in Table 1, together with their default values.

Amazingly, the assumption embodied in Equation (0.1) leads to a major simplification. Because each of the states with M or more phosphates inactivates by binding arrestin at the common rate, μ, it is irrelevant whether any additional phosphates are added once M have been attached. Therefore, we reach the astonishing prediction that the cascade can be truncated at any stage after state M, and yet the predicted kinetics of arrestin binding will be unaltered. In particular, the cascade can be truncated right at stage M, as illustrated for the case of M = 3 in Figure 1C. What this means is that the only phosphorylation rate constants that affect the kinetics of arrestin binding are those for the first M steps leading up to the formation of R*·P(M). This prediction came as such a surprise that we checked it numerically, by running simulations using different numbers N of available sites for phosphorylation, and for each N ≥ M that we tested, we obtained probability distributions (not shown) that were identical, within the noise.

For reaction schemes of the type illustrated in Figure 1C, the macroscopic solution for the impulse response was analyzed by Baylor, Hodgkin and Lamb [17]. For several simplifying cases of the rate constants, they derived the time course of state M, and in our formulation this represents the probability density function for the time to arrestin binding, t_Arr. From their equations we can readily obtain analytical solutions for the predicted probability density p(t_Arr) in three straightforward cases. The expressions that we derive have utility in two distinct ways. First, these equations provide us with multiple checks on the accuracy of the stochastic simulations that we perform. Second, in our “binary R* shut-off” model, they provide conveniently simple analytical expressions for the distribution of R* lifetimes, that can be used in conjunction with the deterministic set of reactions for the downstream phototransduction cascade to provide analytical expressions for the predicted distribution of SPR amplitudes.

Case 1: All rate constants equal

The simplest case occurs when all the non-zero rate constants are equal, in other words, when all the phosphorylation rate constants (at least those up to M−1) are the same as the common arrestin binding rate constant, μ, defined above in Equation (0.1), so that

v ( n ) = μ ,  for  n = 0   ...   M - 1.     (0 .2)

In this case of “all rate constants equal”, the probability density function for the time to arrestin binding, p(t_Arr), is obtained from Equation (44) in [17] as

p ( t Arr ) = μ ( μ t Arr ) M e − μ t Arr M ! .        (0 .3)

Note that, as we discovered above, this solution is independent of the total number of sites, N, and instead is a function of the minimum number of sites, M, required for arrestin binding. Note also that Equation (0.3) is the same “Poisson equation” that was derived for equal rate constants by Fuortes and Hodgkin [18] in their classic study of Limulus photoreceptor responses.

In this Poisson case, it is straightforward to integrate Equation (0.3) to obtain the cumulative probability density for arrestin binding, P(t_Arr), as

P ( t Arr ) = ∫ 0 t Arr p ( t ) d t = 1 − e − μ t ∑ k = 0 M ( μ t ) k k ! .    (0 .4)

One can also obtain the mean time to arrestin binding and the coefficient of variation (cv = SD/mean) of the time to arrestin binding as

mean ( t Arr ) = ( M + 1 ) / μ      (0 .5) cv ( t Arr ) = 1 / M + 1     (0 .6)

Case 2: Common arrestin binding rate constants and common phosphorylation rate constants

When each of the phosphorylation rate constants (at least up to M−1) is common, ν, but different from the common non-zero arrestin binding rate constant, μ, we have

v ( n ) = v ,  for  n = 0   ...   M - 1     (0 .7)

with the arrestin binding rate constants still defined by Equation (0.1). This case was termed “Final rate constant different” in [17], where the solution was given as their Equation (45). That solution may be rewritten here as

p ( t Arr ) = μ ( ν ν − μ ) M { e − μ   t Arr − e − ν t Arr ∑ k = 0 M − 1 ( ( ν − μ ) t Arr ) k k ! }   (0 .8)

In this case, the mean and cv are more complicated and will not be presented.

Case 3: Rate constants arranged arithmetically

A third tractable simplification arises when the rate constants of phosphorylation and arrestin binding happen to be arranged arithmetically; that is, when

v ( n ) = ( M + 1 − n ) μ ,  for  n = 0   ...   M − 1   (0 .9)

with the common non-zero arrestin binding rate constant μ still defined by Equation (0.1). In this model, with M = 3, the four non-zero rate constants would be 4 μ, 3 μ, 2 μ, and μ. This situation was termed “Independent activation” in [17], where the solution was given as their Equation (37), and may be written as

p ( t Arr ) = ( M + 1 ) μ e − μ   t Arr ( 1 − e − μ   t Arr ) M .     ( 0.10 )

In this case, one can apply integration to obtain the mean time to arrestin binding as

mean ( t Arr ) = 1 μ ∑ k = 1 M + 1 1 k    (0 .11)

which, for M = 3, reduces to mean(t_Arr) = 25/12 μ^-1. It is noteworthy that both Equation (0.5) and Equation (0.11) conform to the expectation for a stochastic cascade that the mean time taken to reach a given state is equal to the sum of the mean dwell times for the individual transitions leading up to that state.

Next, we describe three models for the way in which the R* molecule becomes inactivated as it progresses through the reaction steps illustrated in Figure 1 and analyzed above.

Model 1: Conventional case of graded reductions in R* activity

In the conventional model of graded reductions in R* catalytic activity, we follow the formulation of previous studies by specifying the transition rate constants for state R*·P(n) as

v ( n ) = v max exp ( − ω P n ) ,  for  n = 0   ...   N - 1 ;    v ( N ) = 0.     (1 .1)

This equation indicates that the strength of the phosphorylation reactions declines monotonically as the number of attached phosphates increases, as reported in biochemical experiments [8]. We are not aware of any way of obtaining an analytical solution for the probability distributions of the states in this case. Together, Equations (0.1) and (1.1) contain three parameters (ν_max, ω_P, and μ) that need to be specified before stochastic simulation can be undertaken.

For subsequent substitution into the downstream transduction cascade, we need the catalytic activity of each state, R*·P(n). It is convenient to express this in fractional terms, ρ(n), where ρ(0) = 1 for the initial state R*·P(0). As in previous studies, this fractional catalytic activity is presumed to decline exponentially with n:

ρ ( n ) = exp ( − ω G n ) ,  for  n = 0   ...   N .    (1 .2)

Here, we allow for the possibility that the steepness parameter, ω_G, may be different from the steepness parameter ω_P for phosphorylation in Equation (1.1); this was also permitted in [9], whereas in [6] and [11] it was assumed that ω_G = ω_P.

Model 2: “Binary” R* activity

In the second model, we dispense with the assumption that the catalytic activity of R*·P(n) declines with the increasing number of attached phosphates, and instead we assume that the activity ρ(n) is constant. We refer to this as “binary” R* activity, because the R* remains fully active until it binds arrestin.

ρ ( n ) = 1 ,  for  n = 0   ...   N .    (2 .1)

To determine the time course of R* activity in this model, we only need the time t_Arr until arrestin binds, and it was the realization of this fact that provided the main impetus for deriving analytical expressions for the probability distributions of t_Arr in the preceding section.

Of particular importance is the cumulative probability distribution, P(t_Arr), because in the binary activity model the activity of R* is unity until arrestin binds, whereupon it drops to zero. Therefore, the mean R* activity is simply the complement of P(t_Arr), that is,

R * bin ( t ) ¯ = 1 − P ( t ) .    (2 .2)

Therefore, for the case of “all rate constants equal”, Equation (0.4) gives the mean time course of R* activity as

R * bin ( t ) ¯ = e − μ t ∑ k = 0 M ( μ t ) k k ! .    (2 .3)

At very early times, this expression approximates to

R * bin ( t ) ¯ t → 0 ≈ 1 − ( μ t ) M + 1 ( M + 1 ) ! .  (2 .4)

Therefore, in the binary R* model, the initial time course of R* activity does not approximate an exponential decay but is instead “flat-topped”, as shown in the Results section.

The expected form for the probability distribution of SPR amplitudes can be readily determined in the binary R* model, because the solution of the downstream reaction cascade involves only a single stochastic parameter, t_Arr, for which we know the probability distribution in certain simplifying cases. Therefore, if we solve the downstream cascade of reactions for a range of values of t_Arr, we can obtain any desired feature of the SPR as a function of t_Arr; for example, we can obtain the peak amplitude, SPR_peak, as a function SPR_peak(t_Arr) of t_Arr. Then, for any case where we know the probability density p(t_Arr), we can evaluate the probability density of SPR_peak as

p ( S P R peak ) = p ( t Arr ) / ( d S P R peak / d t Arr )    (2 .5)

In practice, the function SPR_peak(t_Arr) can be evaluated at a small number of values of t_Arr and then fit using a spline function and differentiated. This approach is illustrated subsequently in the Results section.

Model 3: Three-state model with low R* activity before arrestin binding

One drawback of the binary R* model turns out to be its inability to predict a sensible SPR waveform in the arrestin knockout phenotype (see Results). This occurs because in that case the R* activity never shuts off. To deal with this issue, we have investigated an alternative model in which the R* activity drops to a lower level after M phosphates have been attached. As illustrated in Figure 1D, we envisage that once M phosphates have attached, the R* molecule can enter a new state according to a rate constant κ. We further envisage that, following this change, arrestin can now bind, again with rate constant µ (Figure 1D). Therefore, instead of Equation (0.1), we have

κ ( n ) = 0 ,  for  n = 0   ...   M - 1 ; κ ( n ) = κ ,  for  n = M   ...   N .    (3 .1)

Therefore, simply by replacing μ with κ, the analytical expressions that we derived above for arrestin binding time now apply instead to the time t_low to reach the conformation with low activity. We additionally need to specify the arrestin binding rate constant as

μ = 0 ,  initially, then  μ = μ low  in the low activity state .   (3 .2)

Finally, we need to specify the R* activity, which we assume drops to a low value, ρ_low, and then drops to zero upon arrestin binding, which we formulate as

ρ = 1 ,  initially, then  ρ = ρ low ,  then  ρ = 0 ,  when Arr binds .   (3 .3)

To derive the mean R* activity deterministically in this three-state case, it is convenient to define the following three variables: P_F(t), the probability that the molecule remains in the fully activated state; P_FL(t), the probability that it is either in the fully activated or in the low-activity state (i.e., the probability that it is not arrestin-bound); and the difference P_L(t) = P_FL(t) − P_F(t), the probability that it is in the low-activity state. In the three-state case, P_F(t) is analogous to R*_bin(t), except with μ replaced by κ, while P_FL(t) is a corresponding expression with the additional rate constant μ. Accordingly, we can express the mean R* activity in the three-state case as P_F(t) plus a scaled-down version of P_L(t) to give

R * tri ( t ) ¯ = ( 1 − ρ low ) P F ( t ) + ρ low P FL ( t ) .    (3 .4)

In the general case, with arbitrary rate constants, it is necessary to evaluate the terms P_F(t) and P_FL(t) in Equation (3.4) by numerical integration. However, for the simplest case, with all rate constants equal (ν = μ = κ), Equation (3.4) reduces to

R * tri ( t ) ¯ = e − μ t [ ∑ k = 0 M ( μ t ) k k ! + ρ low ( μ t ) M +1 ( M + 1 ) ! ] .    (3 .5)

As in the case of the binary R* activity model, we can also predict the form of the distribution of SPR peak amplitudes, although in this case there is an approximation involved. Thus, when we obtain the relationship between SPR peak amplitude and the time t_low to the low-activity state, our approximation is to take the subsequent time to arrestin binding as fixed at 1/μ rather than being stochastic. Our subsequent results show that this provides a remarkably accurate prediction. By analogy with Equation (2.5), the predicted distribution in the three-state model is

p ( S P R peak ) = p ( t low ) / ( d S P R peak / d t low ) .    (3 .6)

Downstream phototransduction cascade

For any of the models above, we can simulate the time course of R* activity, or, in the special cases of simple ratios of rate constants, we know the probability density of arrestin binding times. We can then solve the deterministic reactions for the downstream phototransduction cascade to generate the predicted individual SPRs. To do this, we employed a well-established model of phototransduction, with Ca²⁺-mediated feedback onto cyclic GMP concentration via GCAPs and guanylyl cyclase. We used a distributed model with longitudinal diffusion of cGMP and Ca²⁺ within the outer segment, and as a check we also used a reduced set of these equations for the isotropic case of a well-stirred outer segment.

Photoisomerizations are deemed to occur at discrete spatial locations, that is, at disks, and PDE activity is likewise deemed to be restricted to those disks so that its activity is specified by an ordinary differential equation (o.d.e.). For both cGMP and Ca²⁺, a partial differential equation (p.d.e.) is required to describe the longitudinal diffusion. The three differential equations and their ancillary equations are as follows. Most variables include the spatial (x) coordinate as well as time. The parameters in these equations are defined in Table 2.

d E * ( x , t ) / d t = v RE R * ( x , t ) − k E E * ( x , t )    (4 .1)

∂ c G ( x , t ) ∂ t = φ cG ( x , t ) + D cG ∂ 2 c G ( x , t ) ∂ x 2    (4 .2)

∂ C a ( x , t ) ∂ t = φ Ca ( x , t ) + D Ca ∂ 2 C a ( x , t ) ∂ x 2    (4 .3)

φ c G ( x , t ) = α ( x , t ) − β Dark c G ( x , t )    (4 .4)

φ Ca ( x , t ) = 1 2 f Ca L   j cG ( x , t ) − L   j ex ( x , t ) ℱ B Ca V cyto    (4 .5)

α ( x , t ) = α max 1 + ( C a ( x , t ) / K GCAP ) m GCAP    (4 .6)

L   j cG ( x , t ) = J cG , max c G ( x , t ) 3 c G ( x , t ) 3 + K cG 3    (4 .7)

L   j ex ( x , t ) = J ex,max C a ( x , t ) C a ( x , t ) + K ex    (4 .8)

J ( t ) = ∫ L { j cG ( x , t ) + j ex ( x , t ) } d x    (4 .9)

The membrane currents above are expressed as current density per unit length of outer segment, j_cG(x,t) and j_ex(x,t). However, to simplify comparison with the conventional lumped formulation, they have been written in the form L j_cG(x,t), with units of current over the entire outer segment length, L. Consequently, when these equations are coded numerically, we can conveniently use the same symbols in both the spatial and lumped cases and then use either the integral Equation (4.9) in the spatial case or a simple sum in the lumped case.

Boundary conditions at ends

The boundary conditions at each end of the outer segment are of the reflective (zero flux) kind. They apply for both cGMP and Ca²⁺, and can be written as:

∂ c G ( x , t ) ∂ x = 0 = ∂ C a ( x , t ) ∂ x , at each end of the OS .   (4 .10)

Boundary condition at location of photoisomerization

For these partial differential equations, Equation (4.4) above does not incorporate a term for light-induced PDE activity because the spatial extent of R* activity at each location where a photoisomerization has occurred is deemed to be infinitesimally narrow. Therefore, at each such location we need to apply a boundary condition for cGMP hydrolysis. In the symmetric case of a single photoisomerization delivered at the middle of the outer segment, defined as x = x₀, this is straightforward. Equating the flux in the positive x direction with the rate at which cGMP is hydrolyzed by half of the activated PDE on that disk, we obtain the boundary condition as

∂ c G ( x , t ) ∂ x | x 0 + = E * ( x 0 , t ) β sub L 2 D cG c G ( x 0 , t ) .    (4 .11)

The full derivation of this equation is given in Appendix 1, and it is identical to Equation (B 1) derived by Lamb and Pugh [1]. However, we note that the corresponding Equation (6) in Gross et al. [11] differs in having a factor of 4 in its denominator, but in our view that formulation was incorrect. Gross et al. [11] invoked a second factor of 2, on account of the fraction of the outer segment volume that is cytoplasmic. However, analysis of the equations shows that this is not appropriate. One way of rationalizing this is to note that the volume factor has already been accounted for in the definition of the effective longitudinal diffusion coefficient. In comparing the equations between the two studies, one also needs to note that they used a different parameter for the PDE hydrolytic efficacy that they termed β_idv, which is N_disks times our β_sub, that is, β_idv = N_disks β_sub, where N_disks is the number of disks in the outer segment.

Although this analytical form of the boundary condition is important in estimating the hydrolytic activity underlying the rogue responses, it is not required for numerical simulation. For numerical simulation, we instead used spatial elements of finite width δx, and we took the PDE activity to be distributed over the width of the element where the photoisomerization occurred. Therefore, in any element in which a photoisomerization occured, Equation (4.4) was replaced with the following modified form of the conventional lumped case equation:

φ cG ( x , t ) = α ( x , t ) − { β Dark + β sub L δ x E * ( x , t ) } c G ( x , t ) .    (4 .12)

Ensemble variance for dim-flash responses

One can consider an idealized experiment, where the number of photoisomerizations per trial is Poisson-distributed with a mean of m, and where the individual SPRs have a mean time-course a(t) and variance time-course σ₁²(t). In this idealized case, the expected ensemble variance for a very long series of identical flashes is

σ 2 ( t ) = m { a 2 ( t ) + σ 1 2 ( t ) } .    (5 .1)

Typically, one would expect the ratio σ₁(t)/a(t) to be quite small (around 0.3; see Results section) so that its square would be around 0.1. As a result, the experimentally recorded variance is dominated by the Poisson contribution of the first term in Equation (5.1), and the “between SPR” variability elevates the predicted variance by only about 10%. In practice, the situation is more complicated because of the finite number of flashes that can be delivered, nonlinear summation (especially when the light stimulus is spatially restricted), noise, and lack of recording stability.

Numerical simulation of stochastic reactions for R*

The stochastic reactions for R*, described in the panels in Figure 1, were simulated using Gillespie’s method [19]. We will begin by describing the simplest cases (Figure 1A–C) where just two fates are possible for each state: either phosphorylation or arrestin binding. First, the mean transition rate constants μ(n) and ν(n) were assigned. Then, for each states, two pseudo-random numbers, r₁(n) and r₂(n), uniformly distributed between zero and unity, were used. The lifetime that the molecule remained in state R*·P(n) was simulated as −ln(r₁(n))/{ν(n)+μ(n)}, and the state to which the molecule transitioned was simulated as phosphorylation if r₂(n) < ν(n)/{ν(n)+μ(n)} and as arrestin binding otherwise. If a finite flash width was required, the time of photoisomerization was simulated as uniformly distributed over that flash width. These simulations were very fast, and 10⁶ iterations could be performed in ~1 s on a laptop PC. For the three-state R* activity model, simulation of the full version (Figure 1D) would have required more complicated code, so we chose instead to simulate only the truncated linear configuration (Figure 1E), and in this case 10⁶ iterations again took only ~1 s.

Numerical simulation of stochastic reactions for E*

For most of the simulations described here, the time course of SPR-induced PDE activity, E*(t), was taken to be deterministic, as described by the o.d.e. Equation (4.1). In one case, however, we additionally simulated the generation and decay of E* as a stochastic variable. At the single location of the photoisomerization, E* molecules were generated stochastically at a mean rate given by ν_RE times the instantaneous R* activity (ρ = 1 or ρ = ρ_low), and each generated E* was assigned a stochastic lifetime with a mean of k_E⁻¹.

Numerical integration of equations for downstream phototransduction cascade

The differential equations for the deterministic reactions of the downstream phototransduction cascade were integrated numerically using a Matlab o.d.e. solver, typically “ode15s”. For the case of longitudinal diffusion in the outer segment, the method of partial discretization was employed [20] to represent the system as a set of spatially coupled o.d.e.s, with the outer segment divided into a finite number of spatial elements (50 in the illustrated Figures).

The parameters we chose are specified in Table 2. The values in the “Default” column were used for each model, with the exception that for Model 1 the three values in the column “Model 1” were used. We began by evaluating the dark resting state as follows. For Equations (4.1)–(4.9), we set all the spatial and temporal derivatives to zero and expressed the steady cGMP in terms of the steady Ca using Equation (4.4) with (4.6)–(4.8). We then determined the value of Ca that rendered the net influx of calcium zero by numerically finding the zero for the numerator of Equation (4.5). The resulting dark resting values are listed in Table 2. Note that the resting state was identical in each model because we constrained the ratio α_max/β_Dark to be the same in each case.

The numerical integration of the downstream reactions was orders of magnitude slower than the stochastic simulations, and 10⁵ repetitions of the spatial diffusion case typically took around 5 h on a laptop PC. In practice, as few as 1000 repetitions were sufficient to provide extremely good estimates of the ensemble mean and variance of the SPR. However, to obtain a good histogram of the distribution of SPR peak amplitudes for the illustrated Figures we ran 10⁵ repetitions.

Model 1: Conventional model of R* graded shut-off

We began by verifying that our implementation of the stochastic model of graded R* shut-off in conjunction with our formulation of the downstream reactions predicted SPRs with properties conforming to those recently reported by Gross et al. [11]. To this end, we adopted parameter values that were nearly identical to theirs, as listed in Table 1. One minor difference was that for consistency with our subsequent simulations of Models 2 and 3, we chose the “abrupt” arrestin-binding cut-in described by Equation (0.1), as used recently in [10], rather than the graded approximation to this relation employed by Gross et al. [11] in their Equation (10).

Figure 2 collects the results of our simulations for the graded R* shut-off model. The top pair of panels summarize the behavior of individual R* molecules in 10⁶ repetitions. The remaining panels show SPRs predicted by integrating the downstream cascade reactions: the middle pair of panels each show a handful of SPRs, while the bottom pair plot the ensemble behavior for 10⁵ simulations.

Figure 2A compares the probability distribution for R* integration times, p(τ_R*), that we obtained (blue) with the corresponding curve from [11] (red) that was plotted as the inset in their illustration 6E. The similarity is very close, although our simulations gave marginally fewer short integration times and a marginally greater number of long integration times. Figure 2B compares the time-course of the mean of the simulated R* activity for the 10⁶ repetitions and shows that the traces were barely distinguishable between the two studies. Figure 2C plots the first 50 SPRs simulated and gives a flavor of the individual waveforms; Figure 2D selects some outlier responses and will be considered shortly. Figure 2E plots the ensemble mean and SD for 10⁵ SPRs integrated using the downstream reactions, and the traces are extremely similar to the corresponding curves illustrated in [11]. Figure 2F plots the probability distribution for the peak amplitudes of these 10⁵ SPRs (in blue), and the close resemblance to the distribution obtained in [11] (red circles) is striking.

These simulations generated the following values of interest, which are compared with the corresponding values from [11] in brackets, where available. For the dark steady-state: cG_Dark, 4.1 µM (4.1 µM); Ca_Dark, 322 nM (320 nM); and total outer segment current, 18.4 pA (18.7 pA). For the stochastic responses: R* integration time, mean 39.7 ms (40 ms), and cv 0.547 (0.52); cv of SPR area, 0.459; fractional amplitude of SPR at peak, mean 0.0420, and cv 0.332; and mean time-to-peak of SPR: 119 ms. Based on the closeness of these values and the similarity of the four comparisons in panels A, B, E, and F, we conclude that our SPR simulations conformed very closely to those of Gross et al. [11].

Despite this close similarity, we observed a disturbing feature of a minority of the simulated SPRs that does not accord with our experience of the form of real SPRs. As exemplified by Figure 2C,D, a minority of the simulated responses tended to exhibit a plateau after the initial peak. In Figure 2C, it is apparent that a sub-set of the first 50 simulated responses exhibited a fairly flat region at around half the amplitude of the peak; qualitatively the same phenomenon was apparent in any group of 50 consecutive simulations that we selected (not shown). Investigation of this phenomenon revealed that the plateau resulted from trials that happened to exhibit relatively long durations in the doubly phosphorylated state, R*·P(2). Therefore, to better illustrate these events, we selected outlier responses representing 5% of the total trials based on two criteria, and we present these in Figure 2D. From the first 200 trials we selected the 10 responses that had the longest times to exit from state R*·P(2), and these are plotted in blue. Each of these traces has the appearance of heading toward a fractional amplitude of ~0.042 at times of 400–500 ms. Next, from the same set of 200 trials we selected the 10 responses that exhibited the longest time-to-peak, and these are plotted in red. Most of these red traces reach peak later than 300 ms at a level of roughly 0.037. We observed qualitatively the same kind of behavior when we selected the outliers from any group of 200 consecutive trials (not shown).

We suggest that this tendency for a sizable fraction of trials to exhibit a plateau is an inevitable feature of the graded shut-off model and arises (at least for the parameter values chosen) from trials that happened to exhibit a long dwell time in state R*·P(2). Although each state is broadly analogous to every other state, there are some factors that combine to generate “unusual” events for the doubly phosphorylated state. First, the mean lifetime in this state is moderately long, at 1/(80 s⁻¹ e⁻²) = 92 ms, so that 5% of trials are expected to display dwell times in this state of longer than ~270 ms. Second, the amplitude elicited by a long dwell time in this state happens to asymptote to a level close to the peak amplitude of the mean SPR (see dotted horizontal lines on the right in Figure 2D). As a result, the inevitable events that exhibit long dwell times in state R*·P(2) generate responses that are noticeably unusual. For the other states, the effects of comparatively long dwell times are visually less noticeable, either because the mean lifetime is quite short (states n = 0 or 1) or because the asymptotic response level for that state is considerably smaller than the peak (states n ≥ 3).

We think it very likely that it is the existence of a moderate number of SPRs exhibiting such a plateau that causes the rather “flat-topped” shape for the time-course of the ensemble standard deviation, shown by the red trace in Figure 2E. In addition, we think that it is the existence of responses such as those shown in Figure 2D that gives rise to the sharp peak in the amplitude histogram (between 0.035 and 0.042) in Figure 2F. We are not aware of published experiments that have reported SPRs exhibiting kinetics of this kind, and nor have we observed responses that obviously resemble these events in our own recordings of dim flash responses. Nevertheless, we cannot completely rule out the possibility that SPRs with this shape do occur in real rods. However, our firm inclination is to suggest that this feature predicted by the model of graded reduction in R* activity is inconsistent with the properties of observed rod responses.

In addition to this issue with plateau-like events, another major problem was that the required amplification of the downstream reactions was excessively high, being greater than can be accounted for by the measured biochemical parameters. Thus, the product ν_RE β_sub that was needed both in Model 1 (for ν_RE = 300 s⁻¹ and β_sub = 0.063 s⁻¹) and in [11] was ν_RE β_sub ≈ 19 s⁻². This corresponds to an amplification constant for phototransduction of A = ν_RE β_sub n_cG ≈ 57 s⁻², far higher than the values reported in experiments on mammalian rods (see Discussion).

For comparison with biochemical estimates, β_sub is given by (½ k_CAT/K_m) / (½ πr²L BP N_Av) [1], where k_CAT/K_m is the catalytic activity of the fully activated PDE dimer, BP is the buffering power of the cytoplasm for cGMP, N_Av is Avogadro’s number, and r = 0.75 µm and L = 22 µm are the radius and length of the outer segment. The factor of 1/2 in the numerator reflects the fact that two G* molecules are required to fully activate the dimeric PDE, while the factor of 1/2 in the denominator corresponds to the fraction of the envelope volume that is cytoplasmic. For amphibian rods at room temperature, k_CAT/K_m has been measured as 4.4 × 10⁸ M⁻¹ s⁻¹ [21]. Substitution of these values gives β_sub = 0.019 s⁻¹ / BP. If the cytoplasm exhibited finite buffering power for cGMP (i.e., BP > 1), then the value would be lower than 0.019 s⁻¹. Conversely, if the mammalian enzyme exhibited higher activity at body temperature, then the value could be higher (see Discussion). However, taking these two factors into account, it is difficult to envisage that the biochemical measurements could correspond to a hydrolytic rate constant of more than ~50% above 0.019 s⁻¹, which would set a realistic constraint that β_sub ≤ 0.03 s⁻¹. On this basis, the value of β_sub = 0.063 s⁻¹ required in Model 1 again appears excessive.

An independent estimate of β_sub for mouse rods was recently reported by Gross et al. [13], based on measurements of “rogue” responses representing the failure of activated R* to shut off. However, in our view their estimate needs to be halved for the following reason. As derived in Appendix 1 and noted following Equation (4.11), their analysis incorrectly introduced a second factor of 2 in the denominator of their Equation (6). Accordingly, their analysis of rogue responses actually predicts a value of β_sub = 0.03 s⁻¹ in GCAPs^−/− rods. Again, this is approximately half the value of β_sub required to achieve responses of appropriate amplitude in Model 1.

An additional problem with the graded R* shut-off model is that its prediction of the arrestin-knockout phenotype is not very realistic. When the rate constant for arrestin binding was set to zero (μ = 0), the predicted mean SPR exhibited a tail considerably smaller than observed in experiments, decaying to <4% of the peak after 5 s (not shown). This occurred because in the graded shut-off model the activity of state R*·P(6) is extremely small, at e⁻⁶ ≈ 0.0025, with the consequence that the activity of R* is substantially quenched even in the absence of arrestin binding.

The combination of these problems leads us to conclude that the graded R* shut-off model, in the form described both here and in [11], does not provide a plausible description of the observed SPRs in mammalian rods.

Model 2: “Binary” R* activity

In light of these difficulties with the model of a graded decline in R* activity, we chose to investigate a model in which the R* activity does not decline in a graded manner with successive phosphorylations and instead remains maximal until arrestin binds. We refer to this situation as “binary” R* activity. Such a model has the convenience that the probability distribution of R* lifetimes can be readily predicted in several simplifying cases, as derived in the Theory section. Not only does this provide a means of checking the accuracy of our stochastic simulations, but it also allows us to predict the distribution of SPR amplitudes.

The results of our simulations for the binary R* shut-off model are collected in Figure 3, where the organization of panels is broadly similar to that in Figure 2. The top pair of panels summarize the simulated activity of R* molecules in 10⁶ repetitions; the middle pair show a selection of SPRs; and the bottom pair plot the ensemble behavior for 10⁵ simulations.

For the stochastic shut-off parameters, we chose to set all the rate constants equal, both for simplicity and also so that we could examine whether the results were described by the theory developed in Equations (0.3)–(0.6) and (2.3)–(2.5); we used ν = μ = 60 s⁻¹ (Table 1), and we set the flash duration to zero. For the downstream reactions, we used the default parameters in Table 2. These differ in the following respects from those used in the graded R* shut-off model in [11] and duplicated in Figure 2. First, we reduced the maximum cyclase rate to the value used by Koch and Dell’Orco [22], α_max = 120 µM s⁻¹. Next, to keep the dark resting state unchanged, we held the ratio α_max/β_Dark constant by reducing β_Dark in the same ratio to β_Dark = 3.2 s⁻¹. Finally, we adjusted the magnitude of the PDE hydrolytic activity to generate SPRs of the expected amplitude, ~ 4% of the dark current. The resulting value of β_sub = 0.024 s⁻¹ (Table 2) turned out to be very close to that predicted both by biochemical measurements and by the analysis of “rogue” R* responses [13].

For 10⁶ simulations of the R* stochastic reactions, Figure 3A plots the probability distribution for arrestin binding times p(t_Arr) as the blue histogram. For comparison, the smooth red curve is the prediction of Equation (0.3), and the fit is extremely good. For M = 3 and ν = μ = 60 s⁻¹, its peak is predicted to occur at 3/(60 s⁻¹) = 50 ms, as observed. The mean and the cv of the arrestin binding time (and thus of R* activity) are predicted by Equations (0.5) and (0.6) to be 66.67 ms and 0.5, respectively, and the values we obtained in the simulations were 66.66 ms and 0.4995. In Figure 3B, the blue curve plots the ensemble mean time course of R*(t) for the 10⁶ binary R* activity events, while the dashed red trace plots the prediction of Equation (2.3) and is seen to completely overlie the blue curve for the simulations. The closeness of the simulated activity to the predictions gives us confidence that our simulation algorithm provided an accurate rendition of the stochastic reactions.

Figure 3C superimposes the first 50 SPRs simulated, and two features are apparent. First, in this case with zero flash duration, the traces begin rising along a common curve. Second, there is no hint of any plateau-like responses of the kind seen in Figure 2. Figure 3D plots the SPRs obtained for a series of specified values of t_Arr, ranging from 5 ms to 800 ms, and confirms the common early rising phase. The red circles mark the peaks of these individual traces, with filled red used for t_Arr up to 130 ms and open symbols used for t_Arr equal to or greater than 200 ms. The magnitudes of these peaks are subsequently used to generate the relationship between SPR_peak and t_Arr that is required in Equation (2.5) for predicting the distribution of amplitudes. Note that the probability of t_Arr > 150 ms was 2.2%.

For 10⁵ simulated SPRs, Figure 3E plots the time-course of the ensemble mean and SD and shows a characteristic feature of experimentally measured SPRs, which is that the ensemble SD reaches peak almost 40 ms later than the ensemble mean at 163 ms compared to 126 ms. Figure 3F plots the probability distribution for the peak amplitudes of the 10⁵ SPRs as the blue histogram and compares this with the prediction of Equation (2.5) in red, demonstrating that the prediction is very accurate. The peak amplitudes of the individual SPRs (measured as a fraction of the dark current) had a mean of 0.0413, with a cv of 0.348, considerably reduced from the cv for arrestin binding time of 1/2.

Our interpretation of these results is as follows. First, the good fit of the predictions to the simulated results supports the view that our algorithms accurately simulated the stochastic reactions. Second, the binary R* activity model with the parameters we have chosen generates simulated responses that account for many of the features observed for mammalian rod SPR experiments. Third, the binary R* activity model does not exhibit the first two problems we noted for the graded R* activity model; that is, it does not generate occasional SPRs that exhibit a plateau-like late phase and it only requires a value for β_sub that is consistent with predictions from biochemical measurements. We suggest that adoption of a realistic estimate for the gain product ν_RE β_sub (of approximately 7 s⁻²) necessitates the use of an integration time for R* activity of the order of 65 ms or more to generate SPRs of realistic amplitude. In our view, the proposal of an R* integration time of 40 ms in previous work seems unrealistic and appears to have arisen as a consequence of assuming that R* activity shuts off in a graded manner.

Despite the good description of WT SPRs obtained in Figure 3, there is a problem with the binary R* model when the arrestin knockout phenotype is modeled (by setting μ = 0). As the R* activity then never shuts off, the simulated response simply rises to a maximum and remains there. This behavior is not illustrated but is essentially the same as shown subsequently for the GRK1 knockout (obtained by setting ν = 0). We therefore investigated an extension of the binary R* model in which the catalytic activity declined to a low level after a certain number of phosphates had been attached.

Model 3: Three-state model, with low R* activity before arrestin binding

To deal with the failure to account for the arrestin knockout phenotype, we modified the binary R* activity model to include a penultimate state with lower activity, as shown schematically in Figure 1D. In this “three-state” model, we hypothesize that the R* catalytic activity remains maximal until a change is triggered by the attachment of some threshold number of phosphates. When this change occurs, the catalytic activity drops to a low level, and a new conformation amenable to arrestin binding is formed. As indicated in Figure 1D, we assume for simplicity that the rate constant of entering the low-activity conformation is fixed (at κ) for all states with M or more phosphates so that the scheme simplifies to the linear sequence illustrated in Figure 1E. It is straightforward to simulate this linear scheme stochastically and likewise to obtain analytical solutions in the simplifying cases of ratios of rate constants, as set out in the Theory section.

The results of our simulations for this three-state model are collected in Figure 4, where the organization of panels follows that in Figure 3. The top pair of panels summarize the simulated activity of the R* molecule in 10⁶ repetitions; the middle pair show a selection of responses; and the bottom pair plot the ensemble behavior for SPRs in 10⁵ simulations.

For the downstream reactions, we retained the parameters employed for the binary R* model. For the stochastic shut-off parameters, we retained M = 3, and we set the rate constant κ for the change to the low-activity conformation equal to the other rate constants at κ = ν = μ = 60 s⁻¹. To account for the late phase of the arrestin knockout phenotype, we chose the low-state activity as ρ_low = 0.1; see Equation (3.3).

For 10⁶ simulations of the R* stochastic reactions, Figure 4A plots the probability distributions for both the time to entering the low-activity state (t_low, blue histogram) and the time to arrestin binding (t_Arr, black histogram). For comparison, the smooth red curves plot the predictions of the Poisson expression, Equation (0.3), and in both cases the fit is excellent. In Figure 4B, the dashed blue curve plots the mean time course of R*(t) for the ensemble of 10⁶ R* events in this three-state case. The overlying dashed red curve shows the prediction of Equation (3.5), while the dotted red curve shows the prediction of the binary R* activity model, Equation (2.3). The mean of the R* activity (simulated or predicted) barely differs from the prediction of the binary activity model, indicating that the inclusion of the third state, with low activity ρ_low = 0.1, causes little change in the mean R* activity in this WT case.

Figure 4C plots the first 50 simulated SPRs, and the traces are qualitatively similar to the corresponding simulations in Figure 3C, with no hint of outliers exhibiting unusual kinetics. Figure 4D superimposes the SPRs obtained for a series of specified values of t_low, ranging from 5 ms to 800 ms, with a fixed subsequent delay to arrestin binding; see the Figure caption. The red circles mark the peaks of these individual traces, with filled symbols used for t_low from 5 ms to 130 ms and open symbols used for t_low from 200 ms to 800 ms. The magnitudes of these peaks were subsequently used to generate the relationship between SPR_peak and t_low that is required for predicting the probability distribution of amplitudes. Note that the probability of t_low > 150 ms was 2.2%, and the probability of t_Arr > 150 ms was 5.6%.

For 10⁵ simulated SPRs, Figure 4E plots the time course of the ensemble mean and SD. Just as was the case for the binary model, the ensemble SD reaches peak around 40 ms later than the ensemble mean, at 167 ms compared to 127 ms. Figure 4F plots the probability distribution for the peak amplitudes of the 10⁵ SPRs as the blue histogram and compares this with the prediction in red, obtained using Equation (3.6) with the measurements from Figure 4D. Again, the prediction is very accurate. With all three rate constants set to 60 s⁻¹, the SPR peak amplitudes had a mean of 0.0423 and a cv of 0.339. We therefore conclude that the three-state model generates responses that are barely distinguishable from those of the binary activity model for these parameters applicable to WT mammalian rods.

Stochastic generation and decay of E*

For the simulations above, we have approximated the time course of PDE activity, E*(t), as being deterministic, according to Equation (4.1) and as used in recent modeling [11]. However, because of the modest number of E* molecules activated during the SPR, this approach has its limitations. For example, for a mean R* lifetime of τ_R* ≈ 70 ms and a mean rate of E* activation of ν_RE ≈ 300 s⁻¹, the mean number of total E* subunits activated throughout the duration of the SPR would be n_E* = τ_R* ν_RE ≈ 21, and the number active at any instant would be smaller than this. We therefore now follow Hamer et al. [6] and Reingruber and Holcman [9], who considered stochastic activation and shut-off of E* in their analyses of the graded R* shut-off model (i.e., Model 1), and we simulate the stochastic E* case applied to the three-state Model 3.

Figure 5 shows the results obtained when stochastic activation and shut-off of E* were added to the stochastic R* activity for Model 3. We used the same seed for the random number generator and thereby ensured that the sequence of R*(t) trials was identical to that obtained in Figure 4. Panel A illustrates the simulated E*(t) time course for the first 50 SPR trials, corresponding to the same R* events that were used to generate Figure 4C above, and shows substantial fluctuations in the number of activated E* subunits. Panel B plots the mean (blue) and SD (red) of E*(t) for the entire set of 10⁵ trials. In addition, the dashed black trace plots the deterministic time course (see Equation (4.1)) predicted by convolving the mean R*(t) time course with an exponential decay, exp(−k_E t), and it provides an excellent fit. Panel C plots the photocurrent responses for the first 50 SPR trials and therefore corresponds to Figure 4C but with the addition of stochastic generation and decay of E*. Note that at late time in the simulations (>400 ms), discrete levels of response are visible, corresponding to the presence of the last few E* subunits in each trial. Panel D plots the mean (blue) and SD (red) of the SPRs for the entire set of simulations. For the parameter values listed in Table 2 (and in particular for this value of β_Dark), the waveforms in panels C and D are quite similar to those in panels A and B, indicating that in this case the downstream reactions introduce relatively minor filtering of the E*(t) kinetics. The ensemble of SPRs had a mean peak amplitude of 0.0427 and a cv of 0.382. Panel E plots the relationship between peak SPR amplitude and the time to the low-activity state of R* for the sample traces shown previously in Figure 4D.

Figure 5F plots the probability distribution for the peak amplitudes of the 10⁵ SPRs as the blue histogram. Note that this trace exhibits multiple minor peaks that we interpret as corresponding to discrete numbers of E* molecules being present at around the time of the peak; the first eight such peaks are visually apparent with this bin width. With fewer than 5 × 10⁴ trials or with broader sampling bins, multiple peaks like this were not readily resolvable (not shown). To predict an expected distribution in this more complicated case with stochastic E* time-course, we began with the same approach as in Figure 4 and plotted the predicted histogram from Figure 4F as the dotted red trace in Figure 5F; this curve had its peak at 0.0415. Next, we assumed that stochastic fluctuations in E* numbers would contribute variability with a coefficient of variation 1/√n_E*. Therefore, we convolved the dotted histogram with a Gaussian distribution having an SD of 0.0415/√21 to generate the continuous red trace in Figure 5F. Although we cannot provide a rigorous justification for this prediction, we note that it provides a reasonable description of the histogram obtained by simulation, apart from the absence of the multiple minor peaks. We conclude that the inclusion of stochastic steps for the generation and decay of E* leads on the one hand to SPR responses (Figure 5C) that more nearly resemble the shape of observed dim flash responses but on the other hand leads to only minor changes in the form of the ensemble mean, SD, and amplitude distribution.

Simulation of responses for mutant phenotypes

To investigate the wider applicability of our model, we repeated our simulations of the three-state R* activity case (Model 3) after adjusting the parameters to correspond to those expected for four mutations in the phototransduction cascade: Arr^+/−, Arr^−/−, GRK1^+/−, and GRK1^−/−. For the arrestin hemizygote, we halved μ and for the arrestin knockout, we set μ to zero. For the rhodopsin kinase hemizygote, we halved ν and for the rhodopsin kinase knockout, we set ν to zero. In the case of the three-state model, Figure 6 presents the ensemble mean responses for the WT case and these four mutants. As a check, we employed a deterministic description for R*_tri(t), corresponding to Equations (3.4)–(3.5), as driving functions for numerical integration of the downstream reactions, and the computation time was then extremely rapid. This generated traces that were very similar (not shown), although the simulated responses were marginally smaller at around the peak of the mean SPR, presumably because of nonlinearity caused by greater response compression for larger SPRs than for smaller ones.

The ensemble mean response for the Arr^+/− case (dashed red trace) differs only marginally from the WT response (black), whereas the Arr^−/− response exhibits a final steady level of roughly half the peak amplitude of the WT SPR. The magnitude of this final level in the Arr^−/− simulation depended on the choice of catalytic activity in the “low-activity” R* state, and a magnitude of roughly half the WT peak amplitude was obtained by setting ρ_low = 0.1. The GRK1^+/− mean response is about 30% larger than the WT response, whereas the GRK1^−/− response rises to a plateau at a little over double the peak amplitude of the WT SPR. In qualitative terms, each of these simulated mean responses bears a remarkable similarity to the respective SPR responses reported in experiments on mouse rods; see, for example illustration 8 in [23] with data from [24,25] and illustration 1 in [26], with data from Dr. M.E. Burns. Responses very similar to those recorded from GRK1^−/− rods are also obtained for rods in which rhodopsin has a truncated C-terminus [27,28].

Inadequacy of graded R* shut-off model

In our view, the conventional model, of a sequence of graded steps in the reduction of R* catalytic activity as recently simulated in [11] and duplicated here in Model 1, fails to provide an adequate description of the SPRs. Although the ensemble mean and SD resemble those reported in the literature, there are serious problems with two other aspects of the predictions.

First, the graded R* shut-off model is characterized by a small proportion (5%–10%) of simulated SPRs that exhibit either a plateau-like region or a very late peak (Figure 2D). This result appears to be an inevitable consequence of the assumption of an exponential decline in both the catalytic activity and the rate constant of phosphorylation. These occasional events do not correspond to published reports of mammalian rod SPRs. Their occurrence appears to generate the rather sharp peak in the amplitude distribution histogram of the simulated responses (Figure 2F), a feature that again has not been reported in published experiments.

Second, the required “gain parameter” needs to be set to roughly double the value predicted from other experiments. This gain parameter can be represented either by the product ν_RE β_sub or equivalently by the “amplification constant” A = ν_RE β_sub n_cG [1], where n_cG is the cooperativity of channel opening by cGMP, which is widely accepted to be n_cG = 3. To obtain a peak amplitude for the WT SPRs of ~4% of the dark current, the gain product needed to be set to ν_RE β_sub ≈ 19 s⁻² both in the present work and in [11], corresponding to an amplification constant of A ≈ 57 s⁻². We now compare the magnitude of this gain factor required for the graded R* shut-off simulations with three estimates derived from experiments, two using electrophysiological approaches and one using biochemical measurements.

Electrophysiological measurements of the amplification constant can be obtained by fitting the rising phase of responses to families of flash responses [1], and reported values in the literature range from A = 5–23 s⁻² in experiments on rodent rods [29-33]. It is likely that some of the lower values resulted from excessive low-pass filtering [32], and our best guess for the true value of A in WT mouse rods is A ≈ 24 s⁻², corresponding to ν_RE β_sub ≈ 8 s⁻². Recently, a second completely independent estimate for this gain parameter was obtained by Gross et al. [13] in their experiments on “rogue” SPRs in mouse GCAPs^−/− rods, where they obtained a value that we have corrected to ν_RE β_sub ≈ 9 s⁻² (see Results section for Model 1). Biochemical experiments on amphibian rods [21] can be extrapolated (with caution) from room temperature to mammalian body temperature. For the rate of E* activation per R*, we anticipate an upper limit of ν_RE ≤ 350 s⁻¹, but there are some uncertainties in estimating β_sub. Although it is often thought that the Q₁₀ for biochemical reactions may be around 2, the relevant Q₁₀ for the PDE activity is likely to be smaller than this because the reaction rate is already quite near the limit for aqueous diffusion, and the increase in temperature of ~15 °C (from 22 °C to 37 °C) is unlikely to cause as much as a twofold increase in rate. Accordingly, we think that the amphibian value of k_CAT/K_m ≈ 4.4 × 10⁸ M⁻¹ s⁻¹ [21] will convert to a mammalian value of <8 × 10⁸ M⁻¹ s⁻¹. Substitution into the expression for β_sub (see Appendix 1, Equation A.10) then gives β_sub <0.033 s⁻¹ / BP, where BP ≥ 1 is the cytoplasmic buffering power for cGMP. As a result we think that the biochemistry points to a realistic upper limit for the gain parameter of mouse rods of ν_RE β_sub ≈ 12 s⁻², provided that BP ≈ 1, but it would be lower than this if the cytoplasm exhibits any buffering power for cGMP. In summary, the three experimental estimates for ν_RE β_sub are 8, 9, and <12 s⁻², and these are all considerably smaller than the value of 19 s⁻² that is required in the graded R* shut-off model to fit the observed SPR amplitude.

An alternative possibility that we tested, to try to avoid an excessive magnitude for the gain parameter in the graded shut-off model, was to slow the R* shut-off reactions from the values proposed by Gross et al. [11,13] to increase the mean R* lifetime. When we halved the parameters ν_max and μ (thereby doubling the mean R* lifetime), and reduced ν_RE β_sub to 12 s⁻², the simulated SPRs had approximately the correct mean amplitude, but several other features were unusual and did not resemble physiological responses. The shape of the mean response had an unrealistic, considerably more “triangular” shape, with its width at half-height increased to ~400 ms, the tail phase was too slow, the time-to-peak was slightly increased to 132 ms, and the outlier events became much more prominent, appearing to be double the duration of those illustrated in Figure 2D. Accordingly, we rejected the possibility that the use of slower R* shut-off reactions could avoid the need for an excessive gain parameter in the graded R* shut-off model.

Binary and three-state R* shut-off models

Because of these deficiencies in the graded R* shut-off model, we investigated two variants of a scenario in which R* retains maximal activity until at least three phosphates have attached (referred to here as M = 3). In the first variant, the “binary R* shut-off” model, the R* activity remains at unity until it drops abruptly to zero when arrestin binds. In the second variant, the “three-state R* shut-off” model, we invoked an additional state in which the R* catalytic activity dropped to a low level, ρ_low, before the binding of arrestin. These two variants of abrupt shut-off predicted nearly identical behavior to each other and one that conformed closely to the behavior reported in the literature (Figure 3 and Figure 4) in the case of parameters appropriate for WT mouse rods. However, the binary model was not able to account for the kinetics of the responses recorded from the rods of Arr^−/− mice, and we required the additional low activity state. As shown in Figure 6, the simulated responses of the three-state model provided a satisfactory description of the ensemble mean SPRs simulated for the rods of four genotypes of mouse: Arr^+/−, Arr^−/−, GRK1^+/−, and GRK1^−/−. The height of the late plateau level in the Arr^−/− simulations provided our basis for choosing the low-state activity as ρ_low = 0.1.

Inclusion of stochastic reactions for E* generation and decay

The sample simulated SPR responses for both versions of the abrupt R* shut-off model, illustrated in Figure 3C and Figure 4C, exhibit a smooth time course for each response because all the steps other than R* formation and decay have been assumed to be deterministic. To approximate more closely to reality, we added stochastic simulation of the generation and decay of E* to the downstream reactions to generate the responses illustrated in Figure 5. With the inclusion of E* fluctuations, the sample SPR responses (Figure 5C) more nearly resembled the shape of observed dim flash responses. Incorporation of this stochastic description caused only minor change in the form of the ensemble mean, SD, and amplitude distribution of the SPRs.

R* time-course and integration time

Our modeling and simulations provide some new perspectives on the time course of R* activity and on the R* “integration time”, τ_R*, defined as the area under the time course of fractional R* activity, averaged across the simulated trials. In particular, we find that to achieve a sufficiently large SPR amplitude (~4% of the dark current) with a realistic magnitude of the gain product ν_RE β_sub ≈ 7 s⁻², the R* integration time τ_R* needs to be of the order of 60–70 ms, considerably larger than the value of 40 ms recently reported [11,34] for WT mouse rods. We think that our estimate applies irrespective of the shape of the decline in the ensemble mean R* activity, that is, that the estimate τ_R* ≈ 60–70 ms is dependent primarily on the downstream model and in particular on its “gain” parameters rather than on the stochastic model chosen to describe R* activity.

Ensemble mean R* time course

In the binary and three-state models, the time course of the ensemble mean R* activity does not approximate an exponential decay but instead exhibits an initial flat top as illustrated in Figure 3B and Figure 4B, with the early time course given by 1 – k t^M+1; see Equation (2.4). To the extent that our models are realistic in this regard, this predicted time course provides some noteworthy insights.

First, there is an issue with the extraction of estimates of R* lifetime from “Pepperberg plot” experiments conducted on mutant mice with different expected rates of R* shut-off. In such experiments, the shift of the relationship between saturation time and log intensity has been used to provide estimates of R* lifetime [34]. However, the theoretical basis for that approach is predicated upon the assumption that R*’s mean activity declines exponentially; see [34]. Accordingly, we urge caution in the interpretation of these estimates of R* lifetime, unless one obtains clear evidence to support the assumption of an exponential decline in R* activity.

Another issue relates to the accuracy (or otherwise) of the “delayed Gaussian” prediction of the Lamb and Pugh [1] model for the onset phase of the rod’s response to bright flashes of light. That model ignored shut-off reactions and was therefore explicitly restricted to times shorter than the fastest shut-off reaction. Nevertheless, experiments on electroretinogram (ERG) scotopic a-wave responses in both human subjects and anesthetized mammals demonstrate that the predictions are quite accurate until at least 20–30 ms; see, for example illustrations 4C, D in [35]. If it were the case that the mean R* activity in the rods declined exponentially with a time constant of 40 ms, then the activity would have declined to ~60% after 20 ms, and one might expect the recorded photocurrents to deviate substantially from the Gaussian approximation. In contrast, if the mean R* activity in rods follows a time course similar to that illustrated in Figure 3B or Figure 4B and remains above 98% at 20 ms, then one might indeed expect the recorded photocurrents to conform reasonably closely to the Gaussian approximation until at least 20 ms.

These concepts are examined in Figure 7, where families of bright-flash responses have been predicted under five different assumptions concerning the shut-off reactions. When neither R* nor G*/PDE* shuts off, the predictions are shown in black, with the dashed traces giving the Gaussian approximation and the continuous traces giving the solution of the differential equations, with the R* shut-off rate constant k_R* = 0 and with k_E = 0. The green traces present the solution when G*/PDE* shuts off normally (k_E specified in Table 2) but R* does not shut off (k_R* = 0). Next, the red traces plot the predictions for the three-state model with parameters as in the Tables and in Figure 4. The blue traces plot the predictions for the case of exponential shut-off of R* with a time constant of k_R*⁻¹ = 40 ms.

These calculations show that the Gaussian approximation fails to provide a perfect fit to the “no shut-off” case. They also show very close agreement between the green and red traces, for both of which G*/PDE* shuts of normally and where the only difference is whether R* fails to shut off (green) or whether it shuts off according to the three-state model (red). Importantly, the blue traces for an R* exponential lifetime of 40 ms are seen to deviate markedly from the other traces beyond about 20 ms after the flash. In other words, the observed good fit of the Lamb and Pugh theory to the ensemble of responses is more difficult to account for if the R* activity is assumed to drop exponentially with a time constant as short as 40 ms than is the case where the R* activity is assumed to be “flat-topped”.

Finally, as far as the R* time course is concerned, the previously assumed rapid exponential decline in R* activity has always been difficult to rationalize from a teleological perspective. Given that it took hundreds of millions of years for rod phototransduction to evolve the exquisite sensitivity needed for single-photon detection, it seems difficult to envisage what advantage there could be in causing that very high R* activity to begin declining (on average) from the instant that the R* is formed. Instead, one might envisage that it would be advantageous to maintain that hard-won catalytic activity for as long as possible and then to shut it off rapidly. This, in fact, is approximately what the binary and three-state models achieve.

Physical nature of the postulated low-activity state

To account for SPRs in the arrestin-knockout phenotype, we invoked a low-activity state that the R* molecule enters once several phosphates have been attached, but we did not specify the nature of this state. One possibility would be that once several phosphates have bound, the R* molecule undergoes a conformation change, with this altered conformation exhibiting not only greatly increased affinity for arrestin but also lowered activity toward transducin. Another possibility would be that once several phosphates have bound, the kinase (GRK1) binds more tightly than previously, making access for transducin more difficult. However, we emphasize that our simulations provide no prospect for distinguishing the precise nature of such a postulated low-activity state.

Comparison with biochemical estimates of R* affinity for transducin and for arrestin

The biochemical measurements of transducin affinity (and hence of assumed R* activity) reported by Gibson et al. [8] were obtained using macroscopic experiments in which the mean number of bound phosphates across different preparations was measured. Any interpretation about the probability distribution of rhodopsin molecules having specific numbers of attached phosphates is model-dependent, involving, for example, the rate constants of phosphorylation under the particular experimental conditions. In particular, one cannot assume that a graded relationship between measured affinity and mean number of bound phosphates in any way reflects a graded relationship between R* activity and actual numbers of bound phosphates. Future experiments will be needed to measure the catalytic activity of the molecule in its different phosphorylation states. As far as we can see, our postulated abrupt drop in R* activity is in no way inconsistent with the observed graded decline in biochemical affinity measured as a function of the mean number of phosphates. For the binding of arrestin, Vishnivetskiy et al. [16] measured the abundance of rhodopsins bearing different numbers of phosphates in their different samples, and they then measured the degree of arrestin binding to these samples. They found that high-affinity binding of arrestin required at least two bound phosphates and was maximal with three. Despite that finding, recent stochastic models of R* activity (including ours) have assumed that the binding of just two phosphates is insufficient to cause rapid binding of arrestin and that instead three are required.

Plateau-like SPRs

Plateau-like SPR events have been routinely observed experimentally when phosphorylation of R* is disrupted, as, for example, in GRK1^−/− rods [25], in rhodopsin C-terminus mutants [27,28] and in certain cases of substitution of the Ser and Thr residues in the carboxy terminus [36,37]. In addition, they are also observed in WT rods in a very small proportion (~0.3%) of flash presentations, and as such they have been termed “aberrant” or “rogue” responses [13,38]. Typically, these aberrant SPRs exhibit a plateau level around 50% higher than the normal peak, followed by an abrupt recovery to baseline. To explain these experimentally observed aberrant responses in the conventional model of graded reduction in R* activity, it has been necessary to assume that phosphorylation fails to occur for a small minority of R* events. However, in terms of our three-state stochastic model, an alternative explanation would be that phosphorylation occurs normally but that the subsequent transition to the low-activity state is delayed for the order of seconds. It is difficult to speculate on a possible molecular mechanism, but it might involve either a change to or a degradation of a tiny fraction of the rhodopsin molecules or alternatively a disruption of some component with which R* interacts.

In conclusion, we think that our simulations rule out the conventional graded shut-off model as providing a realistic description of experimentally recorded SPRs. And while we emphasize that our simulations do not provide any direct evidence that the activity of R*·P(n) remains unattenuated up to n = 3, we think they demonstrate that a model of this kind is capable of accounting for many of the properties of experimentally recorded SPRs.

Our model opens up some exciting possibilities for the future. The analysis of additional experimental data may permit refinement of the model and provide improved estimates of the underlying physical parameters. In particular, we think that such future refinements may readily account for the different classes of SPRs recorded by Azevedo et al. [37] using mice with Ser-only and Thr-only C-terminus mutations of rhodopsin. Furthermore, now that we have a deterministic expression (Equations (3.4)–(3.5)) for the mean R*(t) time course in the three-state shut-off model, it should be straightforward to extend the model to account for more complicated aspects of light adaptation, including modulation by recoverin, modulation of channel activity, and so on.