The rate equations for cGMP and Ca²⁺ are:

where β∈0,1,

In mouse rod outer segments, cGMP is synthesized by two membrane guanylate cyclases, ROS-GC1 and ROS-GC2, that are regulated by two neuronal Ca²⁺ sensors, GCAP1 and GCAP2, with different affinities for Ca²⁺ (reviewed in Wen et al., 2014). Their activities are described by the terms in brackets on the right-hand side of eq. (1a). The constants αmin<αmax, in μMs−1 are given by the production rates of cGMP at maximum and minimum Ca²⁺ concentration (theoretically as [Ca2+]→∞ and [Ca2+]→0, respectively). Kcyc in μM is the Ca²⁺ concentration for which the rate of production of cGMP by a ROS-GC is half-maximal. The number mcyc is the dimensionless Hill’s exponent. ROS-GC2 pairs primarily with GCAP2 while ROS-GC1 pairs with either GCAP, but for modeling purposes, it is only important to distinguish the overall Ca²⁺ dependent activities. In the last two terms, khyd and kσ;hyd in μm³ s−1 are the catalytic rates of hydrolysis of cGMP by basal and activated PDE subunits, respectively.

The term JcG is the current carried by the CNG channels. The multiplying constant fCa is the dimensionless fraction of cGMP-activated current carried by Ca²⁺. The constant jcGmax, in pA, is the current for maximal [cGMP] (theoretically as [cGMP]→∞), whereas KcG, in μM is the cGMP concentration corresponding to half-maximal channel opening, and mcG is the dimensionless Hill’s exponent. The parameter KcG is itself a function of Ca²⁺ due to the presence of calmodulin as described by the equation:

where mCaM is the Hill’s exponent for calmodulin activation by Ca²⁺, the constant KCaM is the Ca²⁺ concentration for a half-maximal calmodulin effect and KcG,min and KcG,max are the minimum and maximum binding affinities of each of the cGMP binding sites on the channel, occurring for [Ca2+]→0 and [Ca2+]→∞ respectively.

The term Jex is the current due to Na⁺/Ca²⁺, K⁺ exchange, jexsat in pA is the saturated exchange current occurring at maximum [Ca2+] (theoretically for [Ca2+]→∞) and Kex in μM is the [Ca²⁺] at which the exchange current is half-maximal. The constant η is of the form (BCaF)−1, where BCa is a dimensionless constant that takes into account the Ca²⁺-buffering effects in the cytosol and F=96,500 C mol⁻¹ is the Faraday constant. With these specifications, the units in eq. (1b) are self-consistent.

Photon absorption converts R to R^*, which is a substrate for phosphorylation by RK. Phosphorylation decreases the ability of R^* to activate transducin, but complete quench only occurs after the binding of arrestin-1 (Arr). Denote by [Rj∗] the concentrations of activated rhodopsins in their j^th phosphorylated state, identified with the number j−1 phosphates attached to R∗ up to a maximum of six phosphorylations (for mouse). Newly activated R∗ molecules have no attached phosphates, so they are in the state j=1. Rhodopsins in the state j transition either to the state j+1 at a rate of λj or they bind to Arr at a rate μj, thereby being quenched. Removal of phosphates from R^* occurs on a much slower time scale (Berry et al., 2016) and is disregarded in the model. The quantities λj and μj measure, in s−1, the binding rates of Rj∗ to RK and Arr, respectively. The numbers λj+μj−1 measure the statistical average of the sojourn time of an activated rhodopsin in the j^th state. It is assumed that μj of Arr to Rj∗ depends only on the j^th phosphorylated state. For example, μ_j is zero if all phosphorylated states are removed or prevented, e.g., by knockout of RK. For the model, all μj values were set to zero for j < 4 and set to the maximal rate of arrestin binding, μmax, for j ≥ 4. Similarly, λn+1=0 after Rn+1∗ binds Arr (is quenched). The quantities λj are functions of any biochemical process that either increases or decreases the probability of Rj∗ binding RK. For example, a [Ca2+] drop releases RK from its complex with recoverin, thereby increasing the probability of R∗ binding RK and thus shortening the sojourn times before the next phosphorylation. While other dependencies might be present, we assume that λj=λjCa2+ where the form of this function has to be specified. For j≥2, the [Rj∗] is augmented by λj−1Rj−1∗, that is the [Rj∗] imported from the previous state, and is depleted by λj+μjRj∗, that is by the [Rj∗] that transitions to the next state or is quenched. Thus, flashes produce only new [R1∗], which is depleted by (λ1+μ1)[R1∗], and localized at times of flashes by the Dirac mass δt. The rate equations for [Rj∗] are:

