Analysis of dim-light responses in rod and cone photoreceptors with altered calcium kinetics

Rod and cone photoreceptors in the retina of vertebrates are the primary sensory neurons underlying vision. They convert light into an electrical current using a signal transduction pathway that depends on Ca\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{2+}$$\end{document}2+ feedback. It is known that manipulating the Ca\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{2+}$$\end{document}2+ kinetics affects the response shape and the photoreceptor sensitivity, but a precise quantification of these effects remains unclear. We have approached this task in mouse retina by combining numerical simulations with mathematical analysis. We consider a parsimonious phototransduction model that incorporates negative Ca\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{2+}$$\end{document}2+ feedback onto the synthesis of cyclic GMP, and fast buffering reactions to alter the Ca\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{2+}$$\end{document}2+ kinetics. We derive analytic results for the photoreceptor functioning in sufficiently dim light conditions depending on the photoreceptor type. We exploit these results to obtain conceptual and quantitative insight into how response waveform and amplitude depend on the underlying biophysical processes and the Ca\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{2+}$$\end{document}2+ feedback. With a low amount of buffering, the Ca\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{2+}$$\end{document}2+ concentration changes in proportion to the current, and responses to flashes of light are monophasic. With more buffering, the change in the Ca\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{2+}$$\end{document}2+ concentration becomes delayed with respect to the current, which gives rise to a damped oscillation and a biphasic waveform. This shows that biphasic responses are not necessarily a manifestation of slow buffering reactions. We obtain analytic approximations for the peak flash amplitude as a function of the light intensity, which shows how the photoreceptor sensitivity depends on the biophysical parameters. Finally, we study how changing the extracellular Ca\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{2+}$$\end{document}2+ concentration affects the response.


