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

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, Ca2+ 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 Ca2+ 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.


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 PLOS  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][8][9][10][11][12][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 2+ in the cytosol of the ROS, and (iii) the nonlinear relations linking ionic current to [cGMP] and [Ca 2+ ]. 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 2+ feedback was evaluated by "clamping" [Ca 2+ ] at the dark level. The suppressive effects of locality, nonlinearity and Ca 2+ 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 2+ exchange and the cGMP-gated channels, respectively, are given by In Eq 1, j sat ex is the saturated exchange current (as [Ca 2+ ] ! 1), K ex is the Ca 2+ concentration at which the exchange rate is half maximal, and S rod is the surface area of the lateral boundary of the ROS. In Eq 2, j max cG is the maximal cGMP-current (as [cGMP] ! 1), m cG is the Hill exponent, and K cG is the half-maximal channel opening concentration of cGMP. These formulae are "local" as they provide the current densities in terms of space-time values of [ (x, y, t) at (x, y) 2 S at time t, and the total current j tot , across the whole lateral boundary of the ROS, are where dS is the surface measure on S. The experimentally measured electrophysiological response, is the current suppression relative to its dark value, i.e., 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 first is the total relative charge suppression up to time t, the second is the total charge suppression over the entire time course of the SPR, the third is the current suppression I(t) when the response is maximal, i.e., at time t peak , and the last one t ! j tot (t), as defined by Eq 3, is the dynamic of the total actual current across the outer shell of the ROS along the time course of the process. The CVs for these functionals, with the exception of the last one, have appeared in the literature as a measure of SPR variability ( [7-13, 13, 15-20]). The integral quantities I int (t) and I area are regarded as suitable because "the area captures fluctuations occurring at any time during the response, and thus provides a good measure of the total extent of response fluctuations" ([9]). Pointwise fluctuations are tracked by I(t) and hence I(t peak ) ( [11,15,16]). Which of these best measures the randomness of the SPR, remains elusive. 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 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi R X ðj tot ðt; oÞ À While starting with the current drop I(t), the numerator of this fraction is the standard deviation of j tot (t) and not its relative drop. Normalizing the numerator by the probabilistic mean of j tot (t), gives the CV of j tot (t), i.e., CV j tot ðtÞ � � ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi R X ðj tot ðt; oÞ À Combining these formulae yields a relation between CV[I(t)] and CV[j tot (t)], of the form As t ! 1, the system returns to its dark adapted steady state. Since the models make no provisions for "dark noise" in the phototransduction system there are no further fluctuations, and variability reduces to zero. Indeed, for CV[j tot (t)] as defined by Eq 7, the denominator remains "close" to its dark value whereas the numerator goes to zero, so that the corresponding CV then goes to zero as t ! 1. But for CV[I(t)] as defined by Eq 6, both numerator and denominator go to zero, yielding a non-zero asymptotic value for the corresponding CV (for example, see §, §). 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 2+ clamped virtual experiments.

Diffusion of cGMP and Ca 2+ in the cytoplasm
Following ( [11,23,24]), in the cytoplasm [cGMP] satisfies the diffusion equation Here D cG is the diffusion coefficient of cGMP in the cytosol and D � x denotes the Laplacian with respect to the transversal variable � x ¼ ðx; yÞ only. These are diffusion processes, parametrized with z 2 (0, H), taking place on the homogenized horizontal layers of the ROS cytoplasm. Activation occurs at the level z = z� and d z � is the Dirac mass in z concentrated at z = z�. The coefficient k � s;hyd is the surface hydrolysis rate (in μm 3 s −1 /#) of cGMP by the surface density of [E � ] σ . The parameter ε o is the thickness of the discs and νε o is the width of the interdiscal spaces. The term F is given by The first term in round brackets in the definition of F represents the production of cGMP by GC, which is located on the faces of the discs. Here α max and α min are positive constants with dimensions μMs −1 . The last term in Eq 10 represents cGMP depletion due to hydrolysis of cGMP by PDE at a basal rate β dark . This process occurs at the faces of the discs and it involves the surface concentration of E through a surface hydrolysis rate k σ;hyd . The dynamic of Ca 2+ is described by the equation complemented by the flux condition of the lateral boundary S of the ROS Here D Ca is the diffusivity of Ca 2+ in the cytoplasm, n is the outward unit normal to S, the currents J ex and J cG are defined in Eqs 1 and 2, η (in nm) is a positive parameter and f Ca is the fraction of current carried by Ca 2+ through the CNG channels (see B of S1 Appendix, for the meaning and values of these parameters). 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 −2 ) 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 in D, complemented by homogeneous initial data, and no-flux boundary conditions on the boundary of D. Here k T � E is the coupling coefficient from T � to E � . The constant k E is the rate of deactivation of T � within the T � -E � complex. The constants ν j are the rate of activation of T � by R � j through a successful encounter at time t 2 (t j−1 , t j ]. The times s j = (t j − t j−1 ) are exponentially distributed random sojourn times of R � in its jth phosphorylation state. The method is introduced in detail in ( [12]), along with the corresponding parameters, which are reported in B.1.1 of S1 Appendix, and in B.2.1 of S1 Appendix for salamander.
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][12][13][14][23][24][25][26][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., The variability of [E � ] σ is computed by the variability of the functionals The first is the activity of E � up to time t, the second is the total activity of E � produced over the entire lifetime of the process and the last is the maximum of E � (t) at its peak time t � peak . While these functionals directly express the randomness of the activation/deactivation cascade, to our knowledge, they are not experimentally accessible in intact rods.
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 6 ¼ 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][12][13][14]) (see B of S1 Appendix). The cGMP diffusion coefficient for mouse was taken as D cG = 120 μm 2 /s, close to the reported experimental value D cG � 140 μm 2 /s ( [28]). In our mouse simulations we also tested a theoretical value, D cG = 330 μm 2 /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 2 /s and � 40μm 2 /s. Using D cG = 330 μm 2 /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 2 .

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 ! 1 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 Gðr; y; z; tÞ ¼ def 1 À ½cGMP�ðr; y; z; tÞ ½cGMP� dark ; Here (ρ, θ) are the polar coordinates on the activated disc, and z 2 (0, H) is the longitudinal variable along the axis of the ROS. The locality of E � and that of Ca 2+ feedback onto GC produce strong cGMP concentration gradients, both across the activated disc, and along the axis of the ROS. Fig 2 illustrates the time dynamics of Gðr; y; z; tÞ at the activation site (ρ = 0) and at the rim (ρ = R) of the activated disc as described by the FSR, TWS and GWS models. While mitigated by large diffusivities, which tend to equilibrate the [cGMP], the gradient between activation site and rim persists at all diffusion coefficients 0 < D cG < 1, and across species. The TWS and GWS models do not distinguish between center (ρ = 0) of the activated disc and its rim (ρ = R), thereby missing the dramatic drop in cGMP near the activation site. The magnitude of change is smaller with the GWS model because the cGMP reduction is distributed along the entire length of the ROS. The average axial z-profile of [cGMP](ρ, θ, z, t peak ) depletion at the outer shell of the ROS at the peak of the SPR is shown in Fig 3. These profiles result from tracking locally, the movement of the second messengers within the ROS by the FSR model. The TWS produced a similar z-profile of [cGMP] drop that slightly overestimated the peak suppression. In marked contrast, the GWS model yielded a lower maximal [cGMP] depletion, spread uniformly along the ROS. Deterministic simulations were carried out for short and long lived R � with lifetime τ av and 2τ av , where τ av was computed as the sum of the   ] 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 2 s −1 ( [13]) and D cG = 330 μm 2 s −1 ( [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 2+ 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.

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 2+ ], 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  Table 2. Mouse: δ% is the relative error between current density suppression J(t peak ) at peak time t peak and its linearized version J lin (t peak ) for R � lifetimes τ av and 2τ av , using the FSR, TWS and GWS model. A: Current density J(t peak ) and J lin (t peak ) computed at the rim of the activated disc; B: Total current suppression integrated over the outer shell I(t peak ) and I lin (t peak ). 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 2 /s (Upper Panels) and D cG = 330 μm 2 /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 .

FSR
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 Table 3

. Salamander: First estimations of SPR variability and the errors introduced by linearization and loss of spatial resolution.
Deterministic simulations with R � lifetimes τ av and 2τ av , using the FSR, TWS and GWS models. Δ% is relative difference of the outputs for each pair of runs. Computation of E � area , and the maximum value E � max of E � (t) as defined by Eqs 14 and 15 for each of these two runs. Computation of I(t peak ) using the nonlinear relations of Eqs 1 and 2 and I lin (t peak ) using the linearization Eq 17 for the cGMP component of the current and the nonlinear relation Eq 2 for J ex . E � max , E � area and G S ðt peak Þ are independent of the linearization Eq 17.  Table 4. Salamander: δ% is the relative error between current density suppression J(t peak ) at peak time t peak and its linearized version J lin (t peak ) for R � average lifetimes τ av and 2τ av , using the FSR, TWS and GWS models. A: Current density J(t peak ) and J lin (t peak ) computed at the rim of disc containing R � ; B: Total current suppression integrated over the outer shell I(t peak ) and I lin (t peak ).  [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.

FSR TWS GWS
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 2 /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; y; 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 2+ 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 2+ 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 2+ 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 2+ 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  ([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 2+ in Eq 10, i.e., removing the Ca 2+ 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 2+ -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!1 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 2+ . 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 2+ -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. Locality and nonlinearity reduce variability in the SPR of retinal rods The reason is that the SPR for WT and Ca 2+ -clamped peak at different times. In recordings of "Ca 2+ -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][34][35]). As a consequence the CV of I(t peak ) for WT and Ca 2+ -clamped are comparable only in their own intrinsic time scales τ = t/t peak .  (Fig 6, Upper Panels), was revealed in relative times (Fig 6, Lower Panels). Precisely, clamping Ca 2+ 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 2+ -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 2+ -clamped recordings of ( [33][34][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 2+ -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 2+ -clamped was at all times, absolute or relative, higher than that of WT. Thus the nonlinear and non-constant dependence on [Ca 2+ ] 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 2+ 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 2+ 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 2+ 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 2+ 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 Increased SPR variability with Ca 2+ clamping, as assessed with the j tot (�) functional defined in Eq 3. Time dynamics of CV(j tot (�)) (current) are shown for WT (red traces) and Ca 2+ clamped (dashed blue traces) SPRs. In all cases the FSR model was in force. Panels A1, A2: mouse with D cG = 120 μm 2 /s. Panels B1, B2: mouse with D cG = 330 μm 2 /s. Panels C1, C2: salamander. Upper Panels A1, B1, C1: Time dynamics of t ! (j tot (t)) (current), in absolute times t in s. Lower Panels A2, B2, C2: Time dynamics of τ ! (j tot (τ)) (current), in relative times τ = t/t peak , so that τ = 1 corresponds, in each case, to the amplitude of the SPR.
https://doi.org/10.1371/journal.pone.0225948.g007 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. Collected appendices A-D. Appendices include the weak formulations for 2nd messengers and the G-protein cascade (A), parameters used for mouse and salamander models (B), a calibration for the choice of ν RG when D cG = 330 μm 2 /s (C), and a theoretical calculation linking volumic and longitudinal diffusivities (D). (PDF)