The [RK] and hence, the phosphorylation rates of the activated rhodopsin are a function of the [Ca2+] that changes over the course of a flash response. A description for the case of single step rhodopsin quenching and well stirred conditions will be adopted here for sequential, multiple phosphorylations and final quenching. Only the phosphorylation rates are influenced by the calcium concentration drop, whereas Arr binding rate is assumed to be calcium-independent. By extending equation (A12) of Nikonov et al. (2000):

This equation holds true for a single-step deactivation, for which kR=λ and kR,max=λmax. Thus, it is assumed that the same biochemical relationship holds at each state of R∗. From equation (C1) of Nikonov et al. (2000):

According to equations (C2) and (C4) of Nikonov et al. (2000):

Then the dependence on [Ca2+] comes from the algebraic second order equation (C3) of Nikonov et al. (2000):

Here Ki for i=1,…,4 and M are biochemical parameters defined in Nikonov et al. (2000).

During its average lifetime τR∗, photoexcited rhodopsin R^* activates many transducins (T→T∗) by catalyzing GTP exchange for the GDP bound to their α-subunits. Each molecule of T^* associates, one-to-one, with a catalytic subunit of the effector forming a T^*--E complex. Full activation of the catalytic subunit is assumed, denoted by E^*. A single molecule of E^*, during its average lifetime τE∗, hydrolyzes over 50 molecules of the second messenger cGMP which, dissociating from the cationic channels they keep open, causes channel closure and thereby suppression of the inward current JcG (Pugh and Lamb, 1993 2000). T^* possesses an intrinsic GTPase activity that terminates T^*--E activity. The synthesis of cGMP overtakes hydrolysis and as the levels rise, CNG channels reopen, ending the time in saturation, Tsat.

Let [T]₀ and [T^*]denote the initial, basal concentration of transducin in the rod outer segment (ROS) and the concentration of activated transducin, respectively. Also, let [E]₀ and [E^*] denote the concentrations of subunits of the effector PDE in darkness and of subunits of activated effector PDE^*, respectively. Assuming independent activation of each PDE subunit, [PDE]0=12[E]0 and[PDE∗]=12[E∗], the rate equations for [T^*] and [E^*], for a well-stirred activation/deactivation model are:

The constant kT∗E in μm² s−1, is the rate of formation of the complex T^*--E. It is assumed that activation of a molecule of transducin occurs upon encounter with a molecule of activated rhodopsin. Thus, the first term on the right-hand side of eq. (9a) measures the increase of [T∗] due to activated rhodopsin [R∗]. The second term measures its depletion due to its binding to PDE. The constant kE* = 1/τ_E, where τ_E is the average E^* lifetime dictated by the deactivation rate of the T^*--E complex. It is remarked that eqs. (9a, b) represent a simplified model of the disc cascade widely used in the literature. In order to study the behavior of rod cells under a high saturating illumination, a more detailed model is proposed and used for the numerical simulations, as reported in eqs. (21a-f), explicitly accounting for the RGS9 dynamics.

Suppose that a bright flash Φ (in instantaneous number of isomerizations) applied at time t=0 causes cGMP levels to fall below the minimum required for CNG channel opening. Then the output current drop saturates for a time, Tsat (Figure 1). In practice, there are technical difficulties in measuring Tsat; bright flash responses do not have sharp transitions into or out of saturation on the rising and recovery phases. Although JcG terminates abruptly, the decline in Jex carried by the sodium/potassium/calcium exchanger follows a slower time course and there is noise that is omnipresent in electrical recordings. To circumvent these difficulties, saturation time is typically measured from mid-flash or from flash onset to a criterion recovery of the response. Such a practice inflates Tsat, but the upward shift does not alter the slope of the linear portion of the Pepperberg curve (e.g., Figure 2). Various publications have adopted different criterion recovery levels, so care must be taken in comparing the X-intercepts across studies.