Introduction
Vision in most vertebrates starts with the absorption of light by rod and cone photoreceptors located in the retina of the eye (Ebrey and Koutalos 2001).Rods are highly sensitive to light and sustain vision under scotopic light conditions, whereas cones are much less sensitive compared to rods, but they adapt and maintain vision even under brightest illumination (Burns and Baylor 2001).The perception of light starts with the absorption of light by photopigments located the outer segment (OS) of rod and cone photoreceptors, which initiates a phototransduction cascade that leads to an electrical current response (for reviews see (Burns and Pugh 2010;Arshavsky and Burns 2012;Arshavsky et al. 2002;Pugh and Lamb 2000)).Many of the biochemical processes that participate in this signal transduction pathway are modulated by Ca 2+ feedback (Vinberg et al. 2018;Koch and Dell'Orco 2013;Nakatani et al. 2002;Korenbrot and Rebrik 2002;Rispoli 1998).Ca 2+ feedback not only enables photoreceptors to adapt their sensitivity to increasing light intensities in order to avoid early saturation (Yau and Nakatani 1985;Nakatani and Yau 1988;Matthews et al. 1988;Nakatani and Yau 1988;Matthews 1991;Koutalos and Yau 1996;Fain et al. 2001;Pugh et al. 1999), but this feedback also affects the light response in darkness (Klaus et al. 2021;Burns et al. 2002;Sakurai et al. 2011;Koutalos et al. 1995;Lagnado et al. 1992;Torre et al. 1986).For example, in most species, dark-adapted rods show monophasic responses to brief flashes of light (Hamer et al. 2003;Field and Rieke 2002;Pugh andLamb 2000, 1993).However, when the Ca 2+ kinetics are distorted and slowed down by the application of exogenous buffers, biphasic flash responses have been observed in amphibian rods (Torre et al. 1986;Lamb et al. 1986;Torre et al. 1986;Korenbrot and Miller 1989;Rieke and Baylor 1998), in primate and guinea pig rods (Field and Rieke 2002;Burns et al. 2002;Matthews 1991), and in mouse rods (Burns et al. 2002;Makino et al. 2004).For cones, biphasic responses have even been observed without the application of exogenous buffers (Korenbrot 2012;Holcman and Korenbrot 2005;Schneeweis and Schnapf 1999;Schnapf et al. 1990;Baylor et al. 1987;Schnapf et al. 1987;Matthews et al. 1990;Nunn et al. 1984).Recently it has been claimed that biphasic cone responses are not physiological but generated by experimental conditions that distorted the Ca 2+ dynamics (Cao et al. 2014).
The dynamics of the free Ca 2+ concentration in the OS depends on influx through cyclic GMP (cGMP) gated cyclic nucleotide-gated (CNG) channels, efflux via NCKX exchangers, and buffering reactions (Pugh and Lamb 2000).The impact of buffering has been modelled either using explicit equations for the interactions of free Ca 2+ with a combination of slow and fast buffers (Forti et al. 1989;Tamura et al. 1991;Hamer et al. 2003;Dell'Orco et al. 2009;Chen et al. 2010;Invergo et al. 2014;Astakhova et al. 2015), or by considering a simplified approach with an effective buffering capacity B ca (Sneyd and Tranchina 1989;Tranchina et al. 1991;Nikonov et al. 1998;Pugh and Lamb 2000;Caruso et al. 2005;Gross et al. 2012a, b;Reingruber et al. 2013;Lamb andKraft 2016, 2020) (in Korenbrot (2012) a Ca 2+ -dependent value for B ca has been assumed).Although the latter models are only valid for fast buffering kinetics, it turns our that they are nevertheless often sufficient to reproduce experimental data.
Since Ca 2+ feedback is ubiquitous in biology, the Ca 2+ kinetics and the effect of buffering have been studied both analytically and numerically for many biological systems, see for example (Falcke 2004;Wagner and Keizer 1994;Sneyd et al. 1995).In phototransduction, fitting procedures and numerical simulations have been mostly applied to adjust the Ca 2+ kinetics to experimental data and to investigate the impact of buffering (Forti et al. 1989;Tamura et al. 1991;Hamer and Tyler 1995;Nikonov et al. 1998;Moriondo and Rispoli 2003;Hamer et al. 2003;Korenbrot 2012)(in addition to numerical simulations, some first-order analysis has been performed in (Schnapf et al. 1990;Nikonov et al. 1998)).Although simulations show that changing the Ca 2+ kinetics affects the response dynamics and the photoreceptor sensitivity, without analytic results it remains difficult to obtain a precise conceptual and quantitative understanding.
In this work, we combined numerical simulations with mathematical analysis to systematically study how the Ca 2+ kinetics affects waveform and amplitude of light responses in mouse rod and cone photoreceptors.We started from a parsimonious phototransduction model that comprises the principal transduction processes that are known to be relevant for dark-adapted photoreceptors, and that we have previously calibrated with experimental data from mouse rods and cones (Reingruber et al. 2020;Abtout et al. 2021).Whereas in previous work we assumed that under physiological conditions the free Ca 2+ concentration changes in proportion to the current, we now consider the general case where the free Ca 2+ concentration depends on influx via CNG channels, efflux via exchangers, and buffering.To keep the analysis comprehensible, we on fast buffering reactions and we ignore the impact of slow buffering.This is in line with models that use a buffering capacity B ca .To obtain analytic results, we performed a linear-response analysis for light intensities that are sufficiently dim such that changes of the cGMP and Ca 2+ concentrations, and subsequently of the current, are small compared to their corresponding steady state values in darkness.Because cones are much less sensitive to light than rods (Ingram et al. 2016), this analysis is valid up to much brighter light intensities for cones compared to rods.Thus, the definition of dim light depends on the photoreceptor sensitivity and therefore differs between rods and cones.For cones, dim light comprises a range of light intensities that is much larger than the one for rods.We combined analytic results and numerical simulations to investigate the light response as a function of the Ca 2+ kinetics.With a low amount of buffering, we find that the Ca 2+ concentration changes in proportion to the current, and we recover our previous results from (Reingruber et al. 2020;Abtout et al. 2021).In the low buffering range the waveforms of brief flashes of light are monophasic.As the amount of buffering increases, the Ca 2+ dynamics becomes delayed with respect to the dynamics of cGMP and current, and biphasic waveforms emerge that contain a damped oscillation.A phase space analysis shows that the transition from monophasic to biphasic responses depends on the ratio between the rate μ ca that controls the the Ca 2+ kinetics, and the dark turnover rate of cGMP β d .Because β d changes among species, and between rods and cones (Reingruber et al. 2020;Pugh and Lamb 2000), the same Ca 2+ kinetics does not entail similar waveforms in rods and cones.We use our analytic results to dissect the contributions of the various biophysical processes to the waveform, and to identify the processes that limit the response recovery.We further derive an analytic approximation for the peak amplitude as a function of the flash intensity, which reveals how the photoreceptor sensitivity depends on the underlying biophysical parameters and the Ca 2+ kinetics.Finally, we investigate the response of cones to stimulations with longer steps of light, and we study the effect of changing the extracellular Ca 2+ concentration.

Phototransduction model
We start from the phototransduction model from (Abtout et al. 2021;Reingruber et al. 2020), and we consider here a more general Ca 2+ dynamics where the free Ca 2+ concentration does not necessarily change in proportion to the circulating current.In short, light activates the visual pigment R * with a rate that is proportional to the light intensity φ(t) times the collecting area κ.R * activates the G-protein transducin T * with a rate k act , and deactivates with rate μ rh .T * activates phosphodiesterase (PDE) P * with a rate k tr , and deactivates with rate μ tr .Because the deactivation of T * is linked to the activation of PDE, we have k tr = μ tr (to keep the analysis most general, we distinguish between k tr and μ tr ).P * deactivates with rate μ pde , and hydrolyses the cytosolic second messenger cyclic GMP (cGMP) with rate constant β sub .cGMP gates the opening of CNG channels, in the OS membrane.The rate by which the cGMP concentration c cg is synthesized by guanylate cyclase (GC) depends on the , where r α = α min α max .The GC activity is modulated via an intermediate step where Ca 2+ binds to GCAPs proteins (Mendez et al. 2001;Burns et al. 2002).Ca 2+ feedback not only modulates the cyclase activity, it also affects the deactivation rate of an activated photopigment via binding to recoverin, which inhibits the phosphorylation of the photopigment by rhodopsin kinase (Klenchin et al. 1995;Chen et al. 1995Chen et al. , 2012)).However, in this work we focus on the Ca 2+ feedback to the cyclase, which is the most effective feedback under dimlight conditions (Burns et al. 2002;Sakurai et al. 2011;Koutalos et al. 1995).We assume that the Ca 2+ -dependent modulation of photopigment deactivation occurs on a much slower time scale, which affects adaptation process but can be neglected in first approximation for flash responses.
The Ca 2+ content of the OS changes due to influx via CNG channels and efflux via electrogenic NCKX exchangers.We adopt the convention that inward currents are negative.The CNG current is I ch = I ch,max p ch , where p ch (c cg ) = is the exchanger saturation level, and I ex,sat is the saturating current at high c ca .The exchanger stoichiometry is such that the extrusion of a single Ca 2+ ion leads to the influx of a positive charge e + .Thus, the exchanger current is negative, and the Ca 2+ efflux is − I ex e + .Although the exchanger cooperativity n ex is found to be one (Pugh and Lamb 2000), we keep a general n ex for the analysis to see how this parameter affects the results.But for numerical evaluations we use n ex = 1 (Table 1).Ca 2+ regulates the activity of many proteins in the outer segment (Fain et al. 2001;Nakatani et al. 2002;Korenbrot and Rebrik 2002;Rispoli 1998), which contributes to Ca 2+ buffering.With N b buffer species b i (i = 1, . . .N b ), the Ca 2+ concentrations that are bound to these buffers are c ca,b i .With standard buffer kinetics (Keener and Sneyd 1998;Wagner and Keizer 1994), the equations for the total and buffered Ca 2+ concentrations are (1) c b i ,tot are total buffer concentrations, K b i are dissociation constants, μ b i are dissociation rates, F A = 9.65 × 10 −5 sp A μMμm 3 is the Faraday constant, and V os is the outer segment volume.We assume that the buffer kinetics are fast (the times μ −1 b i are small compared to the time scale for the change in the free Ca 2+ concentration due to the current), in which case the amount of buffered Ca 2+ becomes a function of the free (quasi steady-state approximation) (Keener and Sneyd 1998;Wagner and Keizer 1994).Equation 1 now simplifies to (2) The closed system of transduction equations is where β d is the basal cGMP hydrolysis rate in darkness.The total current is To study the dynamics with respect to the steady state in darkness, we use steady state values c ca,0 , c cg,0 and I 0 to introduce normalized variables and parameters (Sneyd and Tranchina 1989) f ch,ca I ex (c ca,0 ).We label normalized variables and parameters with a hat, and we define ĉcg = N b i=1 B b i .With the new variables P * = β sub P * , T * = β sub k tr μ pde T * , R * = β sub k act k tr μ pde μ tr R * , we obtain from Eq. 3 with the gain ξ = k act β sub k tr μ rh μ pde μ tr and The initial conditions in darkness with φ = 0 are R * = T * = P * = 0 and ĉcg = ĉca = α = pch = pex = 1.The current normalised by the dark current is We further introduce the normalised current î = 1 − Î , which is zero in darkness.
With Eq. 7 we get î is usually used in the literature for the light response.Although î is referred to as a current, one has to keep in mind that î is the current change with respect to the dark current.Î and î have complementary properties, for example, whereas Î decreases after the light is switched on, î increases.Before proceeding with the analysis, we add some remarks: 1) Because we study the light response in darkness, we defined the buffering capacities using the steady state Ca 2+ concentration in darkness.In presence of a background light, we would use the steady state Ca 2+ concentration corresponding to this background light.
2) By definition, the expression (1 is one in darkness, and can be neglected in first order.Hence, for low activation we have d dt ĉca ≈ μ ca pch − pex .3) Due to the normalisation, Eq. 5 depends only on kinetic parameters that govern the dynamics.Equation 5 therefore is convenient to study the response dynamics as a function of the light intensity since it contains less parameters than Eq. 3.For example, Eq. 5 shows that the activation rates k act , k tr and β sub affect the light response only as the product k act k tr β sub .4) Eq. 5 has been obtained assuming that the normalisation values c ca,0 and c cg,0 are the steady state concentrations in darkness.In Eq. 5 the steady state in darkness with P * = 0 is therefore always ĉcg = ĉca = 1.To modify these values, we have to add additional parameters to Eq. 5.For example, to investigate how changing the cyclase activity alters the steady state, we introduce

