Local, nonlinear effects of cGMP and Ca2+ reduce single photon response variability in retinal rods

Summary

The single photon response (SPR) in vertebrate photoreceptors is inherently variable due to several stochastic events in the phototransduction cascade, the main one being the shutoff of photoactivated rhodopsin. Deactivation is driven by a random number of steps, each of random duration with final quenching occurring after a random delay. Nevertheless, variability of the SPR is relatively low, making the signal highly reliable. Several biophysical and mathematical mechanisms contributing to variability suppression have been examined by the authors. Here we investigate the contribution of local depletion of cGMP by PDE*, the non linear dependence of the photocurrent on cGMP, Ca²⁺ feedback by making use of a fully space resolved (FSR) mathematical model, applied to two species (mouse and salamander), by varying the cGMP diffusion rate severalfold and rod outer segment diameter by an order of magnitude, and by introducing new, more refined, and time dependent variability functionals. Globally well stirred (GWS) models, and to a lesser extent transversally well stirred models (TWS), underestimate the role of nonlinearities and local cGMP depletion in quenching the variability of the circulating current with respect to fully space resolved models (FSR). These distortions minimize the true extent to which SPR is stabilized by locality in cGMP depletion, nonlinear effects linking cGMP to current, and Ca²⁺ feedback arising from the physical separation of E* from the ion channels located on the outer shell, and the diffusion of these second messengers in the cytoplasm.

Keywords:

Simulation and modeling – Cytoplasm – Photoreceptors – Mass diffusivity – Photons – Salamanders – Second messenger system – Phototransduction

Introduction

Vertebrate rod photoreceptors accurately detect light and reliably discriminate differences at exceedingly low levels of illumination. The biochemical cascade that transduces photons into integrated electrical signals that lead to changes in neurotransmitter release at the synapse is inherently stochastic. A photon, absorbed by the 11-cis-retinal covalently attached to rhodopsin, converts it into all-trans-retinal, which forces rhodopsin into an active state (R → R*). Active rhodopsin R* can be localized anywhere on the rod outer segment disc membrane (Fig 1), and it continues random diffusion in the membrane, encountering and activating a variable number of transducin G-protein (T) molecules by facilitating the exchange of GDP bound to inactive heterotrimeric transducin for GTP. Binding of GTP to the α-subunit of T induces its dissociation from active rhodopsin and from its βγ-subunit leaving R* free to activate additional molecules of transducin. The GTP-liganded Tα (T*) physically interacts with a cGMP phosphodiesterase and activates it (E → E*). Active E* then hydrolyzes a variable number of cGMP molecules until it is deactivated as a result of the hydrolysis of GTP to GDP by the intrinsic GTPase activity of Tα which is facilitated by both E ([1]) and RGS9 protein ([2]). Transducin deactivation is also a stochastic process, which makes E* lifetime stochastic. The resultant drop in cytoplasmic cGMP concentration leads to the closure of hundreds of cGMP-gated (CNG) ion channels on the plasma membrane, preventing the entry of Na⁺ ions into the rod. Channel closure also prevents the influx of Ca²⁺, whereas its efflux by a Ca²⁺-exchanger continues. This leads to a drop in cytoplasmic Ca²⁺ concentration ([3]). Calcium dissociates from guanylyl cyclase activating proteins (GCAPs) to be replaced by Mg²⁺, which converts GCAPs from inhibitors to activators of guanylyl cyclase (GC) ([4]). Increased GC activity then replenishes cytoplasmic cGMP. Synthesis of cGMP by GC is governed by a feedback system based on Ca²⁺. Activity is slow when Ca²⁺ is high in the dark, but accelerates when Ca²⁺ levels fall during the light response. As a result the channels reopen sooner, with a consequent rise of Ca²⁺, and eventual GC inactivation by Ca²⁺-liganded GCAPs. During response recovery R* is rapidly phosphorylated by rhodopsin kinase ([5]). Active rhodopsin R* can acquire as many as six (in mice and humans) or seven (in cows) attached phosphates before it encounters arrestin-1 by diffusion and is completely inactivated by arrestin-1 binding. Progressive rhodopsin phosphorylation reduces the efficiency of transducin activation ([6]). Thus R* exists for variable times in each of the phosphorylation states with different activities and is stochastically shut off by arrestin-1, typically after acquiring three or more attached phosphates.