Tsat increases linearly with LnΦ, with slope τ_D that is determined by the shutoffs of R^* and/or T^*--E (Pepperberg et al., 1992, 1994). It was shown that the shutoff of T^*--E was slower than that of R^* (Krispel et al., 2006), so assuming that τR∗ <<τE∗, the value of τE∗ can be determined from the slope of the Pepperberg curve. Hence,

for a constant C > 0 that has not been described previously and will be discussed in further detail below.

In deriving the Pepperberg curve equation, a well stirred condition was assumed, i.e., all concentrations were taken to be uniform inside the ROS. Moreover, it was assumed that the rhodopsins were all of one species denoted by R^*. With this stipulation only eq. (4a) was in force and

the average lifetime of R^*. For a saturating flash Ф applied to the system at time t = 0,

It was further assumed that the effects of the intermediate activation of transducin can be lumped into a direct activation of E^* by R^*, then by a modified form of eq. (9b):

Integrating over 0,Tsat where Tsat is a positive number to be defined:

for τR∗<<τE∗. From this:

We notate Echar∗=[E∗](Tsat). It is a fundamental quantity that characterizes an essential feature of the system. To define Tsat and compute Echar∗, refer back to eq. (1a), viewed at times of saturation, originating from the dark-adapted state. Stipulate that the interval (0, Tsat) lies within the saturation regime, and that in such a regime, [cGMP] and [Ca²⁺] are in quasi steady state. Thus, we may set approximately to zero the time-derivative in eqs. (1a) and (1b). Since khyd<<kσ,hyd, the term khydE∗cGMP can be neglected. It then follows from eqs. (1a) and (1b):

This relation holds for all times during which saturation is in force. For the same times, eq. (16b) yields Jex=12fCaJcG. Therefore, the total photocurrent is

Let Jmax be the maximum current that the ROS can output (typically, Jmax=Jdark). Then define Tsat as the time for which

Eqs. (16b) and (18) can be solved for the unknowns [cGMP] and [Ca²⁺] at time Tsat. This permits the determination of Echar∗ by eq. (16a). Inserting it in eq. (15) gives a linear relation between Tsat and LnΦ, with slope τE∗, of the form

The quantity Φo can be identified as the flash strength for which Tsat=0. Measurement of Φo could permit the determination of Echar∗.

The constant C in eq. (10), is mathematically derived above as

It is emphasized that C is proportional to τE∗ and independent of Φ. This form of C implies that the straight lines described by eq. (19) corresponding to different values of τE∗ have a common point (i.e., Φ=Φoand Tsat=0) regardless of their slope. This feature actually occurs in the experimental curves of Krispel et al. (2006), reported in (Burns and Pugh, 2009) and is discussed below.

Outputs from the FSR model and the GWS version converge at very bright flashes. Since usage of the GWS version greatly facilitates numerical simulations, we ascertained whether it would be adequate for responses even to the weakest saturating flashes by computing bright flash responses for mouse rods using both models with the parameter set for mouse in the Supplementary Table S1. Figure 2 compares the Pepperberg curves; the close match of the curves verifies that even for two flash strengths that do not quite saturate the response, any spatial inhomogeneity of cascade activity was not distinguished by the two versions of the model. Hereafter, all simulations were made using the GWS model.

Since experimentally determined values for Tsatare always measured after some arbitrary criterion recovery, they are inflated relative to the true saturation times and the X-intercept is shifted to the left, underestimating the true values for Φo and E^*_char that describe the minimal conditions for saturation (see above). For modeled responses in Figure 2, we calculated that the Tsat values determined at 20% recovery were increased by 0.46 ms, compared to those determined at 2% recovery and that this had the effect of translating the X-intercept along the abscissa from a LnΦ of 3.3 to 5.3, i.e., by a factor of 7.4-fold. Such a correction could be applied to the Pepperberg plots in the figures below for which Tsatwas measured from mid-flash to 20% recovery, to improve the estimates of Φo and E^*_char. The magnitude of the translation would have to be adjusted for determinations of Tsat made at different criterion levels of response recovery.