Linear response analysis
The equations for PDE activation in Eq. 5 are linear and can be solved explicitly.The solution is with the Green's function The double-headed arrows in Eq. 10 signify that the second and third terms in the equation are identical to the first term with the exchange μ tr and μ rh in the second term, and μ pde and μ rh in the third term.The non-linear equations for ĉcg and ĉca (Eq.5) cannot be solved analytically for general φ(t).However, analytic results can be derived for dim light where the changes of ĉcg and ĉca are small compared to the steady state values ĉcg = ĉca = 1.For the analysis we introduce the new variables y = − ln ĉcg and z = − ln ĉca , which are both zero in darkness.With y 1 and z 1, by linearizing Eq. 5 we obtain the first-order system of equations where νu = z and Note that α 0 ≥ 0. The solution of Eq. 11 with P * (t) from Eq. 9 is where are the eigenvalues of the matrix (the double-headed arrows have the same significance as in Eq. 10).
Experimentally, photoreceptors are often studied by stimulating them with flashes or steps of light.To connect our analysis to such data, we compute the photoreceptor response to a light-step with duration t and intensity φ.The number of isomerisations produced by this stimulus is R * 0 = κφ t.From Eq. 13 we obtain y = R * 0 ξ g y , u = R * 0 ξ g u and z = R * 0 ξνg u , with From Eq. 8 we obtain in first order for the current For t → 0 we have 1 With Eq. 14 we obtain that î is a sum of exponentials, and similar expressions apply for the other coefficients.For R * 0 = 1 we obtain the single-photon response (SPR).

Steady state with constant light intensity
With constant φ we have P * = ξκφ.The steady state of Eq. 11 is u = y and y = P * . With the steady state condition pch = pex we obtain for the current in dim light where we used Kch 1.For brighter light, the steady state is computed from Eq. 5 by numerically solving pch ( ĉcg ) − pex ( ĉca ) = 0 with ĉcg = β d α( ĉca ) β d + P * .

Emergence of damped oscillations
We investigate how the current îβ in Eq. 17 changes as a function of the the Ca + kinetics.For this we very the rate μ ca , respectively the parameter r ∼ μ ca β d .For (1 − r ) 2 − 4r ν α 0 > 0 the eigenvalues λ 1 and λ 2 are real and positive.In this range the origin y = z = 0 is a stable node and flash responses are monophasic (Strogatz 2018).For example, with fast Ca + kinetics (low amount of buffering such that r 1) we find λ 1 ≈ 1 + ν α 0 , λ 2 ≈ r − ν α 0 → ∞ and îβ ∼ e −β d (1+ν α 0 )t .In the opposite limit of very slow Ca + kinetics (high amount of buffering) the free Ca 2+ concentration remains constant.For r → 0 we have For (1 − r ) 2 − 4r ν α 0 < 0 the eigenvalues λ 1 and λ 2 are complex conjugate with non-vanishing real part λ re = 1+r 2 .In this range the origin is a stable spiral and îβ is a damped oscillation (biphasic response) (Strogatz 2018).Oscillations occur for r 1 < r < r 2 , where are the values where the origin changes between stable node and stable spiral.By writing λ 1 = λ re − iλ im and λ 2 = λ re + iλ im , with , we find The damping rate is , and the oscillation rate is ω = . Because c β 1 and c β 2 are complex conjugate, the amplitude a and the phase ϕ of the oscillation can be computed from a 2 e iϕ = c Since λ 1 and λ 2 are hyperbolic with non-vanishing real parts, the Hartman-Grobmann theorem states that the analysis of the linearized system in Eq. 11 also faithfully reflects the behaviour of the non-linear system in Eq. 5 locally (Strogatz 2018;Ricardo 2020).For a steady background light φ, Eq. 5 has a single fixed point (stable node or spiral), and no limit cycles and self-sustained oscillations exist.This ensures that the visual perception is driven by the light input.Given that Eq. 5 contains only negative feedback, the lack of limit cycles is not surprising (however, there are special conditions where a system with only negative feedback can exhibit Hopfbifurcations (Reidl et al. 2006)).

Changing the extracellular Ca 2+ concentration
We now investigate how the response changes when the extracellular Ca 2+ concentration is altered by a factor of r ca,ex with respect to the reference concentration that had been implicitly assumed for the previous computations.We assume that the CNG current carried by Ca 2+ changes in proportion to the extracellular Ca 2+ concentration.We checked that this is a valid assumption by using the more general Goldman-Hodgkin-Katz(GHK) equation (Keener and Sneyd 1998) to estimate the current change.We further assume that the CNG current not carried by Ca 2+ remains unaffected by changing extracellular Ca 2+ .Hence, we have Îch,ca = r ca,ex f ch,ca , where is the new fraction of the current that is carried by Ca 2+ .The total CNG and exchanger current is The equation for ĉca reads Thus, for r ca,ex = 1 we have that ĉca = ĉcg = 1 are not the steady state solutions in darkness.The new steady state values ĉca,d and ĉcg,d = α( ĉca,d ) are obtained by solving the equation r ca,ex pch ( α( ĉca,d )) − pex ( ĉca,d ) = 0.With pex = r ca,ex pch we obtain from Eq. 21 for the new dark current We can use the new steady state values ĉca,d c ca,0 , ĉcg,d c cg,0 and Îd I 0 to renormalize parameters.Together with the new value for f ch,ca , we again obtain equations like Eqs. 5 and 6 to study the light response.