The stochastic nature of every step in this biochemical cascade generates significant variability. Yet the rod’s single photon response (SPR) is significantly less variable that one would expect ([7]). Different mechanisms were proposed to underlie this unexpected reproducibility of SPR ([7–13]). In ([13]), a fully space resolved (FSR) biophysical model of phototransduction in rods, populated by experimentally tested parameters ([13, 14]), was used to identify the key factors that generate variability of the SPR and those that suppress it. Variability arises primarily from the random steps governing deactivation of the cascade, which generate differences in the spatial distributions and lifetimes for the activated Transducer-Effector complex, T*-E*. The ensuing diffusion of second messengers is deterministic. Variability of the SPR is reduced by (i) the localized depletion of cGMP within the rod outer segment (ROS) by T*-E*, (ii) the subsequent diffusion of cGMP and Ca²⁺ in the cytosol of the ROS, and (iii) the nonlinear relations linking ionic current to [cGMP] and [Ca²⁺]. The first two factors are intricately tied to the complex shape of the outer segment, which is conserved across species (Fig 1). The ROS houses a stack of disc structures, whose upper and lower membrane surfaces contain receptor rhodopsin R, G-protein transducer T and PDE effector E. While these components are mobile on the disc surface, they do not “hop” from disc to disc. Hence their spatial positions vary over time in transverse directions, but are fixed in the longitudinal direction. In contrast, the soluble cGMP substrate of E* moves throughout the cytoplasm of the ROS communicating E* activity on the surface of the activated disc, to CNG channels on the plasma membrane. The discs pose a barrier to movement of cGMP, but they typically contain one or more incisures that facilitate longitudinal diffusion of soluble second messengers.

Because of the small size of the transversal cross section of the ROS, with respect to its length in most species, it is tempting to assume that transversal cGMP diffusion plays a negligible role in the cascade, all the more in the presence of a large diffusion coefficient D_cG, which would favor rapid transversal equilibration ([15, 16]). The relatively small changes in [cGMP] that occur during the SPR invite a further simplification, linearization of the relation between [cGMP] and ion channel activity ([15, 16]). We disprove such approaches by providing numerical evidence that the full 3-dimensional structure of the ROS and the local nature of the cascade activation-deactivation and nonlinearities in the effects of second messengers, play key roles in suppressing variability. To underscore the role of localization and to extract the contribution of diffusion, the simulations were effected for small diameter mouse rods and for large diameter salamander rods, and for several values of D_cG. The role of Ca²⁺ feedback was evaluated by “clamping” [Ca²⁺] at the dark level. The suppressive effects of locality, nonlinearity and Ca²⁺ feedback on SPR variability change in significance over time, an aspect that is largely ignored by the functionals usually used to assess SPR variability. Here we discuss what information is provided by different functionals that evolve over time and use them to show how variability is distorted by reducing spatial resolution and by linearization of second messenger effects.

Materials and methods

The current densities J_ex and J_cG due to Ca²⁺ exchange and the cGMP-gated channels, respectively, are given by

Measuring the variability of the single photon response

Variability of the SPR may be measured by the coefficient of variation (CV), defined as the standard deviation over the mean of a pre-chosen functional:

The relevance of the information contained in the CV of these functionals may change according to the way the signal is processed downstream of the cascade. At the limit of sensitivity, the bipolar cell synapse ignores low amplitude phototransduction noise but may be saturated by the rod SPR. In this case, variability in SPR amplitude may not have functional consequences. The bipolar cell response is also faster than the rod response; the bipolar cell response is largely complete by the time the rod response reaches its peak, so it would seem that rod response recovery is not significant either. Since the functional I_area is dominated by the recovery phase of the rod SPR, CV(I_area) might not be that important. On the other hand, Field et al ([21]) suggest that saturation at the bipolar cell may occur with more than one photon, in which case rod SPR amplitude I(t_peak) and integration time I_int(t) as well as their variability become important. It is possible that both scenarios are correct under different adaptation regimes; Taylor et al ([22]) may be describing the situation under the most dark adapted conditions, whereas Field et al ([21]) may be looking at it under very slightly light adapted conditions. To complicate matters further, rods are electrically coupled to other rods and cones, enabling the SPR from one rod to spread across the retina. But like ripples in a pond disturbed by a stone dropping into the water, the signal diminishes with distance from the rod generating the SPR. So under slightly light adapted conditions, SPR amplitude and area and their variability again become important. The CV of the functional j_tot(t), to the best of our knowledge has not been used in the literature. Yet it seems to be a relevant functional for the reasons we present below.