Once activated, T relieves an autoinhibition of PDE, allowing it to hydrolyze cGMP. The GTPase activity of T^* terminates the interaction between T^* and PDE. The intrinsic GTPase activity is slow, but is greatly enhanced by the encounter with an RGS9 complex. By showing that an increase in the expression level of the RGS9 complex accelerated flash response recovery and reduced the slope of the Pepperberg relations in rods, Krispel et al. (2006) validated the idea that the collisional delay of RGS9 complex with T^*--E was the major determinant of T^* lifetime and identified the shutoff of T^*--E as the rate-limiting step in the response recovery. Attempts to accelerate R^* shutoff by overexpression of RK did not alter the slope of the Pepperberg plot. Our model captured the Pepperberg relations for the four types of rods described in Burns and Pugh (2009), each expressing a different level of RGS9 (Figure 3). The values for τE∗ from our model were: 792 ms for rods expressing 0.2x the normal amount of RGS9 complex, 381 ms for WT, 209 ms for 2-fold overexpression, and 191 ms for 4-fold overexpression. The values differ somewhat from the population averages of: 246 ms for WT, 108 ms for 2-fold overexpression, 80 ms for 4-fold overexpression obtained from linear fits over the lower range of flash strengths (Krispel et al., 2006) but conserve the overall trend.

Since the timing of T^* shutoff (i.e., τE∗) is too slow to affect the rising phase of the bright flash response, rods of each expression level should have reached saturation at a similar flash strength and all Pepperberg curves should have converged at the same X-intercept, Φo_. The curves nearly converged, although the intersection fell slightly below Tsat = 0. We attribute the deviations of the experimental results from the model to biological variation and to technical issues associated with achieving the same collecting area for every rod during recording.

Burns and Pugh (2009) proposed a model for the experimentally obtained Pepperberg results of rods expressing various levels of RGS9 complex that consists essentially of eq. (9b), where kE is replaced by kfRGS9 for a suitable constant kf to be computed by fitting. Since τE∗ is the reciprocal of kE, the larger [RGS9] the smaller τE∗ and vice versa. A key assumption of that model is that [RGS9] remains constant over the range of light intensities that they analyzed. While the fit is good for lower saturating flash intensities (less than ~3,000 photoisomerizations), saturation times depart from the curves at higher flash intensities, becoming non-linear and exhibiting an upward, convex, super-linear bending (Figure 3). It would appear that a new mechanism dominates the response recovery at these flash strengths.

Because the convexity became more prominent in rods underexpressing the RGS9 complex, and was absent from rods overexpressing it, we propose that [RGS9] does not remain constant throughout the process. While the total mass RGS90 is constant, it dynamically subdivides into a portion [RGS9--T^*--E] bound to the complex T^*--E and a portion [RGS9₀] − [RGS9--T^*--E] available for binding to T^*--E according to the rate equations:

with initial conditions

where T0,E0,andRGS90 are the initial, basal concentrations of the transducer G protein, the effector PDE and GTPase-activating protein RGS9. Here νRT is the catalytic activity of R^* to T^*, and kT*E is the coupling coefficient between T^* and E^*. Also kf (forward) is the rate of association of the complex T∗−−Eper unit mass with available RGS9, and kb (backward) is the rate of dissociation of RGS9 from the complex T^*--E. Finally, kcat is the rate of deactivation of the RGS9--T^*--E complex by the hydrolysis of GTP by T^*.

The dynamics of RGS9 along its time evolution after the flash was traced from its basal values and was observed to lag slightly, the dynamics of T^*--E for flash strengths LnΦ<8 (Figure 4). However, for brighter flashes there is an upward, convex, super-linear bending of the Pepperberg relations that the model attributed to a reduction in RGS9 complex availability. According to the model, T^* shutoff was delayed by the extra time needed for an RGS9 complex to collide with T^*--E under these conditions and in the extreme case, would theoretically approach a regime in which shutoff relied on transducin’s intrinsic GTPase activity, unassisted by RGS9.