Results
Equation 5 shows that the Ca 2+ kinetics are determined by the effective rate . Without buffering (B ca = 0) and parameters from Table 1 we compute μ ca ∼ 1.6 × 10 3 s −1 for a rod, and μ ca ∼ 9.4 × 10 3 s −1 for a cone.For large μ ca we have the quasi steady state approximation pch ≈ pex and Î ≈ pch ≈ pex .Moreover, with Kex 1 we have pex ≈ ĉca and ĉca ≈ Î .Thus, with fast kinetics we have that the Ca 2+ concentration changes in proportion to the current, and in this limit we recover the model and results from (Reingruber et al. 2020;Abtout et al. 2021).
With fast kinetics we observe only monophasic flash responses (Fig. 1A and C), whereas biphasic responses (damped oscillations) emerge as the Ca 2+ kinetics are Fig. 1 Simulations of flash responses for mouse rod and cone.The normalized current î from Eq. 8 is computed with Eq. 5 and parameters from Table 1 with a slow and a fast buffer species.A and B Family of flash responses with duration t = 5ms for a rod with fast and slow Ca 2+ kinetics μ ca , as indicated in the panels.The around 25 times higher buffering capacity B ca in B has been obtained by increasing both buffering concentrations by the same factor.The legend gives the number of isomerisations R * 0 produced by the flash.C-D Similar to A-B but for a cone slowed down (Fig. 1B and D).The simulations in Fig. 1 were performed with Eq. 5 and parameters from Table 1.We used one low and one high affinity buffer.The dissociation constant K b 1 = 3μM for the low affinity buffer corresponds to recoverin (Chen et al. 1995), whereas K b 2 = 0.14μM for the high affinity buffer corresponds to the GCAP proteins (Kawamura and Tachibanaki 2021).Since recoverin is the most abundant Ca 2+ buffer in the outer segment (Pugh and Lamb 2000; Kawamura and Tachibanaki 2021), it follows that under physiological conditions the buffering capacity of recoverin B b 1 is much larger than the buffering capacity of the GCPAs B b2 .However, to show that our results are not biased towards low or high affinity buffers, we use for simulations such that low and high affinity buffers contribute equally to the total buffering capacity B ca = B b 1 + B b 2 (analytic results only depend on B ca ).Thus, for the simulations we altered the value of μ ca by changing the buffering capacity B ca with the constraint that B ca = 0.5.Physiologically this corresponds to changing both buffer concentrations by the same factor.For example, with parameters from Table 1 we obtain for a rod μ ca = 50s −1 using B ca = 32, and for a cone we get μ ca = 200s −1 with B ca = 45.We use μ ca rather than B ca to characterize the Ca 2+ kinetics because μ ca can be directly compared to the other rate constants of the model.

Waveform and dynamics of flash responses
The current waveform is obtained by normalizing a flash response with its peak amplitude.The waveform therefore characterizes the current dynamics.The waveforms of the rod and cone simulations in Fig. 1 change only little up light intensities where the peak amplitude reaches values of the order 0.5 (Fig. 2, black solid lines).In this range, the analytic waveform w = î/ î peak computed with Eq. 16 faithfully agrees with the simulations for rods and cones with fast and slow Ca 2+ dynamics (Fig. 2, black lines versus red dashed line).
Next we used the analytic waveform to study how the the response changes as a function of μ ca .We further introduced waveforms for cGMP and Ca 2+ concentration, w cgmp = g y g y, peak and w ca = g u g u, peak , respectively.Because the exchanger current is small, the overall waveform is determined by the cGMP dynamics (Fig. 3C-D, dashed versus solid curves).We find that waveforms are monophasic with fast Ca 2+ kinetics (Fig. 3A-B, green curves).In this limit the Ca 2+ concentration changes in proportion to the current (Fig. 3C-D, green dotted versus solid curve).Biphasic waveforms (damped oscillations) emerge as the Ca 2+ kinetics are slowed down (Fig. 3A-B).In this case the Ca 2+ concentration becomes delayed with respect to the current (Fig. 3C-D, red dotted versus solid curve).
With Eq. 17 we decompose the current waveform into the contributions of the various biophysical processes, w = w rh +w tr +w pde +w β (Fig. 4).Such a decomposition cannot be obtained from simulations.The waveform decomposition shows that the recovery in a rod with fast Ca 2+ kinetics are limited by PDE decay, w pde ∼ e −μ pde t (Fig. 4A).The recovery in a cone with constant Ca 2+ concentration (μ ca = 0) is limited by β d , w β ∼ e −β d t (see the analysis after Eq. 17 for details).This is also the case for the recovery in a GCAPs −/− cone where the cyclase is not Ca 2+ dependent.Note that the Ca 2+ current is not constant in a GCAPs −/− mutant, contrary to the case with μ ca = 0.However, because the Ca 2+ current is small, flash responses in GCAPs −/− mutant mice vary only little when the Ca 2+ dynamics is distorted, contrary to the Wt case (Burns et al. 2002).Except for such special cases, the recovery of a flash response cannot be approximated by a single exponential decay function.For example, because μ pde and β d have similar values in a mouse rod, the recovery in a GCAPs −/− rod is determined by a sum of two exponentials (see Fig. 4D in Abtout et al. ( 2021)).As another example, with slow Ca 2+ dynamics the recovery is determined by a damped oscillation (Fig. 4C and D).

