Copyright © 2008 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 94, Issue 6, 1954-1970, 15 March 2008
doi:10.1529/biophysj.107.116202
Biophysical Theory and Modeling
The Dynamics of Phosphodiesterase Activation in Rods and Cones
Jürgen Reingruber* and David Holcman*, †, 
, 
* Département de Biologie, Ecole Normale Supérieure, Paris, France
† Department of Mathematics, Weizmann Institute of Science, Rehovot, Israel
Address reprint requests to David Holcman.
Abstract
Phototransduction starts with the activation of a rhodopsin (respectively, coneopsin) molecule, located in the outer segment of rod (respectively, cone) photoreceptors. The subsequent amplification pathway proceeds via the G-protein transducin to the activation of phosphodiesterase (PDE), a G-protein coupled effector enzyme. In this article, we study the dynamics of PDE activation by constructing a Markov model that is based on the underlying chemical reactions including multiple rhodopsin phosphorylations. We derive explicit equations for the mean and the variance of activated PDE. Our analysis reveals that a low rhodopsin lifetime variance is neither necessary nor sufficient to achieve reliable PDE activation. The numerical simulations show that during the rising phase the variability of PDE activation is much lower compared to the recovery phase, and this property depends crucially on the transducin activation rates. Furthermore, we find that the dynamics of the activation process greatly differs depending on whether rhodopsin or PDE deactivation limits the recovery of the photoresponse. Finally, our simulations for cones show that only very few PDEs are activated by an excited photopigment, which might explain why in S-cones no single photon response can be observed.
Introduction
Phototransduction is a multistep process which starts when a photon activates a rhodopsin (respectively, coneopsin) molecule in the outer segment of a rod (respectively, cone) photoreceptor. Upon diffusional encounter on internal disks in rods and on the surface membrane in cones, the activated opsin binds successively to many copies of transducin, a G-protein coupled receptor. Finally, each of the activated transducin binds to a single phosphodiesterase (PDE) effector protein 1,2,3,4,5,6. The set of activated PDE molecules hydrolyze cGMP, a cytosolic diffusible second messenger, which leads to the closure of cGMP-gated ion channels and thus to the photoreceptor hyperpolarization. In rods, physiological studies have revealed that even the absorption of a single photon can be detected 2,7,8,9, while for cones, many quasisynchronous absorbed photons (approximately seven) are needed to generate a signal that overcomes the noisy background 10,11,12. Remarkably, in rods also the single photon response time course is very reproducible (see, e.g., 2). Despite of major progresses, it is still a challenging problem to unravel the precise mechanisms responsible for the accuracy and reproducibility of the single photon response in rods.
A high reliability of the rod single photon response implies a low variability of the number of activated PDE. This condition can be achieved by controlling accurately the amplification process. Several factors are involved in this amplification, such as the lifetime of activated rhodopsin and the rates of transducin activation. The chemical reactions that control the deactivation of rhodopsin depend on rhodopsin kinase, recoverin, and arrestin 4,5,6. Recent studies 13,14,15,16,17 have suggested that the reproducibility of the single photon response might be due to a low variability in the lifetime of activated rhodopsin, achieved through rhodopsin deactivation via multiple phosphorylation steps. However, it is still unclear how many deactivation steps are necessary to reproduce the experimental data. In Field and Rieke 14 it was suggested that at least 12–14 steps are needed; however, numerical studies based on Monte Carlo simulations have shown that already seven phosphorylation sites are sufficient to reproduce experimental data 18,19.
To extract the main principles underlying the variability of the photoresponse, we present here a stochastic analysis of PDE activation for both rod and cone photoreceptors. Our model is based on the well-accepted molecular cascade leading to the activation of the G-protein. To analyze the PDE dynamics and the associated fluctuation, we derive equations for the mean and the variance of activated PDE. Since our approach allows us to compute the time course of the mean and the variance of excited PDE, it complements previous stochastic simulations 18,20,21. We derive analytic expressions for the mean and the variance of rhodopsin lifetime and the number of activated PDE, and provide numerical simulations. Furthermore, we study the influence of various parameters such as the number of rhodopsin phosphorylation sites and phosphorylation and transducin activation rates. We explore the impact of the rhodopsin lifetime on the accuracy of PDE activation. We study PDE response for scenarios representing rods and cones in mice and toads. We show that during the rising phase, the PDE variability is much lower compared to the recovery phase. We also analyze the role of whether rhodopsin or PDE lifetimes limit the recovery of the photoresponse. Our results show that the variability of the PDE response depends most significantly on the transducin activation rates. Finally, we present simulations suitable for cones. We find that in cones only very few PDE molecules are activated, which confirms an earlier suggestion 22.
Theory
Model for PDE activation
The transduction process following the absorption of a photon has been well documented both experimentally and theoretically (for reviews see 1,2,3,4,5,6). The Pugh-Lamb model 1,20,23 was based on the properties of two-dimensional random collisions and it predicts accurately the rising phase of the photoresponse. Based on a Markovian approach and using stochastic simulations of diffusion and chemical reactions, in Felber et al. 21, the mean and the variance of the simulated photoresponse were obtained for different lifetimes of activated rhodopsin. We approximate here the different steps leading to PDE activation by first-order chemical reactions and we neglect the molecular dynamics due to diffusion. This approximation is justified by the large number of molecules and the fast diffusion constant. Thus, the number of activated molecules resulting from diffusional collisions has the same temporal law compared to reaction equations.
The well-accepted scenario for PDE activation embodied in our model is the following (see Table 1 for a guide to the parameters used): After a photon absorption, a rhodopsin molecule, denoted by R, undergoes a conformational modification and changes from an inactive into an active form R*. R* deactivation occurs through multiple phosphorylation steps, catalyzed by rhodopsin kinase (RK), and finally through arrestin binding 15,24,25,26,27. We take into account that the affinity of R* for transducin, RK, and arrestin are altered by sequential phosphorylations 18,28. The number of rhodopsin phosphorylation sites is denoted by Np. After photon absorption, rhodopsin changes into the activated state n=N. The parameter N equals the total number of R* deactivation steps, and we assume that it is given by N=Np+1 13 (number of phosphorylation sites+arrestin binding). When R* encounters RK, with a certain probability a phosphorylation occurs and R* undergoes a transition from the state n to n-1, modeled by a state-dependent phosphorylation rate λn. When R* binds to arrestin, there is a certain probability that R* becomes deactivated, modeled by a transition rate μn from the state n to the deactivated state n=0. In each state n>0, R* activates G-proteins transducin (T) with a rate kact(n). While T* can bind to a PDE with a rate k3 to form a complex denoted by PDE*, the same complex can be deactivated with a rate k4. The reciprocal of the rate k4 is the lifetime of PDE* and depends crucially on the concentration of RGS9 29,30,31,32. We neglect depletion of transducin and PDE because the amount of activated molecules is negligibly small compared to the total pool of available transducin and PDE molecules. The kinetic reactions underlying the model (illustrated in Fig. 1) are summarized as
The model presented in Eq.
(1) is a stripped-down version of a more detailed model, such as the one presented in Hamer et al.
18. To keep the model simple and clear, we did not include backward reaction rates. We decided to focus on conditions that can lead to a minimal PDE* variance and therefore omitted backward reaction rates, which would certainly increase the variability of the activation process, as was already noticed in Field and Rieke
14. However, our analysis can be extended without much effort to include backward reaction rates and intermediate bound states between R* and RK. For example, the reactions used to model R* phosphorylation in Eqs.
(1) in Hamer et al.
18 can be incorporated by using appropriate backward rates
(
n now labels also intermediate states),
| | |
| Symbol | Description | |
|---|
| R*n | Activated rhodopsin in state n | |
| R0 | Deactivated rhodopsin | |
| T | Transducin | |
| T* | Activated transducin | |
| PDE | Phosphodiesterase | |
| PDE* | Activated PDE | |
| Np | Number of rhodopsin phosphorylation sites | |
| λn | Phosphorylation rate in state n | |
| μn | Arrestin binding rate in state n | |
| kact(n) | Transducin activation rate in state n | |
| k3 | PDE activation rate | |
| k4 | PDE* deactivation rate | |
| Rτ | Mean/SD ratio of rhodopsin lifetime | |
| RPs | Steady-state mean/SD ratio of PDE* | |
| RP(t) | Time-dependent mean/SD ratio of PDE* | |
|  | Maximum number of PDE* | |
| | |
In contrast, there is no straightforward method to extend our analysis to incorporate also intermediate bound states between R* and transducin (for example, as modeled by the reactions Eqs. (3) in 18). Indeed, as will be seen later on, our mathematical derivations rely on the assumption that R* deactivation occurs independently from transducin and PDE activation. In this case, the analysis of R* deactivation can be decoupled from the analysis of transducin and PDE activation, which greatly reduces the complexity of the computations. Such a decoupling is justified when the lifetime of possible bound states between R* and transducin is short compared to the lifetime of the phosphorylation states. Since transducin activation occurs much faster compared to R* phosphorylation, we assume that under normal conditions R* deactivation occurs almost independently from transducin and PDE activation.
We now proceed with the analysis of the chemical reactions given in Eq. (1). To describe the state of our model, we introduce three stochastic variables (
) that can adopt the values (n, l, k): the phosphorylation state n of R*, the number l of T*, and k of PDE*. The dynamics of the joint probability P(n, l, k, t), that at time t we find R* in the state n and l T* and k PDE* molecules, satisfies a Master equation 33,34. To derive this equation, we first determine from the system of chemical reactions displayed in Eq. (1) the transition matrix W(n, l, k|n′, l′, k′, t) between two states (n, l, k) and (n′, l′, k′),
where
δx,y is the Krönecker delta. From the general shape of the Master equation
33,
34we obtain, by inserting Eq.
(2) into Eq.
(3) (we suppress the initial indices (
n0,
l0,
k0,
t0)),
where
N=
Np+1. The boundary conditions are
The conditions
λ0=
μ0=0 express that the deactivated state is stable. The condition
λ1=0 accounts for the fact that in the state
n=1 all sites are phosphorylated. Immediately after photon absorption, R* is in state
n=
N and the number of
T* and PDE* are zero. Thus, the initial condition for
P(
n,
l,
k,
t) is given by
Dynamics of activated rhodopsin
We now analyze the dynamics of R*. In particular, we estimate the mean and the variance of the duration until R* becomes deactivated by arrestin binding. We start by computing the probability P(n, t) that a R* molecule is in the state n at time t. Our analysis ends with an estimation of the mean to the standard deviation (SD) ratio of R* lifetime.
State probability P(n, t)
To describe the dynamics of R*, we sum Eq. (4) over the indices l and k. We obtain an equation for the probability P(n, t) to find R* in the state n at time t,
where
βn=
λn+
μn. The initial condition is given by
P(
n,
t=0)=
δn,N. Using the vector notation with
we rewrite this system as
where the matrix
S is given by
For pairwise different eigenvalues
βn we can diagonalize
S and using the eigenvectors of
S we derive explicit expressions for
P(
n,
t). With the notation
pi=
λi/
βi,
P(
n,
t) is given by
Mean and variance of rhodopsin lifetime
To compute the mean and the variance of the random R* lifetime T, we use the probability PR(t) that R* is still active at time t, given by
A direct computation using
P(
n,
t) from Eq.
(9) (and the identities in Eq.
(71) and Eq.
(72)) yields for the mean and the variance of the R* lifetime
Equation
(11) for the mean R* lifetime has an intuitive interpretation: it is the sum of the mean lifetimes in each state
n multiplied by the probability to reach this state before being deactivated by arrestin binding.
Reliability of rhodopsin lifetime
We characterize the reliability of R* lifetime by the ratio of the mean to the standard deviation, denoted by Rτ, which is simply the reciprocal of the coefficient of variation (CV). If Rτ is high (respectively, low), the reliability is high (respectively, low). By using Eqs. (11), and following the analysis in the Appendix , we obtain the following estimate:
The upper limit for
Rτ depends only on number of R* deactivation steps
N. Based on very general considerations, this result has already been anticipated
16,
17. However, by using the explicit formulas we can now study the behavior of
Rτ as a function of the underlying rates. Indeed, the maximum value
is achieved if
βn=
const and
pn=1 for all
n. The condition
βn=
const conveys that all deactivation states need to have the same lifetime, while
pn=1 is fulfilled if the arrestin binding rates
μn vanish for all
n>1. Thus, the latter condition states that arrestin only binds when R* is in the state
n=1 and is therefore fully phosphorylated. This assumption is reasonable, since nonvanishing arrestin binding rates for
n>1 effectively reduce the number of deactivation steps and therefore increase the variance.
Mean and variance of the transducin activation rate
The mean and variance of the transducin activation rate are defined as
In principle, the mean and the variance of the activation rate (and any other quantity that depends only on the phosphorylation state of R*) can be computed by using the probabilities
P(
n,
t) given in Eq.
(9). However, we will use an alternative method of calculation that relies on differential equations and the decomposition of the activation rate (see below). This method does not require us to explicitly compute
P(
n,
t) and facilitates the derivation of differential equations for the cross-correlation terms (see Eq.
(34)). We now derive a differential equation for
by differentiating Eq.
(14) with respect to time and by using Eq.
(7),
where
ST is the transposed matrix of
S. To solve this equation, we decompose the vector
into the sum of eigenvectors
of the matrix
ST,
with
By induction, we obtain for the components of the vectors
The differential equation for the mean of
is given by
Using the initial condition that at time
t=0 R* is in the state
n=
N yields
Note that the expression for
provided in Eq.
(19) can be verified by comparing Eq.
(21) with the result for
obtained by inserting Eq.
(9) into Eq.
(14).
The variance 
is calculated analogously to 
by decomposing the vector 
Mean and variance of activated transducin and PDE
The mean and variance of T* and PDE* are defined as
Differentiating Eqs.
(22) with respect to time and using the Master equation Eq.
(4) yields
The definitions for the variance and the correlations read
The equations for the time derivatives of
and
are given by
To close this system of equations we additionally have to derive differential equations for
and
. This is done by using the decomposition of
kact given in Eq.
(17). We first write
The time derivatives of
and
are derived as
The correlations between
kact and
defined by
are computed by decomposing the vectors
into the sum of eigenvectors
of the matrix
ST,
This yields
In practice, the coefficients
are computed numerically by diagonalizing the matrix
S⊤. However, by using Eq.
(19), also analytic expressions can be derived.
Finally, we define the time-dependent PDE reliability ratio RP(t) as
Summary
The explicit expressions for 
and 
allow us to close the system of differential equations for the variance of PDE*. This system consists of Eq. (21) and Eqs. (24) for the mean and Eqs. (31), Eqs. (35), and Eq. (39) for the variances. The simulation results will be obtained by using this close system of equations.
Mean and variance of the total number of activated PDE
It is usually assumed that reliable R* deactivation entails reliable PDE activation 13,14,15,16; however, it is worthwhile to have a closer look at the connection between R* lifetime and PDE activation. For this we compute the mean and the variance of the total number of PDE* molecules produced during an single photon response (SPR), obtained by setting k4=0, and compare it to the variance of R* lifetime. For vanishing PDE* deactivation rate k4, after R* shutoff, a steady state will be reached that contains all the PDE* molecules activated during the SPR. In the Appendix we derive expressions for the steady-state mean and variance of PDE*,
These expressions share similarities with the ones obtained for R* lifetime (Eqs.
(11)), since merely 1/
βn is replaced by
kact(
n)/
βn. We define the steady-state reliability ratio
RPs, which corresponds to the reliability of the total number of PDE*, as the ratio of the steady-state mean to the SD. The reliability ratio
RPs is the inverse of the coefficient of variation (CV). In the
Appendix we obtain a sharp upper bound for
RPs,
In the case when R* activates many PDE* (
), the upper limits for
and
Rτ are both equal to
However, it is interesting to examine whether maximal values for
and
Rτ can be attained simultaneously. The conditions
kact(
n)/
βn=
const and
pn=1 are required such that
achieves its maximum, whereas
βn=
const and
pn=1 are needed to maximize
Rτ. The condition
kact(
n)/
βn=
const expresses that, in each state, the same amounts of PDE* have to be activated, whereas 1/
βn=
const requires that each state has the same lifetime. By adjusting the activation rates
kact(
n), the condition
kact(
n)/
βn=
const can be achieved even when the rates
βn are very different. In that case,
RPs can be maximal while
Rτ is far from being maximal. In general, because of the transducin activation rates, maximal values for
RPs and
Rτ are not achieved simultaneously, which shows that reliable R* lifetime is neither necessary nor sufficient to achieve reliable PDE activation. Only for constant transducin activation rates,
kact(
n)=
kact, the steady-state results for PDE* are determined by the mean and variance of R* lifetime,
We shall now discuss some aspects of the reliability ratio
RPs, which contains the variability of all the molecular events contributing to PDE activation.
RPs is closely related to the coefficient of variation of the area below the PDE* time response, denoted by
CVareaP. In Hamer et al.
18 it was shown that
It is remarkable that
CVareaP depends only on the molecular details involved in PDE activation and not on the PDE* deactivation rate
k4. Consequently,
CVareaP does not depend on whether PDE* or R* deactivation limits the recovery of the SPR. For large numbers of activated PDE, it follows that
The coefficient of variation
CVarea of the area below the SPR current was introduced in Field and Rieke
14. It measures the integrated variability of the SPR: when PDE* dynamics and the SPR current are related by a scaling relationship, we have that
CVarea≈
CVareaP18, which finally leads to
(for sufficiently large
). Consequently, the measured value of
CVarea can be used to obtain a lower bound of the number of R* deactivation steps (by using Eq.
(43)). Indeed, experimental results for mutated mouse rods
13 show that
CVarea behaves approximately like 1/
, the limiting behavior of
(this is the case under the assumption that the number of deactivation steps
N correlates with the number of phosphorylation sites
Np through
N=
Np+1). We conclude that
CVarea is a close measure of
, but in general, it is not for the coefficient of variation of R* lifetime
CVτ.
Results
Numerical simulations of PDE dynamics
To study the PDE dynamics (mean and variance), we run numerical simulations of Eq. (21), Eqs. (24), Eqs. (31), Eqs. (35)and Eq. (39). Our aim is to examine the influence of the various parameters on the dynamics of PDE activation.
Choice of parameters
It has been well documented that the time course of the single photon response (SPR) differs substantially between amphibian and mammalian rod photoreceptors 13,14,16,17,35. In mouse rods, the maximum of the SPR amplitude occurs at ∼0.1s 13,35,36, whereas in toad rods, the maximum of the amplitude occurs at ∼1.9s 16,17. Hence, we decided to run simulations for two different scenarios called Mouse Rod and Toad Rod (see Table 2). Following recent results 35, the mouse rod scenario is characterized by a R* lifetime of 0.080s and a PDE* deactivation rate of 5s−1. Unfortunately, as far as we know, similar experimental data is not available for toad rods. To match the time course of the toad rod photoresponse, we chose for the toad rod scenario a R* lifetime of 3s and a PDE* deactivation rate of 1s−1. Most important, we chose the toad rod parameters such that, contrary to the mouse rod scenario, recovery is limited by R* lifetime and not by PDE* deactivation. This will allow us to explore the impact of whether rhodopsin R* or PDE* lifetime limits the recovery. To compare simulations, we decided to fix the maximum of the mean number of PDE* 
at a value of 150, as suggested in Leskov et al. 37. It is important to note that for our purpose this value is not critical, because other values will result in a simple scaling. Finally, we chose k3=50s−1 for the T*-PDE* binding rate (in 18,19 the authors use k3=200s−1). However, the exact value for k3 is not very important, as long it is not rate-limiting. Our choice of the parameters is summarized in Table 2.
| | |
| | Mouse rod scenario | Toad rod scenario | |
|---|
| PDE activation rate (k3) | 50s−1 | 50s−1 | |
| PDE deactivation rate (k4) | 5s−1 | 1s−1 | |
| Mean rhodopsin lifetime (τ) | 0.080s | 3s | |
| Max. value of  ( ) | 150 | 150 | |
| | |
Motivated by previous studies 18,28, we consider that transducin activation rates kact(n) decay exponentially with the number of rhodopsin phosphorylations, that is
where
ωact is an adjustable parameter. We also assume that the affinity of rhodopsin kinase for R* decays exponentially with the number of phosphorylations,
where
ωλ is also a free parameter.
Since we are interested in conditions leading to the smallest PDE* variance, we will mostly present simulations for a simplified scenario, where the arrestin binding rates μn vanish unless R* is fully phosphorylated (which is the case when R* is in the state n=1). Such a scenario is optimal to achieve a high R* deactivation reliability. Furthermore, we choose the arrestin binding rate for n=1 equal to 
When ωλ=ωact, this choice ensures that kact(n)/βn=const and therefore maximizes the reliability ratio RPs. Moreover, the choice 
adapts the arrestin binding rate to the phosphorylation rates. In summary, the arrestin binding rates will be chosen as
Although we made an effort to decrease the number of free parameters, we still have to specify
N,
λN,
kact(
N),
ωλ, and
ωact. However, the rates
λN and
kact(
N) are fixed by adjusting rhodopsin's lifetime and
. For given values
N and
ωλ,
λN is determined by rhodopsin lifetime according to Eq.
(11). For given values
N,
ωλ, and
ωact we determine numerically the value of
kact(
N) by fixing
. Thus, the remaining parameters that have to be specified are
N,
ωλ, and
ωact.
Impact of the number of phosphorylation sites
Using the theory developed in the previous section, we now study the impact of the number of R* phosphorylation states on the PDE* response. Such an analysis is particularly relevant, since there are transgenic experiments with reduced number of rhodopsin phosphorylation sites 13,15. We run some simulations for the mean and the variance of PDE* for toad rod parameters and ωλ=ωact=0.1. The condition ωλ=ωact ensures that the ratio kact(n)/βn values are constant and thus leads to a maximal steady-state reliability RPs (see Eq. (43)). This condition has also been used previously for simulating the photoresponse 18,19.
Figure 2a shows that during the photoresponse the mean number of T* is very small. Since each T* binds to only one PDE, the number of transducin and PDE molecules that become activated are equal. However, unlike PDE*, T* does not accumulate due to the large rate k3=50s−1. Contrary to the assumption that during the rising phase the ratio of the number of PDE* to T* is constant 1,3,23, we found by comparing Figure 2a with Figure 2b that the time course of PDE* is not proportional to the time course of T*.
We explore in Figure 2c how the PDE variance decreases with growing number of phosphorylation sites Np (Figure 2c). The maximum and the temporal width of the variance both decrease by increasing the number of phosphorylation sites. Additionally, the variance does peak approximately two-times later than the mean and this feature depends only slightly on the number of phosphorylation sites, as it can be observed in Figure 2d. For six phosphorylation sites, the simulations in Figure 2d are very similar to experimental recordings for the photocurrent presented in Field and Rieke 14.
Figure 2e shows the PDE reliability RP(t) as a function of the normalized mean of PDE*, defined as 
. The value x=1 corresponds to the time to peak of the mean. We decided to plot the reliability ratio RP(t) as a function of the normalized mean, since this provides a better resolution of the rising phase and additionally shows how RP(t) changes as a function of the number of PDE* molecules. The horizontal lines in Figure 2e represent the steady-state values RPs. Our choice of the parameters implies that RPs is maximal and approximately equal to 
, see Eq. (43). It is interesting to note in Figure 2e that during the rising phase, RP(t) reaches values that are much beyond the steady-state value RPs. This apparent paradox is a consequence of the activation dynamics and cannot be anticipated from steady-state considerations.
Finally, in Figure 2f we plot the probability PR(t) (given by Eq. (10)) that R* is activated up to time t. With increasing deactivation steps, R* lifetime becomes less variable and more concentrated around the mean value τ. Moreover, since the decay rate ωλ is small, the lifetimes of the states n are very similar and therefore Rτ is very close to the optimal value 
(data not shown).
Transducin activation rates strongly influence the dynamics of PDE activation
To study the impact of the transducin activation and phosphorylation rates, we present in Fig. 3 and Fig. 4 simulations for toad rods, obtained for various decay rates ωact and ωλ. The number of R* phosphorylation sites is fixed to Np=6, which is the value found in mouse rods and many other species 13,15,28. If R* activity decays only slightly with subsequent phosphorylations (e.g., ωact ∼ 0.1), the PDE* variance peaks nearly twice later than the mean (Figure 3c) and the ratio RP(t) at time to peak is much higher than the steady-state ratio RPs (Figure 3d). In contrast, Figure 4ac, illustrates that the parameter ωλ, which controls the decay of the phosphorylation rates, does not affect much the dynamics of PDE activation, although ωλ strongly influences the reliability of R* lifetime (Figure 4d). We conclude that the behavior of the transducin activation rates is more decisive for the PDE* variance than the variability of R* lifetime.
High activation reliability during the rising phase
We now investigate more closely the time course of the reliability ratio RP(t) during the rising phase. The simulations depicted in Figure 2e and Figure 3d reveal that during the rising phase RP(t) reaches a maximum that can be much higher than the steady-state value RPs. Indeed, this behavior follows from the fact that initially the variance and the mean are almost equal (see Eq. (68)). As a consequence, as long as the variance and the mean are close, RP(t) approximately increases like the square root of the mean and, depending on the number of PDE*, can reach values that are much beyond the steady-state limit. At a later time, the variance becomes much larger than the mean and RP(t) decreases. To show this initial behavior of the variance, we plot in Figure 5a (respectively, Figure 5b) the mean to the variance ratio of PDE* corresponding to the set of parameters used in Figure 2e (respectively, Figure 3d).
To achieve a high reliability ratio RP(t) during the rising phase, it is both necessary that the number of phosphorylation sites Np is large (Figure 2e) and the transducin activation rates are almost constant (Figure 3d). Indeed, the main contributions to the PDE* variance during the rising phase are due to the variability of R* lifetime and the variability of the transducin activation rates. The latter can be reduced by choosing ωact close to zero. Increasing the number of phosphorylation sites reduces the variability of R* lifetime, especially for small times. Figure 2f shows that by increasing the number of phosphorylation sites, a growing initial time-window emerges, during which it is very unlikely that R* becomes deactivated. During this period, the variability of R* lifetime is very low, and in particular much lower than the variability of the whole R* lifetime.
Impact of whether rhodopsin or PDE deactivation limits recovery
In Fig. 6 we present simulations for the mouse rod scenario, where the overall PDE* dynamics is much faster compared to toad rods. In the mouse rod scenario, PDE* deactivation limits the recovery, whereas it is R* shutoff in the toad rod case. By comparing the mouse rod simulations in Fig. 6 with the corresponding toad rod simulations in Fig. 3, we conclude that the dynamics of PDE* activation strongly depends on whether R* or PDE* lifetime limits recovery. Consequently, interchanging PDE* and R* lifetime should affect the overall dynamics, as shown by the simulations in Fig. 7. We now examine more closely these two opposing scenarios (see discussion corresponding to Fig. 6 in 18) where R* deactivation is much faster than PDE* deactivation and then when it is the opposite.
When R* lifetime is much shorter than PDE* lifetime, we can ignore PDE* deactivation during the activation period, and PDE activation and decay occur as two consecutive events. It follows that the peak of the mean number of PDE* occurs at a time when R* becomes deactivated (Figure 6a and Figure 7a) and is given by the steady-state value 
The PDE* reliability ratio RP(t) at the time to peak is given by the steady-state value RPs (see Figure 6d and Figure 7d). Moreover, during the recovery phase, RP(t) is largely constant and close to RPs, which can be understood as follows: Since PDE is first activated and then deactivates, during the recovery phase, the mean and the variance of PDE* are given by a decay process with initial values 
and 
For a timescale smaller than
, we approximate
and thus
The times to peak of the PDE* mean and variance are close (Figure 6c and Figure 7c). If we define the duration of the PDE* response as the time until all the PDE* molecules become deactivated, then, for large 
the mean duration increases logarithmically with the number of PDE* as 
Using some computations, it can be shown that the CV of the duration decreases logarithmically with 
We conclude that the duration of the response becomes more and more reliable with increasing number of PDE*, but the CV of the duration is not zero (see also 18).
When R* lifetime is much larger compared to PDE* lifetime, because PDE activation and deactivation occur simultaneously, unexpected effects are generated. In that case, the PDE* reliability at time to peak can be much higher than the steady-state ratio RPs (Figure 7d). However, since the ratio RP(t) cannot grow faster than 
(see Eq. (70)), the reliability at time to peak remains bounded by 
Consequently, the CV of the amplitude cannot become zero (see 18). Even for constant R* activity kact, a steady state is reached with 
and thus the CV is given by 
If R* deactivation is rate-limiting, the duration of the PDE* response is determined by R* lifetime and the CV of the duration is given by CVτ in Eq. (13), which can be very different from CVareaP. Most of the variability is generated during the recovery phase, which causes the variance to peak much later than the mean (Figure 7c).
Finally, Figure 7d reveals that there is a tradeoff between the reliability during the rising and recovery phase: the higher the reliability during the rising phase, the lower the reliability will be during the recovery phase. To analyze this behavior, we remark that CVareaP is independent of what rate limits recovery (see Eq. (46)) and depends only on the number of phosphorylation sites. Thus, CVareaP is identical in both scenarios presented in Figure 7d. Now, if R* lifetime limits the recovery, the PDE* reliability during the rising phase is high, which implies a low area variability in this phase. Consequently, the reliability of PDE* during the recovery phase has to decrease (which implies a higher area variability in this phase) to ensure the overall value for CVareaP.
Influence of arrestin binding rates
We now examine the impact of linearly and exponentially increasing arrestin binding rates. In the previous simulations, we allowed arrestin to bind only when R* was fully phosphorylated. For a given number of phosphorylation sites, this condition is optimal to minimize the variance. However, experimental results indicate that arrestin already weakly binds before R* is fully phosphorylated. Biochemical data 28 suggested that arrestin binds only to phosphorylated rhodopsin and the affinity increases linearly with the number of phosphorylations. Such a linear behavior was used for photoresponse simulations 18,19. However, experiments 15 have indicated that R* phosphorylation at three sites is needed to trigger arrestin binding with high affinity, which does not imply a gradual increase of the binding affinity. In addition, data obtained from transgenic mice lacking arrestin do not favor a gradual increase of the arrestin binding rates. Finally, there are no specific reasons to favor a linear decay of the arrestin binding affinity, while rhodopsin kinase and transducin affinities show an exponential profile.
To investigate the impact of arrestin binding on PDE* dynamics, we compare in Fig. 8 three arrestin binding scenarios called optimal, exp., and linear, obtained for toad rods with NP=6 and ωact=ωλ=0.1: in the optimal scenario, arrestin binds only when R* is fully phosphorylated. In the linear scenario, the arrestin binding rates μn increase linearly with each phosphorylation step. Finally, in the scenario labeled by exp., the arrestin binding rates increase twofold with every phosphorylation step. To better compare these three scenarios, we chose the arrestin binding rates such that they reach the same maximal rate μN=1.8s−1 when R* is fully phosphorylated. The simulations in Fig. 8 show that a linear increase leads to the highest PDE* variance and the lowest reliability ratio RP(t). This behavior is reasonable, since a linear increase also strongly affects the states before R* is fully phosphorylated. With an exponential increase, arrestin binding rates become predominant when R* is almost fully phosphorylated, while they are relatively weak before. By comparing the simulations in Fig. 8 with the ones in Fig. 2, we deduce that large arrestin binding rates that come up already before R* has been fully phosphorylated have a similar impact to reducing the number of R* deactivation steps. For example, the curves in Fig. 8 for the linear scenario are similar to corresponding ones in Fig. 2 for Np=1. We conclude that for a given number of deactivation steps, linearly increasing arrestin binding rates are not efficient to achieve a high PDE activation reliability.
Only a few activated PDE molecules in cones
In cone photoreceptors, several synchronous photons have to be absorbed 10,11,12 to detect a signal out of the noise. For that reason, it is not possible to estimate experimentally the number of PDE* following a single photopigment excitation. Furthermore, due to experimental difficulties, many fundamental chemical constants are still missing for cones. A modeling approach is thus an unavoidable tool to investigate PDE activation in cones.
The origin of the background noise differs between L- and S-cones 12: in L-cones, a large spontaneous photopigment activation rate 12 constitutes the main source of the noise and this is a direct obstruction of a single photon detection. In contrast, the photopigment of S-cones is very stable and the background noise originates from spontaneous PDE activation 12,22.
Since spontaneous PDE activation is the main source of dark noise in rods and S-cones, we would like to investigate the question of why a single photon response can be observed in rods, but not in S-cones. A possible answer comes from biochemical data 38,39, which suggest that an excited photopigment presumably activates only very few PDE molecules. Biochemical results for carp cones 39 suggest that R* phosphorylation is much faster in cones compared to rods (∼50 times faster), which seems to be caused by a higher rhodopsin kinase concentration and activity. Moreover, experimental data 39 also imply that the transducin activation rates are much smaller in cones compared to rods (∼25 times smaller) and PDE deactivation is several times faster in cones compared to rods. This fast rate can be attributed to the higher RGS9 concentration 40,41.
To estimate the amount of PDE* molecules following a photopigment excitation, we have run various simulations. As expected, we found with no surprise that this amount is a decreasing function of PDE and R* deactivation (Fig. 9). The simulations presented in Fig. 9 are obtained by increasing the PDE* deactivation and R* phosphorylation rates of a toad rod by factors of 5, 10, and 15 (the parameters and simulations corresponding to the toad rod can be found in Fig. 2 for Np=6). We do not alter the transducin activation rates, but smaller transducin activation rates (as suggested in 39) would additionally diminish the amount of PDE* in cones. Figure 9a shows that increasing the PDE* deactivation rate k4 from 1s−1 to 15s−1 decreases the amount of PDE* from 150 (amount for rod) to ∼15. Since R* lifetime is not changed by increasing k4, the recovery of the photoresponse is not affected. In Figure 9b, R* deactivation is enhanced. Compared to Figure 9a, this shows that a faster R* deactivation is less effective in diminishing the number of PDE* molecules. In Figure 9b PDE* deactivation becomes rate limiting (k4=1s−1) since R* lifetime is reduced from 3s to 0.6s, 0.3s, and 0.2s. Finally, in Figure 9c, PDE* and R* deactivation are increased simultaneously, which additionally reduces the amount of PDE*.
From our simulations we conclude that in cones, consistent with biochemical data 39, only very few PDE molecules are activated by an excited photopigment. This result can explain that for S-cones, contrary to rods, many synchronous photon absorptions are needed to produce a signal that overcomes the noise amplitude generated by spontaneous PDE activation.
Discussion
We have studied here PDE activation by a single excited photopigment molecule using a Markov model and obtained explicit equations for the mean and the variance. This approach allowed us to investigate in detail the dynamics of PDE activation, which is indispensable and fundamental for the understanding of the photoresponse in rods and cones. Most experimental recordings are about the photocurrent, and today, unfortunately, there are no direct measurements of PDE activity, which is the main subject here. A full quantitative analysis of the photocurrent will imply to extend the model by including diffusible cGMP.