Reducing kb slows the dissociation of RGS9 from the complex RGS9--T^*--E and hence, for flashes of equal intensity, the time in saturationTsat was shorter. Increasing kf sped up the formation of the complex RGS9--T^*--E, accelerating the shutoff of T^*--E, and hence, for flashes of equal intensity, Tsat was shorter. Opposite effects occurred by increasing kb or by decreasing kf. In all cases, for LnΦ<8, their behavior is essentially that of a straight line, and in all cases, one observes a convex, superlinear bending with brighter flashes. To provide further evidence that such a pattern is due to RGS9 depletion, the red curve in Figure 3 was obtained by artificially augmenting the basal value [RGS9]₀ by a factor of 4 and at the same time by reducing the “forward” parameter kf by the same factor. This way, at incipient phases of the process (t≈0), when only a negligible fraction of [RGS9]₀ has been turned into [RGS9--T^*--E], the product

remains essentially unchanged. Therefore, for low values of [RGS9--T^*--E] the complex [T^*--E] dissociates and the complex [RGS9--T^*--E] is generated at the same rates as in the blue curve in Figure 3. A departure occurs at later times as [RGS9]₀ is depleted, by increasing [RGS9--T^*--E].

For even brighter flashes, LnΦ>15,, the model predicted an asymptotic flattening in the Pepperberg curve caused by activation of all available PDE (Figure 5). Full activation of PDE achieves the greatest reduction of cGMP, hence it sets an upward limit on Tsat. To further explore the role of PDE expression, theoretical changes in PDE level were imposed on the model. A tenfold decrement in PDE content actually prolonged Tsat at all flash strengths, because there was a reduction in basal PDE activity that elevated cGMP level in darkness, extending the time required for E^* activity to return to the dark level after a bright flash (gold trace, Figure 5, explained further in Nikonov et al., 2000 and commentary by Govardovskii et al., 2000). Interestingly, the upward bending from linearity in the Pepperberg curve appeared at a reduced flash strength and rose more steeply. The basis was that elevated cGMP supported a greater fraction of open CNG channels and a greater influx of Ca²⁺. As a result, the content of Ca²⁺-bound recoverin was higher and a greater proportion of RK was not available to initiate shutoff of R^* until there was a sufficient light-induced fall in Ca²⁺. Thus, even though PDE levels were reduced, each R^* was able to generate more T^*--E than normal. Reduced expression of PDE along with an increased number of PDE^* per R^* caused the Pepperberg curve to fully activate all PDE and approach the plateau with less intense flashes than for WT rods. Consistent with this explanation, the upward bending occurred with more intense flashes in mutant mice lacking recoverin (Makino et al., 2004) and occurred with less intense flashes in mutant mice expressing lower levels of RK (see below).

Under these conditions, the model was used to interrogate whether the upward bend reached its maximal rate of rise of 9–10 s, the time taken for the response recovery in the absence of RGS9 complex (Chen et al., 2000; Krispel et al., 2003a; Keresztes et al., 2004). The maximal slope within the range of the upward bend (convexity) was nearly 4 s, indicating that the maximal rate of rise was not achieved. Since the upward curvature rises less steeply for WT rods, we predict that the ratio of expression of RGS9 to PDE is adequate to ensure that T^*--E shutoff would never rely on the T^* GTPase activity free from RGS9 acceleration at any flash strength.A theoretical, tenfold increment in PDE levels had minor effects in the opposite direction (red trace, Figure 5): Tsat was somewhat shorter at every flash strength, the upward bending and the plateau occurred at slightly greater Ln(Φ) values.

Experiments addressing the plateau in the Pepperberg plot are not yet available. Challenging technical issues and slower mechanisms of light adaptation arise with the recording of responses that remain in saturation for such extended periods to extremely bright flashes, but in principle, it could be done. An exciting recent development is that a mutant mouse has been generated with greatly reduced expression of PDE (Morshedian et al., 2022), that could be useful for testing the model predictions.