Peak amplitude of flash responses and photoreceptor sensitivity
The Ca 2+ kinetics affects not only the waveform of a light response, but also the peak amplitude (sensitivity).As the Ca 2+ kinetics becomes slower, the peak amplitude of the SPR increases by a factor up to 1.8 for a cone, and up to 2.3 for a rod (Fig. 5A).A slower Ca 2+ kinetics leads to a smaller reduction of the Ca 2+ concentration at time to  peak, a reduced negative cyclase feedback and a larger SPR.A larger SPR also implies a higher light sensitivity.Because the exchanger current is small, the condition with μ ca → 0 also characterizes a GCAPSs −/− photoreceptor.Thus, Fig. 5A shows that the difference of the SPR between a WT and GCAPSs −/− photoreceptor depends on the Ca 2+ kinetics, and the largest difference is attained when the Ca 2+ kinetics are fast.
In first approximation, the peak current increases linearly with the number of isomerisations R * 0 (Eq.16).The linear approximation is valid up to î peak ∼ 0.2 − 0.3 (Fig. 5B; dashed lines are the peak current computed with Eq. 16, diamonds are simulation results), which corresponds to few isomerisations for a rod, but few hundreds for a cone (for example, with fast Ca 2+ kinetics, the peak SPR amplitude in a rod is ∼ 0.07 and ∼ 4.7 × 10 −4 in a cone).The slight discrepancy in Fig. 5A between the simulation and linear approximation for a rod with slow Ca 2+ kinetics arises because of the large SPR amplitude î peak ∼ 0.16.
To obtain an approximation for the peak current that remains valid up to higher light intensities, we combine the linear results with saturation effects.With y = R * 0 ξ g y and z = R * 0 ξνg u we get ĉcg = e −y = e −R * 0 ξ g y , ĉca = e −z = e −R * 0 ξνg u , pch ≈ e −R * 0 n ch ξ g y , and pex ≈ e −R * 0 n ex ξνg u ≈ e −R * 0 n ch ξ g u , where we used Kch 1 and Kex 1.With Eq. 8 we get î With fast Ca 2+ kinetics we have g y ≈ g u (the solution of Eq. 11 Fig. 3 Flash waveforms for different Ca 2+ kinetics.A Comparison of the analytic current waveform î/ î peak computed with Eq. 16 for a rod and for various values of μ ca as indicated in the legend.B Similar to A but for a cone.C Comparison of the analytic rod waveforms for current (solid line), cGMP concentration (dashed lines) and Ca 2+ concentration (dotted lines) for μ ca = 50s −1 (green color) and μ ca = 2s −1 (red color).The current waveform is î/ î peak , the cGMP waveform is g y /g y, peak , and the Ca 2+ waveform is g u /g u, peak (see Eq. 16).D Similar to (C) but for a cone with μ ca = 200s −1 (green color) and μ ca = 10s −1 (red color) (color figure online) for μ ca → ∞ is y = u) and î peak ≈ 1 − e −R * 0 ξ n ch g y, peak .With slow kinetics the Ca 2+ concentration and the exchanger current at peak time are only little affected (g u (t peak ) ≈ 0) and î peak ≈ 2 f ch,ca +2 1 − e −R * 0 ξ n ch g y, peak .These expressions for î peak for fast and slow Ca 2+ kinetics are now in good agreement with simulation results for rods and cones up to saturating flashes (Fig. 5B, diamonds versus solid lines).The numerical values for g y, peak used in Fig. 5B are computed with Eq. 15 and parameters from Table 1: for μ ca = 2s −1 we have g y, peak ≈ 0.32 for a rod and g y, peak ≈ 0.52 for a cone; for a rod with μ ca = 50s −1 we have g y, peak ≈ 0.17, and for a cone with μ ca = 200s −1 we have g y, peak ≈ 0.28.Since g y, peak is not much different between rods and cones, the higher sensitivity of a rod to a pigment activation is due to the larger transduction gain ξ .With parameters from Table 1 we find ξ rod ξ cone ≈ 250, which corresponds to the sensitivity gap between rods and cones in Fig. 5B.The sensitivity of a photoreceptor increases with slower Ca 2+ kinetics because of reduced negative Ca 2+ feedback at time to peak, in agreement with Fig. 5A.With slow Ca 2+ kinetics the exchanger current does not contribute at peak time, and flash responses Fig. 4 Decomposition of the current waveform of flash responses.The decomposition is obtained with Eq. 17.A The recovery in a rod with fast Ca 2+ dynamics is limited by PDE decay, w pde ∼ e −μ pde t .B The recovery in a cone with constant Ca 2+ concentration is limited by β d , w β ∼ e −β d t .This is also true for in GCAPs −/− mutant cones where the Ca 2+ feedback to the cyclase is genetically removed.C-D The recovery with slow Ca 2+ kinetics are determined by a damped oscillation saturate at 2 f ch,ca +2 (also the saturating flash responses in Fig. 1B,D that do not reach the maximal value of one).This effect is more pronounced in a cone due to a larger exchanger current and higher value f ch,ca .For intermediate Ca 2+ kinetics we do not have and analytic approximation; in this case one has to estimate the peak time t peak , and then compute
If ρ is between ρ 1 = 0.08 and ρ 2 = 16.9, the rates β 1 and β 2 are complex conjugate and îβ is a damped oscillation with damping rate β damp ≈ β d 1+ρ 2 (Fig. 6A blue line) and oscillation rate ω ≈ β d 4 is attained for ρ ≈ ρ 2 2 (Fig. 6B).Although the photoresponse contains a damped oscillation over a wide range of μ ca values, the oscillation can be observed during the recovery phase only if the damping rate β damp is not too large.For example, with ρ 2 = 16.9 we compute that oscillations are present for μ ca ρ 2 β d ≈ 69 s −1 in a rod and μ ca 186 s −1 in a cone.However, the rod waveform clearly exhibits an oscillation only for μ ca 10s −1 (Fig. 3A), and the cone waveforms for μ ca 50s −1 (Fig. 3B).It is difficult to give an estimate for the value of μ ca that characterizes when oscillations start to affect the recovery phase, since this depends on the damping rate, the initial oscillation amplitude, and other rate constants that affect the recovery.

Effect of changing the extracellular Ca 2+ concentration
We investigated the effect of changing the extracellular Ca 2+ concentration by a factor of r ca,ex for a cone.Similar conclusions are obtained for a rod (not shown).In the first (the inset magnifies the region for small μ ca ).The dashed black line shows the asymptotic value  1 Fig. 7 Effect of changing the extracellular Ca 2+ concentration for a cone.A Changes of the dark steady state concentrations of Ca 2+ (black curve) and cGMP (red curve), the dark current (blue curve) and the fraction of the CNG current that is carried by Ca 2+ (green curve) when the extracellular Ca 2+ concentration is altered by a factor of r ca,ex .ĉca,d has been computed numerically by solving the steady state of Eq. 22; ĉcg,d = α( ĉca,d ); Îd is computed with Eq. 23; fch,ca = is computed with Eq. 20.B Comparison of current waveforms computed with Eq. 16 for r ca,ex = 0.1 (solid lines) and r ca,ex = 1 (dashed lines) and for various values of μ ca as indicated in the legend (the dashed curves are as in Fig. 3B; dashed and solid curve overlay for μ ca = 0s −1 ) (color figure online) place, changing the extracellular Ca 2+ concentration alters the fraction of the CNG current that is carried by Ca 2+ (Fig. 7A, green curve), which changes the dark steady state concentrations of Ca 2+ and cGMP (Fig. 7A, black and red curve), and the dark current (Fig. 7A, blue curve).For example, with parameters from Table 1 we find that when the extracellular Ca 2+ concentration is reduced by tenfold (r ca,ex = 0.1), the fraction of the CNG current carried by Ca 2+ and the dark steady state concentration of Ca 2+ are reduced by a factor of around 7.3 and 2.5, respectively.The lower Ca 2+ concentration activates the cyclase, which results in a cGMP concentration that is elevated by a factor of around 1.9, and an dark current that is increased by a factor of around 2.9.
To estimate how the Ca 2+ kinetics governed by μ ca (Eq.6) is affected by extracellular Ca 2+ , we note that the ratio I ex c ca is independent of the Ca 2+ concentration for the physiological values n ex = 1 and Kex 1.Thus, the impact on the Ca 2+ kinetics depends on the change of the buffering capacity B ca .If the low affinity buffer recoverin is prevalent, B ca is not much altered by reducing extracellular Ca 2+ .In contrast, with high affinity buffers the value of B ca increases.For example, for a buffer with dissociation constant K b 2 = 0.14μM, a reduction of the dark-adapted Ca 2+ concentration from 0.3μM to 0.15μM due increases the buffering capacity B b 2 by a factor of 2.3.Assuming that low and high affinity buffers contributed equally with c ca = 0.3μM, it follows that B ca increases by a factor of (1 + 2.3)/2 ∼ 1.65 if the dark-adapted Ca 2+ concentration drops to c ca = 0.15μM, which reduces the rate μ ca by a factor of around 0.6.
Changing the extracellular Ca 2+ concentration strongly affects dimensional current amplitude (measured in units of p A) because the latter is proportional to the dark current.This scaling effect can be removed by normalizing responses with the corresponding dark current.The normalized flash response in Eq. 16 changes only moderately, which is mostly to the modified cyclase feedback α 0 .For example, with parameters from Table 1 we compute α 0 = 1.05 for r ca,ex = 1, and α 0 = 0.33 for r ca,ex = 0.1.The reduced cyclase feedback α 0 increases the normalized SPR amplitude by a factor between 1 and 1.5 (not shown) depending on the Ca 2+ kinetics (the strongest effect is obtain with fast Ca 2+ kinetics).In contrast, the dimensional current is additionally scaled by a factor around 2.9 due to the change in the dark current.
The reduced cyclase feedback α 0 = 0.33 with r ca,ex = 0.1 also affects the current waveform, leading to a longer time to peak, reduced oscillations and a slower recovery (Fig. 7B, solid versus dashed curves).The waveforms in Fig. 7B for r ca,ex = 1 and r ca,ex = 0.1 have been computed with the values of μ ca specified in the figure legend.Thus, if B ca and μ ca are also changed by reducing extracellular Ca 2+ , one has to compare the waveforms with the corresponding values of μ ca .