Let X be the probability space of events of I(t) with probability measure dω. Then from Eqs 3 and 4

A further justification for considering CV of J_tot(t) is that transmission of the SPR at the synapse, downstream of the cascade, is concerned with voltage change, which depends upon j_tot(t) and not the relative current drop I(t). The latter is a contrived way to avoid referring to a decrease in an inward cationic current in response to photon absorption. We will use the dynamics of the functional t → j_tot(t) and its CV in the context of Ca²⁺ clamped virtual experiments.

Diffusion of cGMP and Ca²⁺ in the cytoplasm

Following ([11, 23, 24]), in the cytoplasm [cGMP] satisfies the diffusion equation

Both Eqs 9, 10, 11 and 12 might not have a pointwise mathematical meaning and must be interpreted in a weak-integral form. They also need to be complemented by similar processes on the incisures and coupled equations on the lateral boundary of the ROS. A full, rigorous formulation of the model is in A.2—A.3 of S1 Appendix.

The activation-deactivation cascade

Molecules of T* and E* diffuse on the activated disc by the random walk t → x(t) (dimension μm⁻²) of R* on the activated disc D from which the limiting incisure has been removed. Following activation, R* becomes deactivated after an exponentially distributed random time t_R*, of mean τ_R*. During the random interval (0, t_R*), R* evolves through n molecular states, produced by sequential R* phosphorylation by rhodopsin kinase, R j *, j = 1, …, n, each with transducer-activation rate ν_j, with random transition times t_j−1 < t_j ≤ t_R*. Thus the rate equations for the surface densities [E*]_σ and [T*]_σ are

These equations also need to be mathematically interpreted in a weak form (A.4 of S1 Appendix). The FSR model, its mathematical weak formulation, and its biological validation have been introduced and refined in a series of contributions ([11–14, 23–27]). The main equations are reproduced here in a pointwise form, to stress that the only source of variability in the cascade is the surface dynamics of T* and E* on the activated disc, as expressed by the system of Eq 13. The function [cGMP] experiences randomness only through the random term [E*]_σ, which appears on the right-hand side of Eq 9. Then [cGMP] is computed and put to use sequentially in Eqs 2–5.

Let E*(t) denote the total number of molecules of activated effector PDE* at time t, downstream of the activation-deactivation cascade, i.e.,

These functionals parallel at the activation/deactivation level, the functionals I_int(t), I_area, I(t_peak), and j_tot(t) at the response level, as defined in Eqs 3–5. Notice that t peak ≠ t peak *. Randomness of the experimentally measured current suppression I_int(t) and I_area is indirectly imported from E* through [cGMP], by the Eqs 9 and 2. To separate these two levels of randomness in our simulations we report the time dynamic of the CV ( E int * ( t ) ) alongside with the CV(I_int(t)), to highlight how the latter reflects a variability reduction of the former. Numerical experiments are performed on mouse and salamander rods, for which there are complete and consistent sets of parameters ([11–14]) (see B of S1 Appendix). The cGMP diffusion coefficient for mouse was taken as D_cG = 120 μm²/s, close to the reported experimental value D_cG ≈ 140 μm²/s ([28]). In our mouse simulations we also tested a theoretical value, D_cG = 330 μm²/s suggested in ([29]) and used in simulations for the transversally well-stirred model of ([15, 16, 30]). By the conversion formula of ([24]) (D of S1 Appendix) these correspond to longitudinal diffusion coefficients of ≈ 14 μm²/s and ≈ 40μm²/s. Using D_cG = 330 μm²/s and keeping the remaining parameters unchanged reproduces the experimental SPR presented in ([12]), for the choice ν_RG ≈ 230/s (C of S1 Appendix). This new value yields new catalytic activities ν_j for R* in its jth phosphorylated state by the formula ν_j = ν_RG exp{−k_v(j − 1)} ([12]). For the salamander simulations we include 23 radially equally spaced incisures assimilated to right circular sectors with “base” 15nm set on the rim of the discs, radius 4.64 μm, and angle 0.015/4.64 radians. The total area of the incisures is 0.8 μm².

Virtual experiments

Fully space resolved versus transversally and globally well stirred models