A reduction in RK expression slowed R^* shutoff, indicating that in large part, collision time between R^* and its kinase dominates the process. Hemizygous knockout of RK reduced RK levels by 70%, increased saturation times and shifted the X-intercept to lower flash strengths, however, it did not change the slope of the Pepperberg significantly (Figure 6, see also Sakurai et al., 2011; Gross et al., 2012). Upward bending in the RK+/− Pepperberg curve occurred at a lower flash strength than in WT rods (Sakurai et al., 2011), because slower R^* shutoff resulted in activation of more transducins, and a greater amount of T^*--E that depleted more RGS9 at each flash strength.Overexpression of RK by two-fold or even four-fold had little effect on saturation behavior (Krispel et al., 2006; Sakurai et al., 2011). However, expression of a mutant S561L RK on a WT background reduced saturation time without a change in slope (Figure 6, see also Gross et al., 2012). The total level of both types of RK was 8.7-fold higher than normal and the mutant form contained a sequence that specified geranylgeranylation instead of farnesylation to enhance membrane affinity and quicken R^* phosphorylation. The constancy in slope in the face of changes in the rate of R^* shutoff helped to advance the argument that the rate-limiting step in flash response recovery was the shutoff of T^*--E.

Mouse rods normally express two GCAPs that differ in their K_1/2 for Ca²⁺: 46–47 nM for GCAP1 and 133 nM for GCAP2 (Makino et al., 2008, 2012). According to the model, Echar∗ changes with ROS-GC activity because the balance between cGMP synthesis and hydrolysis by E^* sets the cGMP levels and thereby affects response saturation. Knockout of both GCAPs removes Ca²⁺ feedback onto ROS-GC and increases the size of the single photon response five-to six-fold by allowing for a greater drop in cGMP after the photoisomerizations (Mendez et al., 2001; Burns et al., 2002; Chen et al., 2010; Makino et al., 2012). The model captured the increases in saturation times without a change in slope in the Pepperberg relation, because τE∗ remained the same (Figure 7B). Because the effect of T^*--E was opposed by a lower ROS-GC activity, Φo shifts to less intense flashes. Similarly, knockout of both GCAPs from mice expressing varying levels of RK had the same effect of increasing saturation time without a change in slope (Figure 7C). Importantly, the proportional increase in saturation time was similar for all three types of rods (cf. Figure 6).

Knockout of GCAP2 has a less dramatic effect as the maximal rate of cGMP synthesis at low Ca²⁺ drops to 40% and the overall K_1/2 for Ca²⁺ shifts from 133 nM to ~47 nM (Makino et al., 2008). The slope of the Pepperberg is unchanged, but the X-intercept shifts to a lower value because of the response prolongation due to the reduction in Ca²⁺-dependent ROS-GC activity (Figure 7A).

Mutant mouse rods expressing a mutant CNG channel that lacks the calmodulin binding site have reduced saturation times and a subtle but significant decrease in slope, τ_D, with no change in dark current, sensitivity or dim flash response kinetics (Chen et al., 2010). These experimental findings were not predicted by the model (Figure 8). Ca²⁺/calmodulin lowers CNG channel affinity for cGMP in darkness so in the absence of calmodulin modulation, channel affinity for cGMP increases and dark current increases somewhat. During the response to a saturating flash, Ca²⁺ drops to a minimum so calmodulin no longer reduces the apparent affinity of the channel for cGMP, hence WT and mutant channels would reopen at the same concentration of cGMP. As cGMP levels recover and outer segment levels of Ca²⁺ rise, fewer WT channels will open compared to mutant channels. The steeper trajectory of the mutant rod response recovery would tend to shorten Tsat, but because the saturating response amplitude is larger in mutant rods, measurement of Tsat to 20% recovery would be made after a longer delay. The net result in simulations was a slight extension of Tsat at all flash strengths with the loss of CNG channel regulation by calmodulin, leaving the slope of the Pepperberg relation unchanged. The discrepancy between the experimental results and the predictions of the model indicate that other factors must come into play. The absence of calmodulin binding might indirectly affect regulation of the CNG channel in other ways, such as by affecting its phosphorylation state or by altering its interaction with Grb14. Knockout of Grb14 has the effect of reducing Tsat and τ_D (Woodruff et al., 2014). Another possibility is that levels of free calmodulin in the outer segment might be changed and impact other targets such as CaM kinase.