Responses to steps of light
Finally, we considered the photoreceptor response to steps of light with intensity φ and duration t for a cone.Similar results are obtained for a rod (not shown).For sufficiently long step duration t, the initial rising phase of the response is followed by an intermediate steady state phase (plateau phase), and a recovery phase after the light is switched off.Because the plateau current depends on φ but is independent of μ ca , we define the waveform by normalizing a step response with its plateau current.With the steady state solution g y,ss = g u,ss = 1 β d (1+ν α 0 ) t , the dim-light current waveform of a step response is given by The dynamics during initial and recovery phase depend on μ ca , and oscillations emerge as the Ca 2+ dynamics is slowed down (Fig. 8A).Whereas the flash waveform differs between initial and recovery phase (Fig. 3A-B), there is a symmetry between these phases for step responses (Fig. 8A).Indeed, with Eq. 11 we find î( t +t) = îss − î(t), where îss is the intermediate steady state current.
With the steady state relation î = 1 − pch and the first order result y = R * 0 ξ g y we get for cGMP ĉcg,ss = e −R * 0 ξ g y,ss , and for the current ( Kch 1) We find that Eq. 25 well approximates the results obtained by numerically solving the steady state of Eq. 8 up to saturating light intensities (Fig. 8B, blue versus black curve).For example, with parameters from Table 1 we compute 1.9 × 10 6 photon μm 2 s for a cone, and 132 photon μm 2 s for a rod.Equation 25shows that the steady state current in a GCAPs −/− photoreceptor with α 0 = 0 is by a factor of 1 + ñch α 0 ≈ 4 larger compared to WT.This difference is much larger than a value around 1.8 found for the difference in the SPR amplitude (Fig. 5A, blue curve).
Finally, by combining Eq. 24 with Eq. 25 we obtain for a step response Equation 26 is in good agreement with simulation results up to almost half saturating responses for fast Ca 2+ kinetics with μ ca = 100s −1 (Fig. 8C, dashed versus solid curves), as well as slow kinetics with μ ca = 20s −1 (Fig. 8D, dashed versus solid lines).