Our FSR model takes into account all geometrical aspects of the ROS. The TWS model assimilates the ROS to a segment of length H, thereby disregarding the transversal dynamic of the players, and more importantly the intricate relationship between interior and boundary dynamics. While at times used, in view of its mathematical simplicity ([15, 16, 31]), no justification is provided other than it is “generally accepted” because of transversal “rapid equilibration”. The GWS model removes the geometry of ROS altogether and describes the cascade only as a sequence of mass action relations of the various players. The TWS and GWS models can be derived by Eqs 1–13, and their weak formulations (A.3 of S1 Appendix) by progressively removing the geometrical features of the ROS. The limitations of such lumped models have been presented in ([11, 13, 23, 24]). The GWS model is equivalent to letting D_cG → ∞ in the FSR model, so the effect of diffusion was explored by testing several values for D_cG. The TWS model is approached by the FSR as disc diameter is reduced to zero, so to assess the effect of disc size separately, we simulated mouse and salamander rods, because their disc sizes span the range found in nature. In addition these are the two species for which the most knowledge has accumulated.

Local cGMP suppression and variability by deterministic simulations

Some initial simulations were carried out in which R* deactivation was deterministic to compare the response for short and long R* lifetimes. Photoexcitation of R was always taken as the center of the disc located in the middle of the ROS. The relative, local cGMP depletion and its average on the outer shell are

Analyzing variability by stochastic simulations: Linearization and local cGMP depletion

A standard numerical WT experiment consisted in selecting n sojourn times s_j according to their exponential distributions with means τ_j, putting them in Eq 13, from which, by a finite elements, Matlab based code ([32]), [E*]_σ was computed and fed into Eq 9, and its counterpart for [Ca²⁺], to compute [cGMP] and [Ca²⁺] as functions of space and time, at the lateral boundary of the ROS. Finally, using these [cGMP] and [Ca²⁺] values, current and current-functionals were computed from Eqs 1–5. It should be stressed that in Eqs 1 and 2, and hence the subsequent equations, only the boundary values of [cGMP] and [Ca²⁺] were used, since the lateral boundary of the ROS is where the current is generated. The process was repeated 1,000 times and from these outputs CV was computed for various E int *, I_int, I and j_tot, and their linearized counterparts. Simulations were carried out for the values of volumic diffusion coefficient D_cG = 120 μm²s⁻¹ ([13]) and D_cG = 330 μm²s⁻¹ ([15, 16]).

A set of virtual KO simulations, whose results are labelled by LIN (Linearization), was effected by the same steps and procedures where now the cGMP and Ca²⁺ fluxes Eqs 10 and 12 where linearized. A second set of virtual KO simulations, labelled NLD (Non Local Depletion), was effected by eliminating the local cGMP depletion by E* during deactivation, i.e., by keeping [cGMP] equal to its dark value [cGMP]_dark in the last term of Eq 9. A combination of these virtual KOs denoted by LIN+NLD enforces both effects.

Computations were then repeated using the TWS and GWS models to extract the implications to variability of various degrees of space resolution.

Results and discussion

Localization, nonlinearity and variability

Locality and nonlinearity inform variability at multiple levels. First, variability of E* passes to cGMP by means of the last term on the right-hand side of Eq 9. Second, variability of cGMP depletion inherited from E* gets deamplified by the nonlinear dependence of ion channel activity on [cGMP], by Eq 9. Between these two steps is the redistribution of cGMP, as it diffuses centrally from the outer shell, dampening the randomness of cGMP at the outer shell. Channel closure leads to a local decrease in [Ca²⁺], causing GC to replenish the cGMP at the outer shell, further dampening the randomness of cGMP there.

A first approximation of variability can be gained by comparing SPR simulations for deterministic short- and long-lived R*, e.g., deactivation with lifetimes of τ_av and 2τ_av, as indicated in §. The relative difference of the outputs, Δ%, can be taken as a first variability estimator. For I_int(t) and I_int;lin(t) these values were computed at time t = t_peak. The results are summarized in Tables 1 and 2 for mouse and Tables 3 and 4 for salamander. In all cases, Δ% of the functionals E max * and E area * exceeded that of the corresponding functionals I(t_peak) and I_lin(t_peak) reflecting the incomplete transfer of variability from E* to I(⋅), as explained above.

Mouse dynamics: Linearization and variability

