Introduction

A central goal in neuroscience is to understand how behaviorally-significant computations are implemented by cellular and synaptic mechanisms (see (Churchland and Abbott, 2016)). Such computations typically rely on a diverse collection of interacting circuit mechanisms, many of which are nonlinear (Gollisch and Meister, 2010; Isaacson and Scanziani, 2011; Jadzinsky and Baccus, 2013; Zador, 2000). Circuit outputs - which are often most amenable to experimental measurement - typically do not uniquely identify the role of specific circuit mechanisms. Hence progress towards understanding the mechanistic basis of circuit function requires better tools to manipulate specific mechanisms and identify their contribution to circuit output.

Sensory systems provide a clear example of these issues. Sensory processing is strongly shaped by the properties of the sensory receptors themselves and by post-receptor circuit mechanisms. Yet separating the contributions of receptor and post-receptor mechanisms to circuit outputs can be difficult, and computational models often make untested assumptions about their relative importance. For example, models for retinal processing generally assume that photoreceptor responses are linear or near-linear, despite decades of direct evidence that this is not the case (reviewed by (Burns and Baylor, 2001; Schwartz, 2021)). Adaptation to changes in mean light intensity provides a specific and important example. Both photoreceptors and post-photoreceptor circuits adjust response properties such as gain and kinetics to match the prevailing light inputs (reviewed by (Demb, 2008; Dunn and Rieke, 2006)), yet models for retinal outputs typically start by passing light inputs through a linear spatio-temporal filter; this architecture includes linear-nonlinear (Chichilnisky, 2001), generalized-linear (Pillow et al., 2008) and computational neural network (CNN) (McIntosh et al., 2016; Turner et al., 2019) models. Models with an initial linear filtering stage cannot capture the spatially local adaptation produced by photoreceptors. Tools that separate receptor and post-receptor contributions to adaptation are needed to better understand how this salient aspect of retinal processing works.

Here, we show how we can predictably manipulate the responses of rod and cone photoreceptors to causally probe their role in shaping responses of downstream visual neurons. We build on decades of work identifying and testing biophysical models for rod and cone phototransduction (Nikonov et al., 2000; Pugh and Lamb, 1993; Rieke and Baylor, 1998a; Younger et al., 1996). We show that existing phototransduction models, with appropriate parameters, can account for responses of rod and cone photoreceptors from both primate and mouse to a broad range of inputs. The resulting forward models can be used as a front-end for encoding models, resulting in a better accounting of early time-dependent nonlinearities. We then show that these models can be mathematically inverted - enabling the design of stimuli that will give rise to photoreceptor responses with specific desired properties such as a linear dependence on light input. Direct recordings of responses to these stimuli show that they work as designed. This approach provides a tool that can be used to predictably shape the photoreceptor responses and hence to casually dissect the impact of specific properties of the phototransduction currents on downstream visual signaling and perception.

Results

Biochemical model

We chose a model architecture based on the biochemical interactions that comprise the phototransduction cascade. We chose this model rather than empirical models (e.g. (Clark et al., 2013)) because it was more clearly connected to the mechanistic operation of the phototransduction cascade and because it was exactly invertible (see (Angueyra et al., 2022) for comparison of the biochemical model with empirical models).