Discussion
The signal transduction pathway in the outer segment of rod and cone photoreceptors consists of a series of biophysical processes that transform light into an electrical current.Many of these processes are modulated by Ca 2+ feedback, which affects response dynamics and photoreceptor sensitivity.To obtain conceptual and quantitative insight that goes beyond numerical simulations, we studied the light response with a parsimonious model that allows for a comprehensive mathematical analysis.The model includes the principal transduction features that are known to be relevant for the light response under dark-adapted conditions.It incorporates fast buffering reactions to alter the Ca 2+ kinetics, and negative Ca 2+ feedback onto the synthesis of cyclic GMP, which is the most important feedback for dim light (Burns et al. 2002;Sakurai et al. 2011;Koutalos et al. 1995).To obtain analytic results, we performed a linear response analysis for dim light conditions where the current change is small compared to the dark current.The current response is determined by φκξ (see Eq. 5), and since collecting area κ and transduction gain ξ are much smaller in cones compared to rods (Table 1), the linear response analysis remains valid up to much higher light intensities φ in cones compared to rods.Thus, the definition of dim light depends on the photoreceptor type.For example, dim light for a cone can induce already saturating responses in a rod (Fig. 5B).We combined the analytic results with numerical simulations to obtain quantitative insight about how the various biophysical processes and the Ca 2+ kinetics determine waveform and amplitude of flash and step responses.We further investigated how the light response is affected by changing the extracellular Ca 2+ concentration.The Ca 2+ kinetics with fast buffering is governed by the effective rate μ ca = |I ex | (1+B ca )F A V os c ca (Eq.6).V os is the outer segment volume, c ca is the steady state concentration of free Ca 2+ , I ex is the steady state exchanger current, F A is the Faraday constant, and B ca is the cumulated buffering capacity that depends on the various buffer species.Because I ex is proportional to c ca for K ex c ca (see Table 1), it follows that I ex c ca only depends on intrinsic exchanger properties and their density in the membrane.In contrast, the buffering capacity of B ca depends on the Ca + affinities of the individual buffers and therefore the steady state concentration of Ca + , which fur-ther depends on the background light intensity.The most abundant Ca 2+ buffer in the outer segment is recoverin (Pugh and Lamb 2000; Kawamura and Tachibanaki 2021) with a dissociation constant K ca,rec ∼ 3μM (Chen et al. 1995) that is much larger than the dark-adapted steady state concentration of Ca + around 0.3μM (Table 1).With such conditions we conclude that B ca and the Ca + kinetics are not much affected by light as the Ca + concentration further decreases.But B ca and the Ca 2+ kinetics can be altered by adding exogenous buffers (Torre et al. 1986;Korenbrot and Miller 1989;Rieke and Baylor 1998;Field and Rieke 2002;Burns et al. 2002;Matthews 1991;Burns et al. 2002;Makino et al. 2004).
We define the waveform by normalizing a flash response with its peak amplitude.The waveform therefore characterizes the dynamics of a response (Fig. 2).With fast Ca 2+ kinetics, the Ca 2+ concentration changes in proportion to the current and the waveforms are monophasic (Fig. 3).As the Ca 2+ kinetics are slowed down due to more buffering, the Ca 2+ concentration becomes delayed with respect to the current and biphasic waveforms emerge (Fig. 3).The biphasic shape is due to a damped oscillation generated by the negative feedback interaction between Ca 2+ and cGMP synthesis.Since we only consider fast buffering, this shows that biphasic responses are not necessary an indication of slow buffering reactions.The analysis reveals that the presence of damped oscillations depends on the ratio μ ca /β d , where β d is the dark turnover rate of cGMP (Fig. 6A).This is consistent with results from a parameter sensitivity analysis in cones, which identified β d and exchanger properties (which affects μ ca ) as most important for the presence of biphasic responses (Klaus et al. 2021).With parameters from Table 1 we find that a damped oscillation is present if μ ca /β d is below an upper threshold of around 17, and above a lower threshold of around 0.08.These values depend on the strength of Ca 2+ feedback to the cyclase α 0 , and the channel cooperativity n ch (see Sect. 3.3).Oscillations disappear for μ ca → 0 (Fig. 3A) because in this limit the Ca 2+ concentration remains constant and no negative feedback to the cyclase occurs.Oscillations are also absent in GCAPs −/− mutant mice where the cyclase is not Ca 2+ -dependent (Mendez et al. 2001;Burns et al. 2002).The oscillation frequency ω and the damping rate β damp depend on μ ca (Fig. 6A,B).Whereas the oscillation frequency first increases and then decreases as a function of μ ca (Fig. 6B), the damping rate steadily increases with faster Ca 2+ kinetics (Fig. 6A).Although oscillations are present over a wide range of μ ca values, they affect the recovery phase only if the damping and therefore the Ca 2+ kinetics are not too fast.Because oscillation properties depend on the ratio μ ca /β d , the same Ca 2+ kinetics does not necessarily lead to similar waveforms because β d differs between rods, cones, and animal species (Pugh and Lamb 2000).For example, β d is around 3-fold larger in a mouse cone compared to mouse rod (Reingruber et al. 2020).The difference in β d might be a reason why oscillations have been rarely observed in amphibian compared to primate rods (Tamura et al. 1991), or why oscillations have been more frequently observed in Wt cones compared to Wt rods (Schneeweis and Schnapf 1999;Schnapf et al. 1990;Baylor et al. 1987).
With our analytic results we find that the waveform can be decomposed into a sum of exponentials that reflect the underlying biophysical processes (Eq.17 and Fig. 4).Such a decomposition cannot be performed with simulations only.We used the decomposition to study which processes limit the response recovery.For example, the recovery of flash responses in a WT mouse rod with fast Ca 2+ kinetics are limited by PDE decay (Fig. 4A), in agreement with experimental findings (Krispel et al. 2006;Tsang et al. 2006;Chen et al. 2010).In a GCAPs −/− cone the recovery is limited by β d (Fig. 4B).For these special cases it is possible to extract parameters by fitting the response recovery with a single exponential decay function.However, in general, the recovery of flash responses is not governed by a single exponential.For example, for a GCAPs −/− rod the recovery depends on the PDE decay rate μ pde and the dark turnover rate β d (see Fig. 4D in Abtout et al. (2021)).And with sufficiently slow Ca 2+ dynamics the recovery is is not an exponential decay but a damped oscillation (Fig. 4C-D).
The Ca 2+ kinetics not only affect the waveform, but also the photoreceptor sensitivity characterised by the peak current amplitude evoked by dim flashes.The amplitude of the single-photon response (SPR) increases as the Ca 2+ kinetics are slowed down because negative Ca 2+ feedback is delayed and less strong at peak time (Fig. 5A).This finding explains the larger photoreceptor sensitivity observed with exogenous buffering (Matthews 1991).The sensitivity change between Wt and GCAPs −/− photoreceptors is maximal with fast the Ca 2+ kinetics (Fig. 5A).With parameters from Table 1 we compute a maximal change around 2.3 for a rod, and around 1.8 for a cone (Fig. 5A).These values can change depending on the model for the cyclase.The experimentally observed sensitivity change between GCAPs −/− and Wt mouse rods is found to be around 2-3 (Reingruber et al. 2020;Burns et al. 2002;Mendez et al. 2001), which suggests that the physiological Ca 2+ kinetics are fast and the Ca 2+ concentration changes in proportion to the current, in agreement with (Li et al. 2020;Matthews and Fain 2003).
The peak flash amplitude î peak depends on the number of isomerisations R * 0 = κφ t produced by the flash.In the past, the empirical function î peak = 1 − e −aφ t has been often used to model the flash amplitude as a function of the flash strength φ t, and to estimate the flash sensitivity a by fitting the data (Morshedian et al. 2022(Morshedian et al. , 2017;;Astakhova et al. 2015;Korenbrot 2012;Chen et al. 2010;Tranchina et al. 1991;Schnapf et al. 1990;Nakatani and Yau 1988;Baylor et al. 1984;Schnapf et al. 1987;Baylor et al. 1987).A possible explanation for this empirical a function has been given based on a spatial model where the probability of channel closure depends on the statistical superposition of isomerisations along the outer segment (the total occlusion model) (Lamb et al. 1981), or based on a non-spatial model with pigment bleaching (Hodgkin and Obryan 1977).We showed î peak = 1 − e −aφ t is a valid approximation for the peak amplitude for fast Ca 2+ kinetics (Fig. 5B), and we computed the sensitivity as function the biophysical parameters a = κξ n ch g y, peak : κ is the collecting area, ξ is the transduction gain, n ch is the CNG channel cooperativity, and g y, peak depends on the dynamics and can be computed from Eq. 15.With parameters from Table 1 we estimate a = 0.02 for a mouse rod (a has units μm 2 photons ), which is consistent with 0.014 (Morshedian et al. 2017) and 0.026 (Chen et al. 2010) obtained by fitting experimental data.For a mouse cone we compute a much lower sensitivity a = 0.65 × 10 −5 , compatible with (Ingram et al. 2019).Since the values of g y, peak are not much different between rods and cones, the large sensitivity gap between rods and cones is due to collecting area κ and the gain ξ .Since collecting area and gain vary among species, this leads to variations in the sensitivity a.For example, a sensitivity of 0.03 has been estimated for monkey rods (Baylor et al. 1984), 2.77 × 10 −4 and 3 × 10 −5 for monkey cones (Baylor et al. 1987), ∼ 10 −3 for human cones (Schnapf et al. 1987), and 7.3 × 10 −3 for a bass cone (Korenbrot 2012).
We explored how changing the extracellular Ca 2+ concentration affects the light response in a cone (similar conclusions are obtained for a rod (not shown)).Modifying the extracellular Ca 2+ concentration alters the Ca 2+ influx and the fraction of the CNG current that is carried by Ca 2+ (Eq.20 and Fig. 7A, green line).This changes the steady state concentrations of Ca 2+ and cGMP, and the dark current (Fig. 7A).For example, a tenfold reduction in the extracellular Ca 2+ concentration reduces the Ca 2+ influx by a factor of around 7.3, and the dark-adapted Ca 2+ concentration in the OS by a factor of around 2.5.This activates the cyclase and increases the dark-adapted cGMP concentration by a factor of around 1.9, and the dark current by a factor of around 2.9.But the reduction in extracellular Ca 2+ also slows down the dynamics, increases the photoreceptor sensitivity, and reduces the amplitude of oscillations (Fig. 7B).Similar effects as described here have been observed experimentally for primate cones when the extracellular Ca 2+ concentration has been reduced around tenfold (see Fig. 2A in Cao et al. (2014)).
For a cone we also studied how the Ca 2+ kinetics affects the responses to steps of light (Fig. 8) (similar conclusions apply for a rod, not shown).The Ca 2+ kinetics affects step responses only during initial and recovery phase, whereas the intermediate plateau phase depends on steady state properties that are independent of μ ca .For dim light, the recovery phase is a mirror image of the initial phase (Fig. 8A).Hence, whereas initial and recovery phase of dim-flash responses are different and provide complementary information (Fig. 2), this is not the case for step responses.We derived an exponential approximation for steady state current îss (which is the intermediate plateau current) as function of the background light intensity φ (Eq. 25 and Fig. 8B).With Eq. 25 we estimate that the light intensity φ 1/2 for which îss = 0.5 is φ 1/2 = ln 2 κξn ch ∼ 1.3 × 10 6 photon μm 2 s , where we used parameters from Table .1.With the collecting area κ = 0.013μm 2 / photons we find κφ 1/2 ∼ 1.7×10 4 R * s visual pigment activations per second.Since the value is around tenfold smaller than 2.5 × 10 5 R * s extracted from steady state recordings (Ingram et al. 2019), this indicates that our model saturates too quickly with increasing background light intensity.The most plausible explanation is that this model lacks adaptation processes that affect the steady state at higher background light intensities, for example accelerated photopigment and PDE deactivation, and photopigment bleaching (Fain et al. 2001;Fain 2011).Whereas these processes are important for steady state computations in bright light, they can be neglected for dark-adapted flash responses.
Although in this work we have ignored slow buffering reactions, conceptually it is straightforward to include such reactions in future work, and to generalize the linear response analysis.For example, with a single slow buffering reaction, Eq. 11 will have to be replaced by 3-dimensional system of equations.Although the phase space now becomes more complex, we do not expect to find solutions that are qualitatively very different from what we have described here, consistent with simulations results in presence of slow buffer (Forti et al. 1989;Tamura et al. 1991).
We considered a non-spatial model for the outer segment to be able to derive analytic results.Such models are frequently used to study the photoresponse (Reingruber et al. 2020;Beelen et al. 2021;Astakhova et al. 2015;Invergo et al. 2014;Korenbrot 2012;Chen et al. 2010;Hamer et al. 2005;Moriondo and Rispoli 2003;Nikonov et al. 2000Nikonov et al. , 1998;;Sneyd and Tranchina 1989;Forti et al. 1989).Whereas 3D spatial models (Klaus et al. 2021;Bisegna et al. 2008), or effective 1D longitudinal models derived with the assumption of rapid radial equilibration (Reingruber et al. 2013;Lamb and Kraft 2016;Gross et al. 2012a;Pugh and Lamb 1992), provide a more accurate description of reality, the trade-off is that they are usually studied only with simulations.Nevertheless, we performed simulations with our spatial model from (Reingruber et al. 2013) to check that results are not significantly altered with a spatially extended OS (not shown).
Recently it has been found that the large outer segment of a frog rod is a spatially inhomogeneous compartment (Li et al. 2020;Mazzolini et al. 2015).It remains unclear whether this is also true for the much smaller outer segment of mouse or primate rods and cones.We leave it for future work to investigate how a possible spatial inhomogeneity affects the photoresponse.
In this work we used a deterministic mean-field model to derive analytic results that provided functional insight about how the calcium kinetics affects waveform and amplitude of the light response.We neglected noise generated by the PDE activation cascade, for example due to spontaneous PDE activations or low numbers of activated proteins.Since the equations for PDE activation are linear, our mean-field results faithfully characterize averaged responses of a model with stochastic PDE activations.The analytic results are important to understand the simulations obtained with much more complex stochastic models.In future work, stochastic modelling will show how the calcium kinetics affects the background noise and the variability of individual (not averaged) light responses (see, for example, (Hamer et al. 2003;Caruso et al. 2010;Reingruber et al. 2013Reingruber et al. , 2020))).
We validated our analytic results for different Ca 2+ kinetics using the physiological parameters from Table 1.We did not perform a more general analysis to explore rates of error between our linearized and nonlinearized model (Epstein 2018;Strogatz 2018;Ruelle 2009), for example, due to variability in the parameter values found in the literature.We leave it for future work to investigate the implications of a larger parameter space, see for example the sensitivity analysis in (Klaus et al. 2021).
Although we focused on dark-adapted photoreceptors, the analysis outlined in this work can be generalized to study flash responses in presence of a background light.
To do so, one has to first numerically compute the steady state corresponding to a background light (this can be done with a more complex model that includes adaption processes), and then use the steady state values as input for a linear response analysis.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
of open channels, and K ch is the cGMP concentration where 50% of the channels are open.The CNG current carried by Ca 2+ is I ch,ca = f ch,ca I ch , where f ch,ca is the fraction of the current carried by Ca 2+ .The influx of Ca 2+ via the CNG channels isI ch,ca2e + (inward fluxes are negative).The exchanger current as a function of the free Ca 2+ concentration is I ex = I ex,sat p ex , where p ex (c ca ) =