Tables 1 and 2 shows the relative errors in passing from the nonlinear relation Eq 2 to the linear relation Eq 17. These errors were in all cases at least of the order of 10% and were dramatically larger with lower diffusivity (D_cG) and longer R* lifetime, casting doubt on the validity of the linearization Eq 17. Essentially, the computed SPR amplitude is larger with linearization and the discrepancy increases with the extent of cGMP depletion. Consequently, SPR variability was overestimated by linearizing the relation between current and [cGMP] and/or by increasing the diffusion coefficient D_cG and/or by by progressively disregarding the spatial resolution of the ROS, from FSR to GWS. These trends were explored further and validated by more rigorous determinations of CV. As was suggested by deterministic simulations, the distortions in computing the current suppression I_int(t), by its nonlinear form Eq 2 and its linearized form I_int;lin(t) as given in Eq 17 increased with longer R* lifetimes, hence they gave rise to errors in variability. While mitigated by averages, these errors persisted in the computation of the CV of current suppression. Fig 4 traces the dynamics of t → CV[I_int(t)] for mouse rods, both for D_cG = 120 μm²/s (Upper Panels) and D_cG = 330 μm²/s (Middle Panels). Increased diffusivity raised variability. Linearizing the fluxes Eqs 10–12 also increased variability. Although the deterministic simulations suggested that the effect of linearization on I(t_peak) variability was reduced with higher diffusivity (Table 1), computation of CV shows instead, that the relative error was substantial at all time points including t_peak.

The approximation [cGMP] ≈ [cGMP]_dark in the last term on the right-hand side of Eq 9, overestimates the variability of the SPR, by introducing stronger sources of randomness in the cGMP cascade. This term is the mechanism by which E* passes along its variability to cGMP. If [cGMP] is locally depleted, near the random path t → x(t) of R*, the random contribution of [E*]_σ is mitigated by the local smallness of [cGMP], through the product [E*]_σ[cGMP] in Eq 9, lowering the CV of the photocurrent. If, however, [cGMP] were to be uniformly close to its dark value [cGMP]_dark on the whole ROS and kept constant, the randomness of E* passes along to cGMP, and hence to the current suppression, linearly. Such an effect is local in nature, hence we refer to it as the local depletion of cGMP. The variability transferred from E* to cGMP is highest along the Brownian path of R*, where it is immediately reduced by a local sharp drop of cGMP (the coefficient of [E*]_σ in Eq 9), and progressively mitigated as [cGMP] migrates, and [E*]_σ terminates. Fig 5 shows that irrespective of the values of D_cG, approximating [cGMP] ≈ [cGMP]_dark in the variability source term [E*]_σ[cGMP] in Eq 9, (NLD) yields a larger t → CV(I(t)), particularly near the peak of the SPR, when cGMP depletion would be maximal. Finally if both local depletion and nonlinear effects are removed (LIN+LND), then the variability of the current suppression, after an initial delay due to diffusion eventually surpasses that of E*.