Combined work in biochemistry, molecular biology and physiology has produced a clear understanding of the operation of the phototransduction cascade (Figure 1A; reviewed by (Burns and Baylor, 2001; Rieke and Baylor, 1998b; Schwartz, 2021). In brief, the phototransduction current flows through cGMP-gated channels in the photoreceptor outer segment. Levels of cGMP, and hence the phototransduction current, are maximal in darkness. Light activates G-protein coupled receptors, which activate G-proteins and a cGMP-phosphodiesterase. The resulting decrease in cGMP concentration allows some cGMP-gated channels to close and decreases the current. Levels of cGMP are restored by synthesis of cGMP by a guanylate cyclase; cyclase activity is sped by a calcium-feedback pathway that strongly shapes light responses (Burns et al., 2002).

Model and fitting procedure. A. Phototransduction cascade and differential equations that describe operation of key components. Abbreviations are as follows: Stim: stimulus; R: receptor (photopigment) activity; P: phosphodiesterase activity; G: cGMP concentration; S: rate of cGMP synthesis; I: membrane current; C: calcium concentration. B. Fits (red) of a model with four free parameters (γ, σ, η, KGC) to the responses of a mouse rod (see Table 1). All measured responses are from the same rod, and all model responses use the same values for the free parameters. Responses to different flash strengths or flash timing have been displaced vertically in the left two panels. The weakest flash in the far left panel produced on average ∼3 activated rhodopsin molecules, and each successive flash was twice as bright. The bottom trace in the middle panel was to a light step alone, while the other traces in this panel each had 5 superimposed flashes. The horizontal scale bar for the inset in B is 1 sec.

Parameters and best fit consensus values for biophysical phototransduction models for each cell type (see Methods for fitting details).

The biochemical understanding summarized in Figure 1A has led to several quantitative models for phototransduction (Nikonov et al., 2000; Pugh and Lamb, 1993; Rieke and Baylor, 1998a; Younger et al., 1996). We focus on a set of differential equations that form the core of these models. Our goal was to develop a model that generalized across stimuli rather than one in which the parameters accurately reflected the operation of components of the cascade. Below we show that the model illustrated in Figure 1A, with appropriate parameters, captures the responses of rod and cone photoreceptors to a variety of stimuli. Most importantly for our purposes here, the models are sufficiently accurate to enable us to accurately predict stimuli that elicit desired photoreceptor responses, as tested in Figures 6-9.

Our phototransduction model has a total of eleven parameters (Figure 1A). Two of these (q and Smax) can be expressed in terms of others using steady state conditions - e.g. that dC(t)/dt = 0 in steady state. Several others were fixed based on previous measurements and their relatively minor impact on model predictions; this included k, n, and m for both rods and cones (Koutalos et al., 1995b; Rieke and Baylor, 1996) and β for rods (Field and Rieke, 2002). The dark cGMP concentration was estimated from the measured dark current for each recorded cell. The receptor and PDE decay rates (σ and Φ) had near identical effects on predictions (see (Pugh and Lamb, 1993)) and were constrained to be equal. These considerations allowed us to reduce the number of free model parameters to 4 for rods and 5 for cones. This choice reflects a balance of minimizing the number of free model parameters while providing sufficient flexibility to capture key aspects of the light response.

The free parameters consisted of the gain γ with which photons are converted to photopigment activity, the photopigment decay rate σ, the spontaneous PDE activation rate η, the constant KGC specifying the calcium sensitivity of the cGMP synthesis rate, and (for cones) the rate β of Ca2+ extrusion. This choice of free parameters is not unique given the similar impact of some parameters on model output (such as σ and Φ). But, as indicated below, it provided sufficient flexibility to capture measured responses to a variety of stimuli. We numerically identified optimal values for these parameters separately for each photoreceptor type (see Methods).

We were able to collect responses to a larger set of stimuli from individual rods than cones due to differences in recording techniques. We recorded from rod photoreceptors using suction electrodes, which allow for stable, long-lasting recordings. This allowed us to measure responses to a comprehensive set of stimuli from each individual rod. We recorded from cone photoreceptors using whole-cell patch clamp recordings. Light responses are stable in these recordings for a much shorter time period (due to internal dialysis), and consequently we recorded responses to a single stimulus type from each individual cone. Despite these differences in data collection, we identified values for the free parameters in both rod and cone models by numerically minimizing the mean-squared error between model predictions and measured responses (see Methods for more details).

Fitting rod phototransduction models

We measured responses of mouse and macaque rod photoreceptors to a set of stimuli consisting of a family of brief flashes of different strengths, flashes delivered at several times relative to the onset and offset of a light step, and a “variable mean noise stimulus” consisting of Gaussian noise with periodic changes in mean intensity (Figure 1B; see Methods for experimental details). The variable mean noise stimulus approximates the large and frequent changes in input that are characteristic of natural vision (Frazor and Geisler, 2006). These large (up to 30-fold) changes in mean intensity strongly engage photoreceptor adaptation, and the superimposed Gaussian noise probes sensitivity to rapidly varying inputs.

We recorded responses to this collection of stimuli for each rod in our data set. We fit mouse and macaque rods separately, in each case seeking a common model (i.e. with the same parameters) that fit responses across stimuli. We identified model parameters by numerically minimizing the mean-squared error between the model predictions and measured responses across stimuli (see Methods).

Figure 1B compares the predicted and measured responses of a mouse rod photoreceptor. Predicted and measured responses are not identical - e.g. the model predicts a faster initial rising phase of the flash response than that measured and underestimates the amplitude of flashes delivered shortly after a light step. Nonetheless, a common model specified by four free parameters (γ, σ, η, and KGC) accounts for many aspects of the responses to this set of stimuli. These parameters were well-constrained by fitting the measured responses: best fit parameters varied by ∼10% across different photoreceptors, and the mean-squared error increased significantly for 10-20% changes in parameters (see Methods for a more detailed analysis).

To arrive at consensus models for rod transduction currents, we fit measured responses from multiple rods from either primate or mouse simultaneously. We first used the measured dark current and the previously-measured relation between current and cGMP (Eq. 4 in Figure 1; (Rieke and Baylor, 1996)) to specify the cGMP concentration in darkness for each cell. We then allowed γ to vary between cells to account for differences in sensitivity. The remaining parameters (σ, η, and KGC) were constrained to have the same value across cells. We then numerically identified best values for these parameters by minimizing the mean square-error between model predictions and the set of recorded rod responses. These consensus models provide our best estimate of how an unrecorded cell of a given type will respond; these models provide the foundation for the manipulations of the photoreceptor responses in Figures 5-10. Consensus parameters for mouse and primate rods were very similar (see Table 1).

Fitting cone phototransduction models

Due to the limited duration of the cone recordings, we focused on the variable mean noise stimulus. We focused on this stimulus because it was the most effective in constraining cone model parameters, likely because it strongly engages adaptation and thoroughly probes response dynamics.

We identified consensus cone models by fitting measured responses to the variable mean noise stimulus from multiple cones simultaneously, as described above for the rods. Four parameters (σ, η, KGC and β) were constrained to have the same value across cones, while γ was allowed to vary across cells. Consensus parameters for mouse and primate cones differed substantially, reflecting the ∼2-fold faster kinetics of primate cone responses (flash response time-to-peak is ∼35 ms in primate cones, and ∼70 ms in mouse cones; see Table 1). For primate cones, performance using these consensus parameters was very similar to that from the parameters in (Angueyra et al., 2022); for consistency we used these published parameters for our consensus model.

Model performance

Figure 2 tests the performance of the consensus models for both rod and cone photoreceptors. We focused on the variable mean noise stimulus since it probes responses to rapid stimulus variations and large changes in mean intensity that strongly engage photoreceptor nonlinearities. For each individual cell, we specified the dark cGMP concentration (GD) using the measured dark current and allowed the overall gain γ to vary; the remaining parameters were set to the consensus model values in Table 1. The resulting models (red) captured the measured responses (black) well.

Photoreceptor model and fits. A. Comparison of measured responses (black) with predictions of full (red) and linear model (blue) to variable mean noise stimulus (gray). Full model responses use consensus parameters from fitting responses of multiple cells of each type simultaneously, with the dark current and sensitivity allowed to vary between cells (see Methods). The linear model was generated from fitting the low-contrast responses of the full model (see Methods). Insets expand regions in gray boxes. Linear model parameters (see Eq. 7) were α = 0.31, τR = 10.6 ms and τD = 23.6 ms for the primate cone, α = 0.031, τR = 31.6 ms and τD = 56.7 ms for the mouse cone, α = 5.3, τR = 141 ms and τD = 208 ms for the primate rod, and α = 4.7, τR = 115 ms and τD = 185 ms for the mouse rod, B. Fraction of variance explained for the full model fit to each cell individually (y-axis) plotted against that for the consensus model that has fixed parameters across cells except for the dark current and sensitivity. Means ± SDs were 0.90 ± 0.04 (specific model) and 0.87 ± 0.05 (consensus model) for 6 primate cones, 0.92 ± 0.02 and 0.94 ± 0.02 for 6 mouse cones, 0.94 ± 0.02 and 0.94 ± 0.02 and 0.94 ± 0.02 for 6 primate rods and 0.93 ± 0.03 and 0.93 ± 0.03 for 8 mouse rods.

For each recorded cell, the consensus models captured more than 80% of the variance of the responses (Figure 2B); on average these models captured ∼90% of the variance. We compared the performance of the consensus models with that of models fit to each cell individually. Consensus models and models fit to individual cells performed similarly (Figure 2B). Models also generalized well across stimuli and to cells that did not contribute to the fits (Figure 2 - Figure Supplement 1). This suggests that the consensus models can be used to predict responses of unrecorded cells.

Figure 2 also shows predictions of a linear phototransduction model (blue) (see Methods). Rod and cone phototransduction currents depend linearly on light intensity for low contrast inputs (e.g. (Baylor et al., 1984; Hass et al., 2015; Schnapf et al., 1990)). Hence, the linear model was fit to responses of the full model to low contrast inputs at the mean intensity of the variable mean noise stimulus (25 R*/rod/s for rods and 22,000 R*/cone/sec for cones; see Methods for model construction details); this means that the differences between the linear and full model responses can be directly attributed to nonlinearities in phototransduction. Linear predictions deviated from the measured responses considerably more than did predictions of our biochemical model (Figure 2 - Figure Supplement 2). Specifically, linear models overpredicted responses at high mean light levels, and underpredicted responses at low mean light levels.

Some of these issues could be resolved by adding a free scale factor (e.g. in Figure 2 - Figure Supplement 2) or a time-independent nonlinearity to the linear model (i.e. a linear-nonlinear model); however, responses to many stimuli show clear time-dependent nonlinearities that are not captured by linear-nonlinear models (e.g. Figure 5; see also (Angueyra et al., 2022)). The biochemical model, though not perfect, tracked the measured responses quite well.

Figures 1 and 2 indicate that the known components of the phototransduction cascade, with appropriate parameters, can account for the full time course of the rod and cone photocurrents, including the nonlinearities in the photoreceptor responses. These models will not capture slow forms of adaptation that might become important, for example, at high light levels; for the range of stimuli used here, such slow adaptation contributes little to the photoreceptor currents (e.g. (Angueyra et al., 2022)). Our central goal was to use these models to design stimuli that permit predictable manipulations of the photoreceptor responses. Hence, more important than the accuracy of the forward model is whether the stimulus predictions that model generates are accurate. We directly evaluate these stimulus predictions below.

Model Inversion

In this section, we show that the model illustrated in Figure 1A can be exactly inverted to identify the stimulus that corresponds to a specific desired photoreceptor response. The resulting unique correspondence between photocurrent and stimulus applies to the entire stimulus waveform, including the mean intensity and modulations about the mean. In other words, the result of model inversion is an estimate of the full stimulus in real units (isomerizations per second).

Model inversion requires a one-to-one relationship between inputs and outputs - i.e. for every input there must be a unique output and vice versa. Mathematically, this holds for the phototransduction model summarized in Figure 1A: the model consists of linear elements (Eqs. 1, 2 and 5), static nonlinear elements (Eqs. 4 and 6) and a time-dependent nonlinear element (Eq. 3). The linear elements can be described as linear filters, and inverted by deconvolution. The static nonlinear elements consist of one-to-one mappings of inputs to outputs, and can be inverted correspondingly. The time-dependent nonlinear component can be rearranged and solved analytically (see Methods). The result is an analytical relationship between the phototransduction current and the stimulus - in other words an inverse model that takes the photocurrent as input and identifies the corresponding stimulus.

Figure 3 illustrates the steps in this procedure, using the current response generated by the forward model to the variable mean noise stimulus as an example “target” (far right). The first step is to convert this target current response to corresponding time-varying concentrations of cGMP and calcium using Eqs. 4 and 5 (step 1). In step 2, the time courses of cGMP and calcium are used to determine the time-varying PDE activity using Eqs. 3 and 6. In the final step (step 3), the PDE activity is used to determine the stimulus through Eqs. 1 and 2. In this case, as expected, we recover exactly the stimulus that we started with, confirming that mathematically the phototransduction model can be inverted.

Steps in model inversion. A test current was generated from a variable mean noise stimulus using the full model (far right). Step 1 (right) converts this current into changes in cGMP and calcium using Eqs. 4 and 5 (see Figure 1). Step 2 converts the time course of the cGMP and calcium into that of the PDE activity using Eqs. 3 and 6. Finally, step 3 converts the PDE to the stimulus using Eqs. 1 and 2. The estimated stimulus is identical to the initial stimulus because there is no added noise and the inversion process is exact.

The procedure illustrated in Figure 3 is exact in the context of the phototransduction model - i.e. given the architecture of the model, there is a one-to-one mapping of inputs to outputs. In practice, we would like to use this procedure to predict stimuli that will elicit specific photoreceptor responses. These model-based predictions could fail if the forward model does not accurately capture aspects of the real photoreceptor responses that are important for inversion. The model, for example, may capture some temporal frequencies of the response better than others, and this could limit the accuracy of the model inversion at those frequencies that are not captured well.

To test how well the model inversion process works on real responses, we applied it to measured photoreceptor responses to the variable mean noise stimulus. Noise in the measured responses could cause the estimate to be better at some temporal frequencies than others; this is particularly true at high temporal frequencies that are poorly encoded by the photoreceptors. The limited ability to recover spatial information in deconvolution microscopy provides an example of how noise can limit model inversion even with a good forward model. Noise in deconvolution microscopy is controlled by choosing not to try to recover spatial frequencies at which the noise is too high. We took a similar approach: we limited the frequency content of the stimulus to temporal frequencies that we could reasonably expect to recover (0-60 Hz for primate cones, 0-30 Hz for mouse cones, 0-15 Hz for both primate and mouse rods) and required that the power spectrum of our estimate matched that of the true stimulus (see Methods). These constraints reduced high frequency noise in the estimated stimuli. Figure 4A shows stimulus estimates based on measured responses for each photoreceptor type. The estimates closely approximate the actual stimulus, including recovering rapid stimulus fluctuations (insets). For each photoreceptor type, stimulus estimates based on model inversion captured the majority of the stimulus variance (Figure 4B).

Test of model inversion based on measured responses. A. Stimulus (gray in top panel), measured response (black in lower panels) and estimated stimulus (blue in top panels) calculated by using the measured response as input to the inverse model as in Figure 3. Estimates are able to recover both the periodic changes in mean intensity and the more rapid superimposed stimulus modulations (insets). B. Variance explained for stimulus estimates based on the average response across multiple stimulus trials compared to that based on individual responses. Since the model captures only the deterministic part of the response, noise in the individual responses lowers the accuracy of the estimates and causes the points to fall below the unity line. This effect is modest but systematic. Means ± SDs were 0.91 ± 0.02 (single epochs) and 0.92 ± 0.02 (mean) for 6 primate cones, 0.66 ± 0.08 and 0.72 ± 0.06 for 6 mouse cones, 0.74 ± 0.03 and 0.81 ± 0.04 for 6 primate rods, and 0.68 ± 0.04 and 0.77 ± 0.05 for 8 mouse rods.

Variance explained, as in Figure 4B, mixes low temporal frequencies for which we might expect inversion to work well with high temporal frequencies for which it is expected to do less well. To explore the accuracy of the inversion across temporal frequencies, we compared the power spectrum of the stimulus with that of the residuals given by the difference between the stimulus and estimate (Figure 4 - Figure Supplement 1). At low temporal frequencies, stimulus power was much greater than that of the residuals, consistent with accurate stimulus estimation at these frequencies. This difference decreased with increasing frequency, and the two power spectra converged near 30 Hz for primate cones, 10 Hz for mouse cones, and 2-3 Hz for both primate and mouse rods. This provides an estimate of the frequency range over which stimulus inversion is possible for each photoreceptor type.

The ability to convert photoreceptor responses to input stimuli is useful in several ways (see Discussion). Our focus in the remainder of the Results is using it to design stimuli that alter the photoreceptor responses in specific ways, such as negating the impact of adaptation and shaping response kinetics.

Light-adaptation clamp examples

Figure 6 illustrates how we can use the model inversion procedure to predictably manipulate photoreceptor responses, in this case using the cone model and (for simplicity) a sinusoidal stimulus. Unlike Figure 4 where we invert measured responses, we now start with a desired or target response. For the specific applications described below, the aim was to remove nonlinearities in photoreceptor responses and hence we generate the target response from a linear phototransduction model. As in Figure 2 above, the linear model used to generate the target response was obtained by fitting responses of the full phototransduction model to low contrast noise stimuli at a specified mean light level (see Methods). The output of this linear model provides an estimate of how the photoreceptors would respond if we could eliminate their nonlinear properties. The response of the linear model to a sinusoidal stimulus is, as expected, another sinusoid (Figure 6, second panel from left). The response of real cones, and of the full cone model, deviates strongly from a sinusoid (black trace in the far right panel; (Angueyra et al., 2022)).

Adaptive nonlinearities in phototransduction cause the responses to sinusoidal light inputs to deviate from sinusoids. These time-dependent nonlinearities operate sufficiently rapidly to shape each cycle of the response, causing responses to decreases in light intensity to be larger than those to increases and causing the responses to light-to-dark transitions to be slower than those to dark-to-light transitions. To eliminate the impact of these nonlinear response properties, we seek a stimulus that will cause the response of the full model to match the target response - i.e. a stimulus that “clamps” the cone response to the target. The model inversion process identifies such a stimulus directly (red trace in second panel from right). Our primary interest is not in the shape of the stimulus itself but instead in the response that the stimulus produces. The right panel in Figure 6 compares measured cone responses to the original (black) and modified (red) stimuli; responses to the modified stimulus are considerably more sinusoidal than those to the original stimulus, as quantified in more detail below.

The approach outlined in Figure 6 is exact in principle, but could fail due either to inadequacies of the model or to cell-to-cell variability in the photoreceptor responses. Hence, a key step is testing whether the calculated stimuli produce the expected responses. Figures 7-9 test this approach quantitatively for several stimulus manipulations.

Errors in stimulus estimation. Power spectra of stimulus and residual (stimulus - estimate) for example primate and mouse rods and cones (top) and averages across cells (bottom). The stimulus and residual power spectra diverge strongly at low frequencies where the estimates are close to the actual stimulus, and converge at higher frequencies when the estimate becomes poor. The convergence point depends on photoreceptor type, as expected for the different kinetics of rod and cone responses.

Light-adaptation clamp procedure. Starting with an initial stimulus (left), we generate a target or desired response (second panel from right). In this case, the target was chosen to be the response of a linear phototransduction model; for a sinusoidal input, a linear model produces a sinusoidal output. We use the (linear) target response as the input to the inverse phototransduction model and identify the stimulus required to elicit that response (red in third panel from left). Substantial stimulus modifications are required for the full model to produce a sinusoidal output. Finally, we confirm that the modified stimulus works as designed in direct recordings, in this case from a primate cone (right panel). Linear model parameters used to compute the target response (see Eq. 7) were α = 0.31, τR = 10.6 ms and τD = 23.6 ms.

Compensating for increment/decrement asymmetries in responses to sinusoidal stimuli. A. Photoreceptor responses (bottom, black) to original sinusoidal stimuli (top, black) and modified stimuli (red). Dashed lines in the bottom panels are best-fitting sinusoids for reference. Original and modified stimuli are shown above the recorded responses. Linear model parameters used to compute the target responses (see Eq. 7) were α = 0.31 (20000 R*/s) and 1.65 (5000 R*/s), τR = 10.6 and 15.2 ms and τD = 23.6 and 18.1 ms for the primate cone, α = 0.031 and 0.13, τR = 31.6 and 26.7 ms and τD = 56.7 and 50.3 ms for the mouse cone, α = 5.3, τR = 141 ms and τD = 208 ms for the primate rod, and α = 4.7, τR = 115 ms and τD = 185 ms for the mouse rod. B. Mean-squared-error between measured responses and best fitting sinusoids for modified stimuli plotted against that for original stimuli for each recorded cell. Points with error bars are means ± SDs. These values are 0.15 ± 0.05 (modified) and 0.02 ± 0.02 (original) at 20,000 R*/s and 0.08 ± 0.06 and 0.03 ± 0.02 at 5,000 R*/s for 6 primate cones, 0.10 ± 0.01 and 0.04 ± 0.03 (20,000 R*/s) and 0.06 ± 0.03 and 0.02 ± 0.01 (5,000 R*/s) for 5 mouse cones, 0.14 ± 0.08 and 0.05 ± 0.01 for 6 primate rods and 0.16 ± 0.04 and 0.03 ± 0.01 for 7 mouse rods.

Compensating for adaptation produced by changes in mean light level. A. Responses to a brief light flash delivered before and during a step in light intensity. For the original stimuli, flashes delivered before and during the step are identical, and the resulting responses decrease in amplitude ∼2-fold (summarized on x-axis of bottom panels). Red traces show responses to stimuli designed to compensate for the adaptation produced by the change in light intensity (following approach in Figure 6). Linear model parameters used to compute the target responses (see Eq. 7) were α = 1.65, τR = 15.2 ms and τD = 18.1 ms for the primate cone, α = 0.13, τR = 26.7 ms and τD = 50.3 ms for the mouse cone, α = 5.3, τR = 141 ms and τD = 208 ms for the primate rod, and α = 4.7, τR = 115 ms and τD = 185 ms for the mouse rod. B. Summary of gain changes (amplitude of response during the step divided by that of response prior to step) for responses to modified (y-axis) and original (x-axis) stimuli. Open circles show the mean ± SD; values are 1.2 ± 0.2 (modified) and 0.52 ± 0.06 (original) for 10 primate cones, 0.8 ± 0.1 and 0.32 ± 0.05 for 5 mouse cones, 1.3 ± 0.3 and 0.6 ± 0.1 for 7 primate rods, and 1.0 ± 0.1 and 0.47 ± 0.08 for 7 mouse rods.

Manipulating kinetics of photoreceptor responses. A. Responses of primate cones and mouse rods to brief flashes (black) and stimuli designed to slow down responses slightly (left) and more substantially (right) (red). Traces show the mean ± SEM (4 cells for upper left, 3 for upper right, 5 for both bottom panels). Thin traces show the target responses used to generate the modified stimuli. B. As in A for stimuli designed to speed up responses (measured traces show mean ± SEM for 4 cones in the top panel and 3 rods in the bottom panel). C. Measured change in time-to-peak plotted against predicted change for a variety of manipulations of kinetics as in A and B. Linear model parameters used to compute the target responses (see Eq. 7) were α = 1.65, τR = 15.2 ms and τD = 18.1 ms for the primate cone, α = 0.13, τR = 26.7 ms and τD = 50.3 ms for the mouse cone, α = 5.3, τR = 141 ms and τD = 208 ms for the primate rod, and α = 4.7, τR = 115 ms and τD = 185 ms for the mouse rod.

Sinusoidal stimuli

We start with the sinusoidal stimuli illustrated in Figure 6. As predicted by the model, responses of both rod and cone photoreceptors to high-contrast sinusoidal light inputs are strongly nonsinusoidal: responses to decreases in light intensity are larger than those to increases, and responses to dark-to-light transitions are faster than those to corresponding light-to-dark transitions (Figure 7A). These nonlinearities in the photoreceptor responses are clear from the substantial deviations between the measured responses and the best-fitting sinusoids (dashed lines in Figure 7A). These response asymmetries are important in interpreting responses of downstream visual neurons to similar stimuli. For example, asymmetries in signaling of On and Off ganglion cells are often attributed to differences in On and Off retinal circuits (e.g. (Stockman et al., 2014)), and the possibility that they arise at least in part in the photoreceptors themselves is rarely considered.

Responses to the modified stimuli (red traces in Figure 7A) were much closer to sinusoidal than responses to the original sinusoidal stimuli (black traces). Specifically, asymmetries between increment and decrement responses and light-to-dark vs dark-to-light transitions were substantially reduced for the modified stimuli. We quantified these deviations of the measured responses from a sinusoid for both original and modified stimuli by computing the mean-squared error (MSE) between the recorded responses and the best-fit sinusoid. Sinusoidal fits to responses to the modified stimuli had considerably lower MSE than those to the original stimuli in each rod and cone tested (Figure 7B). Neither the sinusoidal stimuli nor the specific cells that contribute to Figure 7 were used in fitting the phototransduction models. Hence Figure 7 tests the ability of the model-based inversion process to generalize across photoreceptors and across stimuli.

Figure 7 demonstrates that the adaptive nonlinearities in phototransduction that distort responses to sinusoidal stimuli can largely be compensated for, resulting in near-sinusoidal responses. These or similar modified stimuli could in turn be used to test the contribution of nonlinearities in phototransduction to asymmetries in responses of On and Off retinal neurons to sinusoidal stimuli such as drifting or contrast-reversing gratings.

Steps and flashes

As a second example of the light-adaptation clamp, we generated stimuli that compensate for the adaptive changes in gain of photoreceptor responses produced by changes in mean light intensity. These gain changes cause responses to a fixed strength flash to become smaller and faster when the mean light intensity increases.

The control stimulus consisted of identical flashes delivered at two mean light levels. As expected, responses to this stimulus (black traces in Figure 8A) showed clear effects of adaptation. Specifically, the response to the flash decreased ∼2-fold at the higher mean light level (summarized in Figure 8B). We then passed this original stimulus through a linear model to estimate the photoreceptor response without adaptation. This linear prediction provided the target response in the stimulus identification process illustrated in Figure 6; the resulting modified stimuli are shown in red at the top of Figure 8A. None of the cells used to test these modified stimuli contributed to the model fitting, so this again tests the ability of the models to generalize across photoreceptors.

The red traces in Figure 8A (bottom) show responses to the stimuli designed to compensate for adaptation. Responses to the test flashes had similar amplitudes at the low and high mean light level. This similarity is held across cells (the ratio of the response amplitudes is plotted on the y-axes in Figure 8B). The change in gain produced by the background was reduced in all photoreceptor types. Responses to the modified stimuli in some cases (e.g. the mouse cones) show a systematic dependence on background, likely due to differences in sensitivity between the cells used to fit and test the model. Nonetheless, the modified stimuli compensated for much of the reduction in response gain produced by adaptation.

Altering kinetics

As for adaptation, separating the contributions of phototransduction and post-transduction circuit mechanisms to the kinetics of responses of downstream cells is difficult. Hence we sought to use the model inversion process to predictably alter the kinetics of the photoreceptor responses, which could help elucidate the impact of photoreceptor kinetics on downstream signaling. For example, compensating for the speeding of photoreceptor responses that occurs with increasing light intensity would isolate the kinetics of the post-transduction circuitry and determine if they change independently from the changes in kinetics of the photoreceptor responses.

Slowing down the kinetics is conceptually straightforward - we slow down the kinetics of the stimulus itself and that will also cause the response to slow down. We apply this to responses to brief flashes. To test if we can do this in a predictable manner, we stretched the time axis of the linear responses to a brief flash by factors of 1.25-4 by scaling the time constants in Eq. 7. We used the resulting response as the target for the model inversion as in Figure 6. We omitted mouse cones because these experiments required longer-lasting recordings with stable kinetics, which were difficult to achieve with these cells.

The top panels in Figure 9A show original (black) and modified (red) stimuli designed to slow down the photoreceptor response slightly (left) and more substantially (right). The bottom panels compare measured (thick traces) and predicted (thin traces) responses to the original flash and the modified stimuli. Changes in time-to-peak of the measured responses to the modified stimuli follow the predictions well for both rod and cone photoreceptors (summarized in Figure 9C).

Speeding up the responses is similarly straightforward in principle but more difficult in practice. Stimuli that speed up the response consist of a brief increment followed by a decrement - which together cause all but the initial rising phase of the response to cancel. In practice, this approach is limited because the contrast of the decrement cannot exceed 100%. Nonetheless, speeding up responses by 20-30% is possible, and stimuli predicted to achieve this are shown in Figure 9B. Measured responses sped as predicted, although the falling phase of the response, particularly for mouse rods, was faster than the predicted response. This again reflects differences in the specific cells used to fit and test the model. As for slowing responses, the time-to-peak of the measured responses closely follows predictions (Figure 9C).

Impact of cone adaptation on ganglion cell responses

The overall goal of this project is to separate the contributions of phototransduction and post-transduction circuit mechanisms to downstream visual signaling. Here, we show two examples of how the ability to manipulate photoreceptor responses can reveal adaptational mechanisms at different locations in retinal circuits.

First, we compare responses of primate cones and downstream retinal neurons to the step plus flash stimulus used in Figure 8. Figure 10A shows responses of a primate cone, a horizontal cell, and an On parasol ganglion cell to the original step plus flash stimulus (black) and the modified stimulus that negates the cone adaptation (red). The stimuli used for the three cells illustrated were identical. As shown in Figure 8, the modified stimulus effectively negates adaptation in the cone transduction currents, and the responses to the cone responses to the two flashes are similar. Negating adaptation in the cones also eliminated most or all of the adaptation in downstream responses (summarized in Figure 10B). Ratios exceeding one are likely due to differences in sensitivity between the cells used to fit and test the model. Previous work reached a similar conclusion by quantitatively comparing adaptation measured in cones and retinal ganglion cells (Dunn et al., 2007). Negating cone adaptation provides an alternative approach to this issue and similar questions about the contributions of photoreceptor and post-photoreceptor adaptation.

Cone and post-cone adaptation. A. Left: schematic of retinal circuit. Right: Responses of a cone (top), horizontal cell (middle) and On parasol ganglion cell (bottom) to the step and flash protocol for both original and modified stimuli. Stimuli were delivered from a computer monitor at a refresh rate of 60 Hz and hence appear different from those in Figure 8. B. Summary of experiments like those in A, plotting the change in gain for the modified stimuli (the ratio of the amplitude of the response to the flash during the step that before) against that for the original stimuli. Means ± SDs were 1.2 ± 0.2 (modified) 0.52 ± 0.06 (original) for 10 cones and 1.2 ± 0.1 and 0.6 ± 0.2 for 5 parasol RGCs. C. Left: horizontal cell responses to a sinusoidal stimulus (black) and stimulus modified to generate sinusoidal cone responses (red, see also Flgures 6 and 7). Right: trajectories of responses to one stimulus cycle. The response to the original stimulus follows different trajectories during the depolarizing phase compared to the hyperpolarizing (c→d→a) phase. D. Summary of results from three horizontal cells and three cone bipolar cells. Time-dependent nonlinearities were measured from the area enclosed by response trajectories as in C. Means ± SDs were 0.08 ± 0.05 (area modified) and 0.22 ± 0.08 (area original) across all 6 cells. Linear model parameters used to compute the target responses (see Eq. 7) were α = 1.65, τR = 15.2 ms and τD = 18.1 ms

Second, we recorded responses of primate horizontal cells and cone bipolar cells to sinusoidal stimuli and stimuli designed to produce sinusoidal cone transduction currents (as in Figures 6 and 7). Horizontal responses to both stimuli showed clear nonlinearities (Figure 10C, left); specifically, responses to the increment and decrement phases of the stimuli were highly asymmetric. The depolarizing and hyperpolarizing phases of the responses, however, appeared more symmetrical for the modified stimuli (red). The change in shape can be visualized by plotting the trajectories of one cycle of the responses vs the stimulus (Figure 10C, right). For the original sinusoid, the depolarizing and hyperpolarizing (c→d→a) response phases differed considerably. These differences were much smaller for the modified stimuli. This difference in response shape held across cells (Figure 10D). In this case, compensating for nonlinearities in cone phototransduction reveals an additional, largely time-independent, nonlinear process shaping horizontal and cone bipolar responses. Further, it disentangles contributions of phototransduction and post-transduction processes to time-dependent and time-independent nonlinearities in the horizontal and bipolar responses.

Figure 10 shows two relatively simple examples of how the stimulus design approach developed here can be used. We return to other possible applications in the Discussion.

Discussion

Rod and cone photoreceptors play an essential role in the responses of downstream visual neurons and visual perception. Refining this picture to identify the specific contributions of photoreceptors and post-photoreceptor circuitry to the computations that underlie vision has been difficult. Here, we have created and tested a tool that allows predictable manipulations of the photoreceptor responses to reveal their role in downstream signaling. Direct tests indicate that this tool is effective, and that it generalizes well across photoreceptors and stimuli. We envision this being useful to causally test how responses of downstream cells or perception are shaped by specific aspects of the photoreceptor responses such as adaptation. We have described in Methods procedures that should help translate this tool to other laboratories.

Limitations

Model validation

The procedure that we use to invert the phototransduction cascade is only as good as the model upon which it is based. Thus a focus of the experiments presented here was to test a range of manipulations for stimuli and photoreceptors that did not contribute to the model fits. We did this by fitting consensus models to recordings from one set of photoreceptors of a given type, and then testing the stimuli generated by these models on different photoreceptors from different retinas. Discrepancies between the predicted and measured responses to these new stimuli reflect both failures of the model inversion and cell-to-cell variability. Systematic differences between predicted and measured responses were relatively small, and most of the differences appear to reflect variability between cells (e.g. due to cell-to-cell differences in sensitivity).

Focus on phototransduction

Our model describes the relationship between light inputs and phototransduction currents. We restricted our consideration to photocurrents rather than photoreceptor synaptic output because the latter is shaped by several mechanisms (e.g. electrical coupling between photoreceptors (Copenhagen and Owen, 1976; Detwiler and Hodgkin, 1979; Schwartz, 1976) and horizontal feedback (Baylor et al., 1971)) that are less well understood than phototransduction. A similar approach could be extended beyond phototransduction to include electrical properties of the photoreceptor inner segment and the output synapses when additional quantitative information is available. The models we develop here will be important in such future work, as they will allow isolation of mechanisms operating post-transduction (e.g. in Figure 10C).

Photopigment bleaching

The model that we used does not consider photopigment bleaching and regeneration via the pigment epithelium. We omitted these aspects of photoreceptor function because the experiments to which we fit and test the model require that we remove the pigment epithelium. We restricted total light exposure in these experiments to minimize bleaching, and this limits the range of light levels over which the model is applicable. This is less of a concern for rods, where rod saturation occurs for light levels at which a small fraction of the rod pigment is bleached, and correspondingly our models cover most or all of the range of rod signals. For cones, the model is limited to light levels < 50,000 R*/s (similar to typical indoor lighting conditions).

Speeding vs slowing responses

Not all manipulations of the photoreceptor responses are equally easily achieved. Speeding responses is particularly difficult because it requires increasing the amplitude of high temporal frequencies, and the ability to do that is limited by the requirement that light intensities not assume negative values. This means that low contrast stimuli can be sped more than high contrast stimuli.

Calibrations

The accuracy of the model predictions depends critically on accurate light calibrations. We have detailed our procedure in Methods and supplied our calibration code (https://github.com/chrischen2/photoreceptorLinearization). With consistent calibrations, the approaches we use here should translate directly across laboratories and to in vivo physiological or perceptual studies

Applications

“Front end” for encoding models for downstream responses and perception

Two broad classes of encoding models have been used to describe responses of retinal ganglion cells and cells in downstream visual areas. Empirical models take light stimuli as inputs and convert these, through a series of linear and nonlinear elements, to predicted responses (Chichilnisky, 2001; Pillow et al., 2008). Time-dependent nonlinearities are particularly hard to capture in such models, and few existing models account for them. Computational neural network (CNN) models similarly take light input and convert it, via the learned CNN weights, to predicted neural responses (McIntosh et al., 2016; Turner et al., 2019). Like empirical models, time-dependent nonlinearities such as adaptation are generally not well described by such models. As a result, these models work best when stimulus parameters such as mean light intensity do not change.

The photoreceptor models described here convert light inputs to photoreceptor responses, capturing time-dependent nonlinearities in this process. Empirical or CNN models could then be used to describe the conversion of photoreceptor signals to downstream neural responses. This architecture - a phototransduction model front end followed by an empirical model - should decrease the demand placed on the empirical model and improve the ability to capture responses to stimuli that strongly engage photoreceptor adaptation. Hybrid models of this type indeed show improved performance in predicting retinal ganglion cell responses across light levels (Idrees et al., 2024). Models for visual perception could similarly incorporate a photoreceptor front end, and by doing so directly test which aspects of perception can be explained by phototransduction and which are due to downstream processing.

“Back end” to decoding models

Approaches to decode neural responses and estimate stimulus properties could similarly benefit from the photoreceptor models. Current decoding approaches either empirically fit response-stimulus relationships or invert encoding models to compute the likelihood of particular stimuli given the neural response (Bialek et al., 1991; Brackbill et al., 2020; Wu et al., 2023). In either case, decoding in the context of stimuli that strongly engage time-dependent nonlinearities has proven difficult.

Incorporating the inverse photoreceptor model into decoding approaches should improve decoding performance. Specifically, existing decoding approaches could be used to estimate the photoreceptor signal from downstream neural responses (e.g. those of retinal ganglion cells), and then the inverse phototransduction model would convert the estimated photoreceptor signals to estimated stimuli. This again could decrease the demand placed on the decoding model. For example, under conditions in which much of retinal adaptation is largely accounted for by adaptation in the photoreceptors, the model used to estimate photoreceptor responses from retinal ganglion cell responses would not need to incorporate adaptive nonlinearities. lnstead, this model could be fixed even as mean light levels change, and the inverse photoreceptor model would account for photoreceptor adaptation. A ganglion cell→photoreceptor model should be simpler to fit and hence perform better than a full ganglion cell→stimulus model.

Complementing genetic approaches

Genetic manipulations provide another approach to alter photoreceptor responses and characterize the impact on responses of downstream visual neurons or behavior. Such approaches have the advantage of associating specific molecular components with alterations in vision - e.g. the role of arrestin (Burns et al., 2006; Xu et al., 1997) and rhodopsin kinase (Arshavsky, 2002; Wilden et al., 1986; Zhao et al., 1995) in specific forms of stationary night blindness (Dryja, 2000; Zeitz et al., 2014). But interpretation of these genetic manipulations can be limited by compensatory changes.

The approach that we introduce here complements genetic manipulations in at least two ways. First, it would be used in conjunction with genetic manipulations - e.g. to restore normal kinetics in photoreceptors in which the genetic alteration changes response kinetics. Second, the stimulus design approach provides an alternative that provides less detailed mechanistic information, but which allows more specific functional manipulations to be made and does not require genetic access.

Manipulations of phototransduction currents and identifying nonlinear circuit properties

The models that we introduce here provide a new tool to causally test the impact of alterations in photoreceptor responses on downstream responses and perception. For example, asymmetries in how increments and decrements in light intensity are processed have been well studied, including responses in retina, cortex and perception (Bowen et al., 1989; Lu and Sperling, 2012). Parallel On and Off visual pathways are initiated at the cone output synapse, and asymmetries in downstream responses or perception are often attributed to asymmetries in the On and Off retinal circuits conveying photoreceptor signals to the ganglion cells (e.g. (Stockman et al., 2014)). The implicit assumption is that photoreceptor inputs to On and Off retinal circuits are symmetrical. While this is likely the case for low contrast stimuli, adaptation in phototransduction means that high contrast stimuli will often not produce symmetric input to those circuits, and On/Off asymmetries may originate at least in part from asymmetric photoreceptor responses to light increments and decrements (Angueyra et al., 2022; Clark et al., 2013; Endeman and Kamermans, 2010). The ability to invert the phototransduction model permits the design of stimuli that minimize such asymmetries in the photoreceptor responses to test how they contribute to On/Off differences in downstream responses (e.g. (Yu et al., 2022)).

More generally, the ability to shape photoreceptor responses in predictable ways provides a needed tool to isolate the effects of photoreceptor and post-photoreceptor circuits to shaping responses of downstream neurons and perception. Nonlinear circuit properties are a particular challenge since their impact depends on the input stimulus, yet they are also particularly important since they account for the bulk of interesting computation. Advances in modeling, particularly CNN-based approaches, can be used to fit circuit outputs and reveal how specific computations may be implemented. But such approaches rarely identify a unique explanation. Another approach is to record from each of the relevant circuit elements. This is technically challenging. Manipulating photoreceptor responses provides another tool, as illustrated in Figure 10. A key feature of this tool is the ability of our photoreceptor models to generalize across stimuli, including accounting for nonlinear properties of the photoreceptor responses to stimuli that did not contribute to the model directly.

Methods

Recordings

We performed electrophysiological recordings from primate (Macaca fascicularis, nemestrina and mulatta, either sex, 2-20 years) and mouse retina (C57/BL6 or sv-129, either sex, 1-12 months) in accordance with the University of Washington Institutional Animal Care and Use Committee. Primate retinas were obtained through the Tissue Distribution Program of the University of Washington Regional Primate Research Center. Primate recordings were all at > 20° eccentricity.

Rod responses were recorded with suction electrodes (Field and Rieke, 2002). These recordings were sufficiently stable that we could record responses to all three test stimuli in Figure 1B. Periodic standardized flashes tested for changes in kinetics or amplitude of the light response, and recordings were terminated if such changes were apparent. Most recordings were stable for at least 10-15 min before changes in kinetics were apparent.

Cones were recorded with whole-cell patch clamp techniques in slices (mouse; (Ingram et al., 2019)) and whole-mount preparations (primate; (Angueyra and Rieke, 2013)). We focused on M and S cones in mouse and L and M cones in primate. Cone responses in these recordings run down quickly due to intracellular dialysis of the cell; hence responses to only one of the stimuli in Figure 1B were recorded from each cone. All data reported here were collected within 2-3 min of patch rupture and the onset of intracellular dialysis.

Light calibrations

Optical power was measured at the preparation with a calibrated power meter (Graseby). Power readings were converted to isomerizations per second (R*/sec) using the photoreceptor spectral sensitivities, the LED spectral outputs, and collecting areas of rods and cones (1 and 0.37 μm2 for primates, 0.5 and 0.2 μm2 for mouse). This conversion proceeded in three steps: (1) the total calibrated power was converted to a power density (W/μm^2) using the size of the illuminated area in the microscope image plane; (2) LED emission spectra were scaled such that their integrals matched the calibrated power density; (3) the resulting calibrated LED emission spectra were converted to R*/sec by taking the dot product with the photoreceptor spectral sensitivities (Baylor et al., 1987; Nikonov et al., 2006).

Light stimuli from the LEDs were focused on the preparation via the microscope condenser or objective. Stimuli uniformly illuminated a 600 μm diameter spot. LED outputs were carefully checked for linearity. For rod recordings, we used an LED with peak output at 470 nm, near the peak of the rod spectral sensitivity. For cone recordings, we used an LED with peak output at 405 nm, which produced near-equal (within 10%) activation of L and M cones in primate and S and M cones in mouse. Figure 10 used an OLED computer monitor (eMagin) to deliver stimuli; the monitor outputs were calibrated as above for the LEDs.

Forward model

The phototransduction cascade was modeled with the set of differential equations illustrated in Figure 1A. This follows previous modeling work in both rods and cones (Angueyra et al., 2022; Nikonov et al., 2000; Pugh and Lamb, 1993; Rieke and Baylor, 1998a; Younger et al., 1996). Our model has 11 parameters. In the first step of the model, the light stimulus (Stim) activates opsin molecules, converting inactive opsin (R) to active opsin (R*). Active opsin decays with a rate constant σ (Figure 1A, Equation 1). Next, phosphodiesterase (PDE) molecules are activated by active opsin molecules via the G-protein transducin; the resulting PDE activity (P) decays with a rate constant Φ (Figure 1A, Equation 2). We assume that the delay caused by transducin activation is negligible and hence omit this step for simplicity. The opsin decay rate σ and the PDE decay rate Φ were interchangeable (Pugh and Lamb, 1993), and the model output depended only on the smaller of these (i.e. the slower process). Hence σ and Φ were constrained to be equal.

The concentration of cGMP in the outer segment (G) is determined by the balance of cGMP hydrolysis mediated by PDE and synthesis (S) through guanylate cyclase (GC) (Figure 1A, Equation 3). For physiological conditions, the outer segment current (I) through cGMP-gated channels depends on a power of the cGMP concentration (Figure 1A, Equation 4; (Rieke and Baylor, 1996)). Calcium enters the outer segment through open cGMP channels and is removed by an exchange protein, with a rate constant β (Figure 1A, Equation 5). The calcium concentration (C) regulates the rate of cGMP synthesis (S); this dependence is modeled as a Hill curve (Figure 1A, Equation 6).

Two model parameters were fixed by steady state conditions and prior measurements. First, the constant q that relates current and changes in calcium can be expressed in terms of the dark calcium concentration, dark current, and rate constant β for calcium extrusion:

Second, the maximal cyclase rate (Smax) can be written in terms of Φ, η, KGC, h, and the dark calcium (CD) and cGMP concentrations (GD):

We fixed the constants k and n that determine the relation between cGMP and current at values measured for rods (e.g. (Rieke and Baylor, 1996)); modest changes in these parameters produced little or no change in model performance due to compensation by other model parameters. Using these constants, the cGMP concentration in darkness was estimated from the measured dark current using Eq. 3 (e.g. the response to a saturating light flash) for each recorded cell. We also fixed h for both rods and cones and β for rods based on previous measurements (Field and Rieke, 2002; Koutalos et al., 1995a; Rieke and Baylor, 1996). We assumed a dark calcium concentration (CD) of 1 μM (Gray-Keller and Detwiler, 1994; Sampath et al., 1999); the model was insensitive to this value, as the dependence of the synthesis rate on calcium (KGC) compensated for any change.

This left a total of 4 free parameters for rod models (γ, η, σ, KGC), and 5 for cone models (γ, η, β, σ, KGC). These parameters were determined by minimizing the mean-squared difference between the measured responses and model predictions while the remaining model parameters were held fixed. Some combinations of model parameters can trade for each other, resulting in similar model performance (e.g. CD and KGC). Our goal here was to identify a model that captures photoreceptor responses across a range of stimuli, rather than a model in which the fit parameters were predictions of the underlying biochemical parameters.

We identified optimal values for the free parameters (see Table 1) using the MATLAB fminsearch routine, which employs the Nelder-Mead simplex algorithm. The time step in these calculations was 0.1 ms for cones and 1 ms for rods. The objective function minimized the mean squared error (MSE) between the measured and predicted responses. The goodness of fit was assessed using the fraction of variance explained, calculated as 1 - (SSE / SST), where SSE is the sum of squared errors between the model and the data, and SST is the total sum of squares (variance) of the data. The model’s performance was validated using independent data sets.

Parameters of consensus models were determined by simultaneously fitting measured responses across all recorded cells of a given type. To evaluate the model’s ability to generalize across different cells and stimuli, we fixed all parameters except for GD and γ. GD was determined from the measured dark current (ID) for each individual photoreceptor. γ was allowed to vary between cells to account for differences in sensitivity, while the remaining parameters were constrained to be identical. For example, consensus parameters for a set of 5 cones would have common free parameters of η, β, σ, KGC and 5 individual γs (one for each cell). Consensus parameters are given in Table 1. Fitted γ values were 8 ± 3 (mean ± SD) for primate cones, 3 ± 1 for mouse cones, 5.3 ± 0.7 for primate rods and 7.4 ± 2.4 for mouse rods.

We evaluated the sensitivity of the model outputs to each fit parameter by identifying parameter ranges that doubled the mean-square error (MSE) of the fits. Table 2 includes these parameter ranges. We also tested sensitivity to combinations of parameters using sloppy modeling (see further details below and (Brown and Sethna, 2003)). This approach identifies specific combinations of parameters that have maximal and minimal impact on model error - in our case the mean-square error. For cone models, parameters controlling the calcium feedback (β and KGC) could compensate for each other to have relatively minor effects on fit quality. Fit quality was particularly sensitive to changes in the parameters controlling the PDE activation (σ and η). Rod models were relatively insensitive to changes in σ and η and were particularly sensitive to changes in KGC and γ.

Parameter variations that double the mean-squared-error of the fits. Error was measured while a single parameter was systematically varied; other parameters were held fixed at the consensus model values.

Sloppy Modeling

Sloppy modeling enables us to understand how a model’s cost function depends on the model parameters individually and in combination. We started by computing the Hessian matrix H at a location in parameter space given by the consensus fit parameters

Here θi and θj represent different parameters and C is the cost function. The Hessian provides a measure of the curvature of the cost function around the best-fit model parameters. We numerically approximate the Hessian by taking incremental changes in pairs of parameters around the best fit and measuring the changes in MSE. The eigenvectors of the Hessian matrix identify directions in parameter space in which the cost function changes slowly or rapidly.

Linear Model

The linear models used to generate the target response were obtained by fitting a parameterized linear filter to the response of the full model to low-contrast Gaussian noise at a specific mean light level. The linear filter is defined as:

Where α is a scaling factor, and τR and τD are the rising and decay time constants. This model provides good fits to measured low-contrast responses of rods and cones (Angueyra and Rieke 2013). Photoreceptor responses to low contrast stimuli depend linearly on the stimulus (Hass et al. 2015), and hence this approach ensured that the linear model matched the full model for such stimuli. The linear filter was convolved directly with the light stimulus to obtain a linear estimate of the responses.

Inverse model

The differential equations that comprise our model consist of several time-independent nonlinearities (Eqs. 4 and 6 in Figure 1), several linear differential equations (Eqs. 1, 2, 5 in Figure 1), and one nonlinear differential equation (Eq. 3 in Figure 1). The time-independent nonlinearities can be inverted directly as a look-up table - e.g. the cGMP depends directly on the current as G = (I/k)1/n. The linear differential equations can be solved directly by deconvolution - e.g. the Fourier transform of the calcium concentration from Eq. 5 in Figure 1 is obtained directly from the Fourier transform of the measured current (I(f)) as

where f is the temporal frequency. The nonlinear differential equation (Eq. 3 in Figure 1) can be reexpressed as an equation for the PDE activity P in terms of the cGMP concentration G and its derivative dG/dt as

This equation can be solved given G (obtained from the current via Eq. 4 in Figure 1), and S (obtained from the calcium concentration C(t) and Eq. 6 in Figure 1). These steps in the inversion process are illustrated in Figure 3.

This process results in an exact mapping between responses and stimuli - i.e. every response I(t) has a corresponding unique stimulus Stim(t). In practice, noise in the measured responses can corrupt the estimated stimuli. This is a general problem in deconvolution of noisy data – e.g. in microscopy. As indicated in the text, we controlled this noise when necessary by imposing that the power spectrum of our estimated stimulus be equal to the power spectrum of the true stimulus. This constraint was imposed by generating the estimated stimulus as described above, and then reweighting the power spectrum of the estimated stimulus to match the power spectrum of the true stimulus. This constraint was used in generating the estimated stimuli in Figure 4.

Predictably altering photoreceptor responses

Model inversion as described above allowed us to identify stimuli that would produce a desired phototransduction current (which we refer to as a target response). For the applications in Figures 6-10 we generated this response using the linear model of the transduction cascade described above (Eq. 7). We passed an initial stimulus (e.g. a sinusoid in Figures 6 and 7) through the linear model to generate the target response. The model inversion process for the full (i.e. nonlinear) phototransduction cascade model then determined the stimulus which, when delivered to the full model, would produce the linear target response. In particular, this allowed us to identify stimuli that negated the impact of adaptation on the photoreceptor responses.

Acknowledgements

We thank Shellee Cunnington for excellent technical support. We are grateful for assistance from the WaNPRC staff, especially Chris English, for access to primate tissue. Shane Boniec, Greg Horwitz and Saad Idrees provided invaluable feedback on a draft of the manuscript. This work was supported by NIH grant EY028542 (FR) and the Air Force Office of Scientific Research under award number FA9550-21-1-0230 (special thanks to the AFOSR’s Cognitive & Computational Neuroscience Program).

Tests of the ability of the model to generalize across stimuli and across cells. A. Generalization across stimuli. Responses of primate and mouse rods were fit to the variable mean noise stimulus only and model performance was evaluated for responses to the flash family and flashes and steps, using the fraction of variance explained to evaluate the fits (see Figure 1). This generalization performance to stimuli not included in the fit (y-axis) is compared to performance when the model is fit to all stimuli (x-axis). Means ± SDs were 0.92 ± 0,05 (generalize) and 0.94 ± 0.03 (fit all) for 6 primate rods and 0.90 ± 0.04 and 0.92 ± 0.02 for 8 mouse rods. B. Generalization across cells. Models were fit to the variable mean noise stimulus. The performance of each cell was evaluated for a model in which all cells were fit (x-axis), or in which the test cell was excluded from the fit (y-axis). γ was allowed to vary for the held-out cell while other parameters were determined by fitting the remaining data as for the consensus models. Means ± SDs were 0.88 ± 0.04 (hold out) and 0.87 ± 0.05 (fit all) for 6 primate cones, 0.90 ± 0.02 and 0.91 ± 0.02 for 6 mouse cones, 0.96 ± 0.02 and 0.96 ± 0.02 for 6 primate rods, and 0.95 ± 0.02 and 0.95 ± 0.02 for 8 mouse rods.

Full models systematically outperform linear models. Linear models were generated by fitting low contrast responses to the full model for each cell. Full models and linear models were optimized with a single scale factor for each cell (γ for the full model and α for the linear model in Eq. 7). The fraction of variance explained was computed for the variable mean noise stimuli as in Figure 2. Means ± SDs were 0.63 ± 0.15 (linear) and 0.87 ± 0.05 (consensus) for 6 primate cones, 0.68 ± 0.04 and 0.91 ± 0.02 for 6 mouse cones, 0.86 ± 0.02 and 0.96 ± 0.02 for 6 primate rods and 0.86 ± 0.04 and 0.95 ± 0.02 for 8 mouse rods.