Fig. 2
Fig. 2 Flash waveform: simulation and analytic result.(A-D) The waveforms of the simulations in Fig. 1 are obtained by normalizing the responses to unit amplitude (black solid curves).The dashed red curves show the analytic waveform î/ î peak computed with Eq. 16 (color figure online)

Fig. 5
Fig.5Peak amplitudes of flash responses.A Change of the SPR peak amplitude as a function of μ ca /β d for a rod (black curve and diamonds) and a cone (blue curve and diamonds).The peak amplitude is normalized to the value for fast Ca 2+ kinetics.Diamonds show values extracted from simulations, solid lines are computed with Eq. 16 and R * 0 = 1.B Peak amplitude as a function of the number of isomerisations R * 0 for slow (μ ca = 2s −1 ) and fast Ca 2+ kinetics (μ ca = 50s −1 for rod and μ ca = 200s −1 for cone).Black and red colors are for rod, blue and green for cone.Diamonds show values obtained from simulations.Solid lines are computed with î peak = 1 − e −R * 0 ξ n ch g y, peak for fast Ca 2+ kinetics, andî peak = 2 f ch,ca +2 1 − e −R *0 ξ n ch g y, peak for slow kinetics (see text for explanations).The dashed lines show î peak computed with the linear approximation in Eq. 16 (color figure online)

Fig. 6
Fig. 6 Phase space analysis for Ca 2+ oscillations.A The normalized exponential decay rates β 1 = β d λ 1 (black line) and β 2 = β d λ 2 (red line) for the current component îβ in Eq. 17 are shown as a function of μ ca β d

Fig. 8
Fig. 8 Responses to steps of light for a cone.A The dim-light current waveform from Eq. 24 for steps with duration t = 1s and for various values of μ ca as indicated in the legend.B The steady state current from Eq. 25 (blue curve) as a function of the light intensity φ is compared to results obtained by numerically computing the steady state with Eq. 8 (black line with markers).C and D Simulations of step responses (solid lines) are compared to the corresponding analytic result computed with Eq. 26 (dashed lines) for fast C and slow D Ca 2+ kinetics.The light stimulations in C and D are identical (color figure online) . Steady state values can either be computed with Eq. 3, or experimental values can be used.We replace the parameters α max , I ch,max and I ex,sat (which are not well known) with the steady state quantities c ca,0 , c cg,0 and I 0 .For example,

Table 1
Parameter values for the mouse rod and cone model Damped oscillations occur for 0.08 < μ ca β d < 16.9 (seeSect.3.3).In this range the eigenvalues λ 1 and λ 2 are complex conjugate, and the damping rate is β damp ≈ β d The normalized frequency ω of the damped oscillation as a function of μ ca