The deterministic simulations indicated that linearization increased variability as assessed by I_area. Fig 4 extends the analysis of CV to t → CV(I_int(t) dynamics for the FSR, TWS and GWS models. In all cases variability increased with decreased space resolution. In particular the GWS model passes along the variability of E* one-to-one to that of the response. In all cases by keeping [cGMP] = [cGMP]_dark in Eq 9 and linearizing the nonlinear fluxes of Eqs 10–12, or both, increased the variability of the SPR. With the combined effect, CV[I_int(t)] matched CV[E*(t)] after a delay and since cGMP diffusion was no longer relevant, it was the same for both values of D_cG and for all models. While mitigated by an increasing degree of space resolution from FSR to GWS, localization and nonlinearity then appear as stabilizing factors to the randomness of the response.

Simulations of ([15]), with a TWS model, a deterministic single-step deactivation, and a relatively high longitudinal diffusion coefficient that favors rapid equilibration (D_cG ≈ 40 μm²/s), give a small relative [cGMP] suppression and z-profiles comparable to Fig 3 Upper panels. It is inferred in ([15]) that because of the smallness of this [cGMP] drop, the rate equations can be linearized, the cGMP drop passes one-to-one to current suppression, and that, while no variability simulations were presented, the local cGMP depletion does not contribute to variability suppression. The simulations and arguments presented here disprove these conclusions (Fig 4). Such effects evidenced by the global/integral functional t → I_int(t) persist when tracing the CV of the local in time current drop t → I(t), as reported in Fig 5. In all cases the variability suppression seems to act right from initial times and persists up to its own asymptotic limit. The use of these time-dependent functionals and their corresponding time-dependent variability is one of the novel points of this investigation with respect to those of ([12, 13]).

Salamander dynamics: Linearization and variability

SPR kinetics are slower in salamander than in mouse in part because ν_RG activity is lower (≈ 185/s versus ≈ 330/s for mouse), R* deactivation is 5 times slower (τ_av ≈ 0.4 s versus τ_av ≈ 0.08 s for mouse) and cGMP has to travel farther from the activation site to the outer shell plasma membrane (radius R ≈ 5.5 μm versus R ≈ 0.7 μm for mouse). The relative [cGMP] suppression at the activation site ρ = 0, is comparable to mouse and it remains large for longer times (Fig 2). Such a dramatic local drop is recorded by the FSR model and not detected by the TWS and GWS models. As expected, the cGMP suppression at the rim is much smaller in salamander than in mouse, with the TWS model overestimating it. The z-profiles of G ( R , θ , z , t peak ), i.e., the S-integrated relative cGMP suppression, are much smaller than in mouse (Fig 3), suggesting a linearization of the form Eq 17 might be valid for salamander. Table 3 computes the variability of the functionals I(t_peak), and I_area from the nonlinear relation Eq 2 and their counterparts I_lin(t_peak) and I_area;lin from the linearized Eq 17, using the FSR, TWS and GWS models, for deterministic simulations with R* lifetimes τ_av and 2τ_av. The functionals E max *, E area * and G S ( t peak ) are independent of the linearization. The distortions introduced by linearization were considerably smaller in salamander than in mouse. The relative errors between I(t_peak) and I_lin(t_peak) for the same deterministic runs did not exceed 3.11% for the FSR model and 5% for the TWS model, although the cumulative error over time with I_area was larger, of the order of 11%. Similarly, the effect of decreasing spatial resolution was somewhat less pronounced than in mouse. Thus in random trials with randomly long- and short-lived R*, the errors due to linearization were not expected to increase the coefficient of variation of the current suppression as much as in mouse. Simulations with stochastic R* shutoff substantiate these predictions (lower panels of Fig 4). For the FSR model, linearization of the J_cG-cGMP relation Eq 2 increases the CV of the current suppression only slightly, whereas the increase is more significant for the TWS model. Roughly speaking, in the TWS model, the variability of cGMP, inherited from E*, is passed along to I_int(t) with minor suppression. The variability of E* however is significantly suppressed as it is transmitted to cGMP because of the depletion term [E*]_σ[cGMP] in Eq 9. If [cGMP] were to be approximated to [cGMP]_dark, the coefficient of [E*]_σ would be relatively large and constant, and hence the randomness of E* would be transmitted one-to-one to cGMP, and passed to I_int(t) essentially unchanged. The sharp local drop of cGMP (Fig 2), makes the coefficient of [E*]_σ relatively small, resulting in a lower transmission of variability from E* to cGMP. The latter is further reduced by its migration by diffusion to the outer shell, and ultimately transmitted to the current drop. Hence, this local depletion, present in mouse and salamander, emerges as a key variability suppression mechanism.

The effects of Ca²⁺ feedback on variability

Molecules of cGMP located distant from E*and those generated by GC during the course of the SPR roughly speaking, act as first responders that rush by diffusion to the activation site to replenish cGMP being randomly depleted of by E*. In doing so, they dissipate “globally” the randomness of E* on the activated disc, dampening the transfer of randomness to the nearest CNG channels at the outer shell by which variability in the current is generated. With Ca²⁺ feedback intact, GC is stimulated close to the CNG channels that are experiencing the greatest light-induced fall in [cGMP], to further dampen “locally” the transfer of randomness. With Ca²⁺ clamped, the GC mediated production of cGMP remains constant; local dampening is absent, leaving global dampening to act alone, thereby elevating the variability of the current. Experimental results of ([16, 33]) suggest that mutant mouse rods lacking GCAPs (GCAP⁻/⁻) essentially behave as if Ca²⁺ were clamped. CV[I(t_peak)] increased to 0.42 from a baseline of 0.34 for WT ([16]). Simulations of ([16]) with a TWS model, using the same stochastic deactivation mechanism as in ([12, 13]), report an increase of CV[I(t_peak)] to 0.38 for GCAP⁻/⁻ rod vs 0.32 for WT. These interpretations were extrapolated from Table F of Fig 6 of ([16]). As a one time point value, CV[I(t_peak)] does not capture the variability over the time course of the SPR. According to an FSR model, clamping Ca²⁺ appears to have a negligible effect on CV of the cumulative/integral functional I_area ([12, 13]). The variability of I_area is considerably higher than the variability of I(t_peak) ([11]), but as indicated above, the former compiles the fluctuations of I(t) over the entire time course of the response. We repeated the simulations of ([12, 13]) for several values of D_cG, using FSR, TWS, and GWS models applied to salamander as well as mouse rods to see whether these results hold across species. We mapped the CV for several functionals including I(t) over the duration of the SPR and extended simulations to salamander.

The simulations are reported in Tables 5 and 6 and Fig 6. In all models and irrespective of D_cG, clamping Ca²⁺ in Eq 10, i.e., removing the Ca²⁺ feedback, had no effect on the CV(I_area), confirming the results of ([12, 13]). CV of the integral current suppression I_int(t) for Ca²⁺-clamped was initially slightly lower than the corresponding one for WT, for later times was slightly higher, and asymptotically it was the same. Since I_area = lim_t→∞ I_int(t) integrates the variability over the entire time course of the SPR, the lower variability at short, initial times was compensated by a slightly higher one at longer times, thereby explaining why the CV(I_area) did not detect a rise in variability due to clamping Ca²⁺. These conclusions are made possible by tracing the variability of the response by means of the time-dependent functional I_int(t) as opposed to the time-global functional I_area as in ([12, 13]).

Early during the rising phase of the SPR, CV of the “pointwise” current suppression I(t) for Ca²⁺-clamped lay below the corresponding CV for WT, irrespective of diffusion coefficient and across species (Fig 6, Upper panels), seemingly at odds with the values of Tables 5 and 6. The reason is that the SPR for WT and Ca²⁺-clamped peak at different times. In recordings of “Ca²⁺-clamped” GCAP⁻/⁻ mouse rods, the SPR reaches its maximum at a t_peak 2-3 fold later than in WT, and rises to an amplitude ≈ 4 fold higher than that of WT ([33–35]). As a consequence the CV of I(t_peak) for WT and Ca²⁺-clamped are comparable only in their own intrinsic time scales τ = t/t_peak.

Fig 6 Lower Panels follow the dynamics of the CV of the pointwise functionals τ → I_int(τ) and τ → I(τ), along the relative dimensionless time for mouse [Lower Panels A (D_cG = 120 μm²/s) and B (D_cG = 330 μm²/s)], and salamander [Lower Panel C]. A “crossing of variability curves” while invisible for CV[I_int(t)] in absolute times (Fig 6, Upper Panels), was revealed in relative times (Fig 6, Lower Panels). Precisely, clamping Ca²⁺ augmented the variability of the SPR with respect to WT, for a time course up to roughly 2 fold longer than its peak time for mouse, and roughly 1.3 fold longer than its time to peak for salamander. The dynamics of CV[I_int(⋅)] for Ca²⁺-clamped remained above that of WT at all relative intrinsic times. For this functional crossing of the variability curves disappeared on the relative time scale.

The “crossing of variability curves” is an epiphenomenon arising as an artifact of the choice of functionals by which one measures variability. In § we indicated that the CV of I(t) might be ill defined since the numerator of Eq 6 is the standard deviation of j_tot, which, as such, needs to be normalized by the statistical mean of J_tot(t) and not that of I(t). In formula Eq 6 the denominator is the statistical mean of the current drop. Now by the Ca²⁺-clamped recordings of ([33–35]) the SPR amplitude is ≈ 4 fold higher than that of WT. Assuming a comparable difference continues to hold along the time course of the SPR, the denominator in Eq 6 is considerably larger than the corresponding one for WT. This accounts for the drastic difference between the t → CV[I(t)] curve for WT and Ca²⁺-clamped. It also suggests measuring variability by the J_tot(t) functional. Fig 7 traces the dynamics of t → CV[j_tot(t)] as defined in Eq 7. By this measure the variability for Ca²⁺-clamped was at all times, absolute or relative, higher than that of WT. Thus the nonlinear and non-constant dependence on [Ca²⁺] of cyclase activity that leads to a light-induced up-regulation of cGMP synthesis as appearing in Eq 10, acted as a variability suppressor of the SPR.

It is worth stressing that the integral defining j_tot(t) in Eq 3 is extended on the outer shell S, whereas variability originates elsewhere. For long times, channels on S near the activated level z_* keep closing and only relatively few are permitted to reopen due to lack of cGMP production by GC. Thus for long times, under Ca²⁺ clamped conditions, [cGMP] is small at S, and the randomly long lived R* generate smaller random fluctuations of cGMP on S, and hence by formula Eqs 2–5 smaller current fluctuations. The reduction in current is made more dramatic by the nonlinear relation Eq 2. Indeed as [cGMP] → 0 on S, the corresponding current J_cG goes to zero much faster, precisely as [ cGMP ] m cG with m_cG ≈ 4. However, again by the limited cGMP production by GC, the system returns only very slowly to its steady state and hence the statistical average of j_tot(t) is not yet close to j_dark. These information in formula Eq 6 suggest a lowering of CV[I(t)] for times t > 2t_peak (Fig 6, Lower Panels).

In the WT situation, closed channels on S are continuously reopened by the input of cyclase mediated cGMP. As a consequence, randomly long lived R* cause continued random generation of cGMP on the outer shell and hence random current and current drop. However again by the Ca²⁺ stimulated continued cGMP production by GC, the system rapidly equilibrates and after a time t > 2t_peak the statistical average of j_tot(t) is close to its dark value j_dark, generating, by formula Eq 6, a large CV.

For this reason the curves of t → CV[I(t)] for WT and clamped Ca²⁺ in Fig 6 cross invert their variabilities at times larger than ≈ 1.5t_peak. Fig 7 tracks the current and hence cGMP, as a direct consequence of the relation Eq 2 and shows that after roughly the same time threshold of ≈ 1.5t_peak, the variability of j_tot(t) decreases for clamped Ca²⁺ much faster than for WT. These results underscore that the notion of “variability” is not an absolute one and depends on the choice of functionals used to measure it, each extracting different information depending on time and location of the response.

Conclusion

Animals have two different types of photoreceptors, ciliary (e.g., vertebrates) and rhabdomeric (e.g., insects) ([36, 37]). Some animals (e.g., marine rag-worm Platynereis) have both types ([38]). The phototransduction cascades in these two types of photoreceptors are quite different, yet one thing they share is multiple layers of photopigment-containing membranes. It is widely believed that this type of organization mostly serves to increase photon catch; the chance that the photon will be absorbed is greatly increased by having the light pass through many photopigment-containing membranes. In rhabdomeric and in ciliary cone photoreceptors, phototransduction occurs within the confines of restricted compartments to prevent second messengers from getting diluted in a large volume of cytoplasm. However, the design changed for ciliary rods for which there are multiple layers of cytoplasm partially separated by membranous discs containing R. Our modeling suggests a reason for this geometry.

The loose segmentation of the cytoplasm by multiple photopigment-containing membranes ensures that light-evoked changes in the concentrations of second messengers retain some locality. This appears to serve an important purpose. Change produced by the activation of just a few molecules of signal transducers and effectors in response to light capture by a single photopigment molecule are “diluted slightly” in a larger volume of photoreceptor cytoplasm. The limited dissipation of these changes suppresses the biochemical effect of inevitable random variation of the number of the effectors activated by a single photopigment, thereby reducing the variability of SPR.

Calcium clamp experiments show that if the calcium concentration is kept constant, the resulting variability at early time points is higher than in WT rods. Importantly, the variability is lower before the peak response is reached, i.e., at the time when bipolar cells respond to rod activation ([39]). Calcium clamping decreases the variability at later time points (Fig 6, Lower Panels). However, this late effect appears to be irrelevant for fully dark adapted bipolar cells, and therefore would not contribute to the reliability of vision. In summary, strict localization of both the cGMP response to a photon and calcium feedback, which is ensured by complex photoreceptor geometry with slivers of cytoplasm separated by discs forming diffusion barriers, serves to suppress SPR variability, which is inevitably generated by the randomly activated and deactivated biochemical cascade. Thus, by using photoreceptors with multi-disc structure nature apparently devised a way to overcome the limitations of biochemical reactions. Rods generate the response to single photons that is significantly less variable, and therefore more reliable, than pure biochemistry would allow. This feature is captured by FSR model, but missed by TWS and GWS models of phototransduction. For more information, see S1 Appendix.

Supporting information

S1 Appendix [a]
Collected appendices A-D.