Phototransduction in vertebrate photoreceptor cells represents a paradigm of signaling pathways mediated by G-protein-coupled receptors (GPCRs), which share common modules linking the initiation of the cascade to the final response of the cell. In this work, we focused on the recovery phase of the visual photoresponse, which is comprised of several interacting mechanisms.
We employed current biochemical knowledge to investigate the response mechanisms of a comprehensive model of the visual phototransduction pathway. In particular, we have improved the model by implementing a more detailed representation of the recoverin (Rec)-mediated calcium feedback on rhodopsin kinase and including a dynamic arrestin (Arr) oligomerization mechanism. The model was successfully employed to investigate the rate limiting steps in the recovery of the rod photoreceptor cell after illumination. Simulation of experimental conditions in which the expression levels of rhodospin kinase (RK), of the regulator of the G-protein signaling (RGS), of Arr and of Rec were altered individually or in combination revealed severe kinetic constraints to the dynamics of the overall network.
Our simulations confirm that RGS-mediated effector shutdown is the rate-limiting step in the recovery of the photoreceptor and show that the dynamic formation and dissociation of Arr homodimers and homotetramers at different light intensities significantly affect the timing of rhodopsin shutdown. The transition of Arr from its oligomeric storage forms to its monomeric form serves to temper its availability in the functional state. Our results may explain the puzzling evidence that overexpressing RK does not influence the saturation time of rod cells at bright light stimuli. The approach presented here could be extended to the study of other GPCR signaling pathways.
Phototransduction is the biochemical process by which a visual stimulus initiates a neuronal response in the photoreceptor cells of the retina. The capture of a photon by a visual pigment molecule triggers a G-protein-coupled receptor signaling cascade that leads to the closure of cGMP-gated ion channels, which decreases intracellular Ca2+ and causes the hyper-polarization of the cell membrane. Changing Ca2+ concentrations in the 0.1 – 0.6 μM range contribute to at least three feedback processes, which serve to regulate tightly the recovery of a photoresponse . One particular means of Ca2+ feedback is exerted through recoverin (Rec), a protein which binds and regulates the action of rhodopsin kinase (RK). When Rec binds two Ca2+ ions, it undergoes a conformational change, exposing a myristoyl group . This “relaxed” (RecRCa) form prevails in the dark, at higher intracellular Ca2+ concentrations, and increases the protein’s affinity for the disk membrane due to an overall augmented hydrophobicity . In the relaxed form, Rec binds RK and prevents it from phosphorylating the activated photopigment rhodopsin (R*) [4–9]. As Ca2+ concentration decreases during a photoresponse, Rec reverts to its tense form (RecT), becoming more hydrophilic, and ceasing to bind RK. The free RK may then phosphorylate R*, leading to the deactivation of R* by arrestin (Arr) and the termination of the response.
One mechanism of the phototransduction cascade whose function remains unclear is the formation of homo-oligomers by Arr. It is known that in dark-adapted rod cells, Arr exists in a steady-state equilibrium of monomers, dimers and tetramers and that during light-adaptation the oligomeric forms dissociate, increasing the concentration of monomeric Arr available to quench R* [10–12]. It is generally thought that the dimers and tetramers together function as storage forms and exist to maintain the ideal concentration of available, monomeric Arr, however the kinetic implications of Arr oligomerization have not been fully explored. Supplementary to biochemical and transgenic animal assays, detailed modeling of the phototransduction system could provide key insights into the dynamic effects of this and other mechanisms.
Mathematical simulations of biochemical networks are essential for gaining a proper understanding of the complex dynamics underlying them. As perhaps one of the best-understood GPCR-mediated signaling pathways, the phototransduction system has seen a rich history of modeling. This has led to comprehensive models of the amphibian phototransduction cascade, which include a large portion of the known reactions that occur during phototransduction [1, 13–15]. Despite the intrinsic limits brought in by the well-stirred approximation, such comprehensive models have proven realistic and widely useful. For instance, systems-level analysis of the dynamics underlying phototransduction has been instrumental in predicting novel functions of signaling proteins in the network including rhodopsin (ref. ) and in unraveling some molecular mechanisms associated to disease (ref. ). Different modeling approaches account for spatial dynamics of the processes occurring in specific cell compartments [17, 18]. While they provide better physical description of the processes, of great relevance to single photon response dynamics, to date these approaches have been applied to limited subsets of reactions.
Recently, a systems biology approach was used to build a comprehensive model of phototransduction in rod cells based on most of the known biochemical information on the signaling processes . The model relied on the quantitative data collected over 40 years of biochemical and biophysical research on amphibian rod photoreceptors, and included all of the modular components of the signaling cascade, namely: a) the activation of the receptor rhodopsin by light stimuli of different duration and intensity; b) the signal amplification steps occurring upon transferring the information from the receptor to the G protein (transducin) and from the G protein to the effector phosphodiesterase 6 (PDE); c) the Ca2+ −mediated feedback mechanisms on the regulation of processes occurring both on discs and on the plasma membrane; and d) the deactivation of both receptor and effector for normal recovery of dark-adaptation cell conditions.
The model successfully reproduced a number of experimental data collected on photoreceptors from different species, stimulated by light ranging over five order of magnitude in intensity [15, 19]. Interestingly, while the quantitative comparison between simulations and experimental data showed striking consistency for amphibian cells, the model proved also able to predict with great accuracy the qualitative dynamics observed in mouse photoreceptors . This latter finding proved the great potential of the system-level analysis, which apparently is able to capture the evolutionarily important topological and dynamic constraints in the same signal transduction pathways in different species. Moreover, the model reproduced and predicted dynamic behaviors observed with the widely used approach of producing transgenic mice with genetic manipulations aimed at investigating the specific role of genes and their products in regulating the cascade. The development of this model has resulted in a framework in which the phototransduction process is broken down into its fundamental reactions, which are then modeled largely according to the law of mass action. This resulting model is effectively modular in design, allowing relatively isolated modifications to a mechanism without great disruptions to the others [20, 21].
We have taken advantage of this modular structure in order to incorporate the latest knowledge on the kinetics of Ca2+-mediated Rec feedback on R* shutdown into a new iteration of this model. We then further extended it to include a dynamic Arr oligomerization mechanism. The modified model provides an unprecedented ability to probe this mechanism’s influence on the dynamics of R* shutdown. Simulations performed with the present model allowed full reproducibility of a number of diverse experimental results, including experiments performed with photoreceptors from animals carrying genetic modifications in the key components of the shutdown machinery.
Using the present model, we have replicated several published experiments, not only to verify that the model can accurately reproduce the experimental results but also to summarize and explore the current knowledge of phototransduction recovery mechanisms and the contributions of these two mechanisms to the process. We found that the model gives qualitatively accurate reproduction of experimental data (accounting for species differences between the amphibian model and the experiments, which were largely performed on mice). In particular, the results suggest an important modulatory role for the process of Arr oligomerization. Lastly, we note that to perform such extensive experimentation with animals in a laboratory would be extremely expensive and time-consuming, pointing to a beneficial role of large-scale system modeling in the experimentation process.
A comprehensive, system-level description of the signaling cascade
We have improved and extended the model of Dell’Orco et al. (ref ), making significant steps forward in the comprehensive, biochemistry-based description of the phototransduction cascade. We herein focused on the still largely unclear mechanisms constituting rate-limiting steps in the kinetics of cell recovery after illumination. First, we reconsidered the affinity relationships of RK and Arr with phosphorylated R*. Previously, these affinities varied exponentially with the number of phosphates attached to R*. Based on experimental data and considerations for goodness-of-fit versus parameter sensitivity, we replaced RK’s exponential decrease in affinity and Arr’s exponential increase in affinity with respective linear relationships (see Methods).
We next replaced a representation of the Rec-mediated Ca2+ feedback on RK based on a quasi-steady state assumption with a realistic process consisting of three fundamental reactions. This permitted the removal of another exponential parameter (w in the model of Dell’Orco et al. (ref )) as well as the explicit modeling of Rec as a significant Ca2+ buffer. Finally, we implemented a dynamic Arr homo-oligomerization mechanism, in which Arr self-associates to form homodimers and homotetramers. The implementation of the Rec-mediated calcium feedback on RK and the dynamic Arr oligomerization mechanism was realized as illustrated in the Methods session. The new model thus resulted in an extended network of interactions (Figure 1), which we used to explore, by numerical simulation, currently inaccessible molecular scenarios, in which the expression level of key proteins were varied in order to gain insights into the regulatory mechanisms acting in the recovery phase of the phototransduction cascade.
We verified the present model by recapitulating simulation results published previously by Dell’Orco et al. , in which illumination intensities extended over five orders of magnitude. The new model features were found not to alter greatly the system’s dynamics for a wide range of experimental conditions. Responses to series of brief flashes, leading from few rhodopsin photoisomerizations up to saturating levels, or to series of prolonged pulses of light of increasing intensity both closely matched experimental data and previously simulated results (Figure 2A & B). Notably, when exposed to prolonged pulses (60s) of light of stronger intensity, this model shows decreased time-to-plateau recovery after the rising phase at all light intensities. Indeed, inspection of the experimental data reveals a similar pattern, with slower recovery transients preceding the plateau .
The performance of the modified model in reproducing light adaptation rather than dark-adapted behaviors was also checked (Additional file 1: Figure S1). Both reduced light sensitivity and accelerated recovery when a saturating stimulus is delivered in the presence of a non-saturating background of light were accurately reproduced by the model (Additional file 1: Figure S1).
The modified model was also tested in its capability of reproducing photoresponses by genetically modified animals (See Methods and Additional file 1). Simulations generally reproduced the experimental data well. However, unlike the previous implementation of the dynamic model, the simulation of RK overexpression more accurately reproduced experimental results . Indeed, in the simulated results of Dell’Orco et al. , RK overexpression resulted in shorter saturation times (Tsat; the duration of time that a response remains at at least 90% of the maximal amplitude after a saturating flash stimulus) after bright flashes (Figure 3), which differs from experimental results; it has previously been shown that Tsat of rods overexpressing RK is not significantly different than that of wild-type (WT) rods . The new model largely resolves this discrepancy.
The model was then employed to elucidate the complex deactivation dynamics of phototransduction. Many experiments have been performed previously to explore these processes through the use of transgenic animals or other means of altering the expression levels of key proteins to determine their roles. Using the present model, we recapitulated several of these experiments. Because much current research on the phototransduction system is performed using transgenic mice, we were precluded from performing quantitative comparisons against the amphibian-based model. However, as previously demonstrated by Dell’Orco et al. , the model with its set of kinetic constraints is capable of generating results that are qualitatively comparable to experimental data, despite species differences, across a wide range of experimental conditions (Figures 2, 3 and Additional file 1: Figure S2-S3).
We will now present in detail the most salient features of the simulated photoresponses obtained with the extended model and make direct comparisons with the known experimental data.
Increased RK availability accelerates recovery from saturation
RK plays a major role in R* shutdown by regulating its phosphorylation, which finally determines the capping by Arr, therefore it is not surprising that much experimental effort in recent years has been made to elucidate the quantitative effects of alterations in the expression level of RK on the overall cascade. In this respect, in a recent work Sakurai et al. generated lines of transgenic mice with different expression levels of RK (two-and three-fold overexpression compared to the wild type) . For mice overexpressing RK, they found results that conflicted with those of Krispel et al. (ref. ), who had produced similar transgenic mice and found that RK overexpression had neither an effect on Tsat for saturating flashes nor on the recovery time constant (τrec; the time constant of a single exponential function fit to the second half of the recovery phase of a non-saturating response). Sakurai et al. found that, for their animals, while the two- and three-fold overexpression of RK did not have an effect on Tsat, in agreement with the previous findings, it did, in fact, lead to a reduction in τrec for non-saturating responses. When we simulated their experiments, we found the same pattern of a decrease in τrec with increasing RK quantity (Figure 4A′ & B′). When the model was tested with saturating flashes, the model showed a slight decrease in Tsat with RK overexpression compared to wild-type (WT) for similar flash intensities, in contrast to the experimental results (Figure 4C′).
Furthermore, Sakurai et al. found that the dominant time constant of recovery from a saturating response (τD; measured as the slope of Tsat over logarithmically increasing stimulus intensities, visualized in a so-called Pepperberg plot as in Figure 4B and B′ ), did not significantly change with RK expression levels. The interpretation given by the authors for this result is that the measurement of τD occurs after RK-mediated rhodopsin (R*) shutoff has been completed and thus only reflects the rate-limiting step of RGS-mediated effector shutdown. The measurement of τrec, on the other hand, may still be influenced by the rate of R* shutoff [24, 26] This stands in contrast to the results of Krispel et al., who found that τD slightly increased with RK two-fold overexpression . In simulations, we found that τD did not vary greatly with RK expression (Figure 4D, Table 1). Since τD is relatively constant with increasing RK expression, the rate-limiting step of recovery is unaffected by RK activity in the model.
τDValues for simulated experiments
0.3x RK UX
0.4x RK UX
2x RK OX
3x RK OX
4x Rec OX
0.2x RGS UX
2x RGS OX
4x RGS OX
6x RGS OX
0.4x RK/6x RGS
3x RK/0.2x RGS
τD values, or the change in saturation time with logarithmically increasing stimulus intensities, for the various experiments simulated herein. Normalization is done against the wild-type value. (WT = wild-type, UX = underexpression, OX = overexpression, KO = knock-out). Only RGS mutants have significant effects on τD.
Intuitively, an overabundance of available RK due to the complete absence of Rec inhibition acting on it would be expected to show similar dynamic behavior to an overabundance of RK due to its overexpression. However, experiments using Rec knock-out mice prove otherwise . When Rec is absent, phototransduction signaling shows decreased sensitivity and shorter saturation times (Figure 5A & B). Unlike the case of RK overexpression, Rec knockout animals not only have shorter Tsat than WT for similar stimulus intensities, but they also show a slight reduction in τD. Moreover, the shape of the decreasing phase of the photoresponses appears to differ from that of the control (compare Figure 5A and B), as in the case of Rec−/− three phases are clearly visible at each light intensity (Figure 5B). When these conditions were simulated, the decreasing photoresponse remained substantially biphasic, moreover we saw only a moderate shift in Tsat, similar to the result of RK overexpression, and no significant change in τD (Figure 5C; Table 1). Thus, whereas experimentally, mice lacking Rec required approximately 9.7-fold more light to reach similar saturation times as WT (ref. ), the model required only approximately 1.6-fold more light to reach WT saturation times. Therefore, the model in fact fails to disentangle the effects of RK overabundance due to RK overexpression and Rec knockout mutations.
Decreased RK availability slows recovery
Decreasing the expression level of RK allows one to directly test whether R* shutdown is the rate-limiting step for the recovery of the cascade. Moreover, measuring changes in Tsat and τD would allow one to assess quantitatively the effects over the whole illumination range. Indeed, it was found in transgenic mice which underexpress RK (0.4x) that the response dynamics are slowed significantly [24, 28]. Reduced RK expression resulted in a slightly increased τD as well as a consistently increased τrec. Similarly, when amphibian rod outer segments were dialyzed with an excess of Rec, response recovery was found to be significantly prolonged . Simulating both of these cases resulted in response curves in agreement with the experimental data. For example, the recovery phase slowed down both for 40% underexpression of RK (compare Figure 6A and B) and for four-fold overexpression of Rec (compare Figure 7A and B). RK underexpression resulted in a negligible change in τD relative to WT, while Rec overexpression caused a large increase in this measurement (Table 1; note that the difference in timescale between simulated and experimental values arises from species differences).
To determine if R*-shutdown was indeed the rate-limiting process in rods underexpressing RK, Chen et al. produced double mutants that also overexpressed RGS (Figure 6C) . If R*-shutdown were truly rate-limiting, then increased RGS activity should not have a significant effect on the response dynamics. However, it was found that the double mutants have faster shutdown dynamics relative to both WT and RK underexpressing rods (Figure 6D), though they are slower relative to rods which only overexpress RGS (Figure 6C). These results could be replicated in our simulations (Figure 6G and H) and further confirm that RGS-mediated effector shutdown is the rate-limiting step in phototransduction deactivation, as further demonstrated by the following simulations.
Effector shutdown is rate-limiting
According to current knowledge, RGS-mediated shutdown of the effector is the rate-limiting step of response recovery. The rate-limiting step is expected to primarily determine τD, or in other words, how Tsat changes with increasing stimulus intensities. Thus, decreasing RGS activity is expected to greatly affect Tsat and τD. Furthermore, since τD is ultimately affected by the rate-limiting step, complementing a scenario of lowered RGS activity with increased RK should not have any additional effect on τD. In fact, when such a hypothetical experiment was simulated, we saw that the dominant time constants were approximately equivalent for rods that only underexpress RGS and double-mutants which additionally overexpress RK, (Table 1). However, we saw a consistent shift in Tsat, for the double-mutant, which required approximately one half log step stronger stimuli to achieve the same Tsat,as the single-mutant (Figure 8).
It has recently been reported that, as the key protein of the rate-limiting step, RGS’s expression level is extremely important to the recovery dynamics of the phototransduction system . We took advantage of the model to verify these results. Despite the difference in species between this experiment and our model, the simulated results were qualitatively in agreement with the experimental data (Additional file 1: Figure S3). τD varies greatly with RGS expression levels in a manner consistent with the published experimental data. It should be noted that a simple model of RGS activity was included by Burns and Pugh to quantitatively explain RGS dynamics relying on Michaelis-Menten kinetics . We show with the present model that we can arrive at the same conclusions by representing the system with only mass-action kinetics of the fundamental reactions involved therein (Additional file 1: Figure S3).
Arrestin concentration does not affect recovery dynamics
If Arr availability due to its oligomerization is important to recovery dynamics, then one might suppose that its concentration also may have a significant influence. However, Gross and Burns showed that underexpressing Arr (0.5x) in mice does not perturb the kinetics of photoresponses over a broad range of light stimuli (~5 - 5,000 R*/flash) . We simulated both the 0.5x under- and 2x overexpression of Arr. In line with the experimental results , underexpressing Arr does not perturb the photoresponse (compare Figure 9A and 9B). Moreover, our model predicts that doubling Arr levels would also be ineffective (Figure 9C). Thus, it would appear that Arr’s affinity for itself works to maintain an ideal concentration of the monomeric form, even under conditions that would change its expression level.
A unique advantage of the model proposed here is that it allows the observation of the time-evolution of the concentrations of individual molecular species in the signaling network, thus revealing interconnections between the underlying dynamics. We traced the temporal evolution of selected molecular quantities, following a saturating flash stimulus leading to 118,000 R*, namely: i) R* with zero to three phosphorylations (Figure 10A); ii) R* with four to six phosphorylations (i.e. the maximum number in the model; Figure 10B); iii) R*-Arr complexes in which R* has been phosphorylated one to three times (Figure 10C); and R*-Arr complexes in which R* has been phosphorylated four to six times (Figure 10D). We also simulated the change in molecular quantities of Arr in monomeric and oligomeric form (Figure 10E and F, respectively). Results compared WT (black solid lines in Figure 10) and 2.4-fold RK overexpression conditions (black, dashed lines in Figure 10). Interestingly, it can be seen that most R* molecules under these conditions reach four to six phosphorylations before binding Arr. It should be noted, however, that it has been shown that three phosphorylation events are sufficient for monomeric Arr to bind R* (ref. ); indeed, we did find that some R* was bound to Arr with less than four phosphorylations (Figure 10B). When RK is overexpressed, more R* molecules unsurprisingly attain at least four phosphorylations (Figure 10C) and they are bound by Arr slightly more quickly (Figure 10D). Under both WT and RK overexpression conditions, Arr binding with less than four phosphorylations is completed within approximately one second (Figure 10B), while R* molecules with four or more phosphorylations require approximately three seconds to be completely bound by Arr (Figure 10D).
Observing the time evolution of Arr molecules in monomeric or oligomeric form supports the notion that Arr oligomers serve to maintain the amount of monomers that are available to bind with R*. After an initial decrease, the quantity of available Arr monomers quickly recovers to nearly the dark steady-state number (Figure 10E). The number of Arr molecules in an oligomeric form, on the other hand, steadily decreases and does not recover within the observed time frame, owing to the slow release of Arr from R* (Figure 10F).
We next investigated the kinetic effect of Arr oligomerization on R* shutdown by repeating the simulations with the rate of Arr self-association reduced to 0 s-1. With Arr oligomerization disabled, we found that a much larger percentage of R* binds Arr before four phosphorylations have been acquired (red traces, Figure 10A, B). Furthermore, Arr binding to R*, regardless of the level of R* phosphorylation, proceeds at a faster rate such that all R* is bound to Arr within approximately one second (red traces, Figure 10B, D). The rate is even faster under RK overexpression conditions. Thus, it is likely that Arr oligomerization is the primary mechanism that can account for the failure of the model of Dell’Orco et al.  to reproduce the lack of an effect of RK overexpression on the time spent in saturation after a bright stimulus.
The modeling approach used here differs from the classical approaches to modeling the phototransduction pathway by its focus on implementing each reaction in a fundamental form described by mass-action kinetics rather than via more abstract kinetics. By basing the model on the principal reactions, preferring linear combinations of parameters, realistic dynamics arise simply through the fundamental nature of the complex molecular interactions rather than through any introduced mathematical complexity. Thus, the model uses a “bottom-up” modeling approach, in which the process is understood through the emergent dynamics of a complex system, as opposed to a “top-down” approach, in which the observed dynamics are distilled down to mathematical relationships which can accurately reproduce them.
We have made significant steps forward in the comprehensive, biochemistry-based description of the phototransduction cascade, focusing on the still largely unclear mechanisms constituting limiting steps in the kinetics of cell recovery after illumination, by describing the Rec-RK interaction in the model with a representation quantitatively accounting for the experimentally determined kinetics and then by extending it to include Arr oligomerization (see Methods). The final model represents perhaps the most comprehensive molecular model of visual phototransduction to date.
Importantly, the new Rec-RK implementation did not significantly alter the dynamics of the model, while its simpler implementation compared to the previous version of the model allowed easier manipulation of these two proteins to clarify their own influence on the dynamics. Meanwhile, the novel Arr oligomerization mechanism allowed us to better understand the discrepancy between simulated and experimental responses of RK overexpression. Our simulated results confirm that the RGS-mediated effector shutdown is the rate-limiting step of phototransduction deactivation, as has been previously asserted. Furthermore, our approach allowed us to investigate this in finer detail than what may be feasible in vitro or in vivo.
We find that for the much faster process of R* shutdown, Arr is rate-limiting due to its homo-oligomerization. This has been alluded to in previous experiments which showed that saturation time after exposure to a bright stimulus is unaffected by an increase in RK activity due to its overexpression [23, 24]. With our model, we could demonstrate that this effect is at least partially mediated by the kinetics of Arr oligomer dissociation. On the other hand, our model fails to capture the fundamental difference between an overabundance of available RK due to its overexpression and that due to the lack of Rec regulation of RK in Rec knock-out animals. This discrepancy merits further investigation.
Our simulated results also suggest that Arr oligomerization has a significant effect on the timing of R* shutdown and on phototransduction recovery dynamics in general, particularly at brighter stimulus intensities. While it ultimately serves to maintain sufficient concentrations of monomers available to shut down R*, the transition of Arr from its oligomeric storage forms to the monomeric form slows down availability for and ultimately its activity in the shutdown of the activated receptor. Thus, RK overexpression has no effect on saturation time because the rate of R* shutdown is ultimately determined by Arr availability. When Arr oligomerization is disabled, on the other hand, high concentrations of Arr are immediately available to shut down R* as soon as it is sufficiently phosphorylated by RK; when RK is overexpressed, this occurs at a faster rate and thus Tsat is reduced.
Because Arr’s three oligomeric states are expected to exist in equilibrium concentrations in the dark (ref. ), some monomeric Arr is always available to quench R* after a weak stimulus. In our model under dark conditions, approximately 13% of Arr is present in monomeric form (~3.8e7 monomers = ~63 μM), 27% is present in dimeric form (~4e7 dimers = ~67 μM) and 60% is present in tetrameric form (~4.5e7 tetramers = ~75 μM), given 3e8 total Arr molecules (Table 2). It is known that the dark steady-state concentration of monomeric Arr varies widely between species, implicating Arr’s affinity for itself as a primary evolutionary target for the tuning of the regulation of R* shutdown . For example, while the formation of tetramers is highly cooperative in bovine photoreceptors, in mice the dissociation constants of dimer and tetramer formation are nearly the same. The resulting estimated concentrations of monomeric Arr in mice, cows and humans are 53 μM, 28 μM and 15 μM, respectively . Thus, the dark concentrations of the three Arr forms that we find in our model have realistic values. However, because the dissociation constant of Arr self-association has not been measured in any amphibian species, the true values may vary from those estimated here.
Steady-state protein concentrations
RecRCa · RK
Steady-state concentrations of Arr, RK and Rec in their various modeled states when no stimulus is presented. The model tracks molecular counts of each of these species. Equivalent micromolar concentrations were calculated using a rod photoreceptor outer segment volume of 1pL.
Finally, it is known that Arr undergoes translocation from the photoreceptor inner segment to the outer segment in a light-dependent manner . This occurs in the time-span of minutes, much longer than the duration of experiments described herein, thus we should not expect to see its effect in such experiments [33, 34]. Furthermore, many electrophysiological experiments, such as the experiments discussed here, use exclusively the rod outer segments, precluding any translocation. Nevertheless, the inclusion of Arr translocation in the model may result in a large portion of it being sequestered and completely unavailable during experiments on dark-adapted cells. As a result, the true kinetic parameters of Arr oligomerization may in fact be different than those that were herein estimated without translocation, depending on the time-course of the true concentrations of all Arr forms present in the outer segment during a photoresponse.
We have demonstrated that the present model, particularly after the inclusion of Arr oligomerization, can accurately simulate many complex dynamics of visual phototransduction, even under mutant conditions, including simultaneous overexpression and down-regulation of genes involved in the regulation of the recovery kinetics. The advantages of such a model are clear: by breaking the process down into fundamental steps, largely modeled using mass-action kinetics, it is simple to realistically simulate the effects of perturbations to the system. Because the model is composed in a modular manner, adding new mechanisms or replacing existing mechanisms with more realistic representations based on new experimental data can be done without greatly disrupting the rest of the model. Because the model is largely comprehensive, it is possible to test the long-reaching effects of the complex interactions and feedback loops present in this system. It is hoped that this model can serve as both a useful tool in future research by providing mechanistic insights as well as a unified compendium of knowledge of the visual phototransduction process.
While the phototransduction dynamics simulated by the present model are largely accurate, there remain some aspects that deserve further attention. Most of the primary mechanisms known to occur during the phototransduction process have been included in the model. However some, such as light-induced translocation of Arr as described above, as well as that of Rec (ref. ) and G (ref. [34, 36–39]), have yet to be integrated. Other processes, such as the action of phosducin or calmodulin, or the dynamic influence of the phosphorylation and dephosphorylation of phototransduction proteins other than R*, remain omitted from the model due to lack of mechanistic information. There remain processes which are present in the model in forms which assume dynamics other than those based on the law of mass action. For example, the mechanism of Ca2+-mediated activation of guanylate cyclase by GCAPs is not explicitly modeled. Improving its representation in the model would allow more flexible experimental manipulation of its parameters when probing its simulated dynamics.
Finally, we took advantage of the inherent modularity of the system to formulate and test the hypothesis that RK-mediated R* shutdown, particularly when RK is overexpressed, is modulated by Arr availability in its monomeric form. Integrating two new reactions into the model required relatively little compared to the time and costs involved in generating and maintaining transgenic or cross-bred animals. With the knowledge provided by the model, focused laboratory experiments to formally test this hypothesis can be undertaken with more confidence than had they been performed without the data provided herein. Indeed, this is one of the largest over-arching goals of the field of systems biology.
A previously published model of the phototransduction system was used as a template on which our changes were introduced , after first incorporating a R*-G pre-coupling mechanism described using the same model . The network of reactions comprising the system was modeled in a deterministic manner, with each reaction being represented by an ordinary differential equation (for the full network of reactions see ref. ). Except where otherwise noted, the law of mass action was used to describe the kinetics. The change with time in the number of molecules of each chemical species was simulated by summing all of the reaction rates which produce the species and subtracting all of the rates which remove the species from the system. The model simulates the reactions taking place in a well-stirred volume, hence the spatial structure of the photoreceptor outer segment is not taken into account. The reactions building up the network were implemented as described in Additional file 1: Table S1 of Dell’Orco et al  and in Dell’Orco and Koch , and the parameters that needed tuning for the present implementation are reported in Table 3.
Changes in parameter values relative to the model of Dell’Orco and Koch . Reaction numbers prefixed with “DO-” refer to reaction numbers listed in the Additional file 1 of ref. ; otherwise they refer to reactions in the present work.
Before the described modifications were made to the model, the relationship between the level of R* phosphorylation and its affinity for other proteins was reconsidered. Previously, the affinity of R* for Gt and RK was modeled to decrease exponentially with each additional phosphorylation while the affinity of R* for Arr increased exponentially. Experimental evidence suggests that the affinity of R* for Gt does, indeed, decrease exponentially with additional phosphorylations . The affinity of R* for Arr, on the other hand, has been shown to increase linearly with the first four phosphorylations [40, 41]. In the model this relationship was changed from exponential to saturating linear. In accordance with experimental evidence, the affinity of Arr for R* was modeled to increase linearly for one-to-four phosphorylations and to remain constant for five or six. Thus, the on-rate was determined as
where n is the number of phosphorylations, k Arr is the rate of binding of Arr to singly-phosphorylated R* and mArr is the rate of increase in affinity with each of the first four additional phosphates.
The effect of phosphorylation on the relationship between R* and RK is less well-understood. It has been suggested that the affinity of RK for R* depends both on its auto-phosphorylation as well as on that of R*, such that its affinity decreases with each phosphate added to either protein . However, only the extreme cases of no phosphorylation and full phosphorylation were measured. We found that the model’s dynamics were strongly influenced by the previous exponential R*-RK relationship, particularly with regards to the RK overexpression experiments. We thus replaced it with a linear decrease in RK affinity for R* with increasing R* phosphorylation, which fit the model equally well. The rate of binding of RK to R* was defined as follows:
where mRK was set to kRK10/5. Interestingly, it was found that setting the denominator of the slope to 5 led to the best fit, meaning that RK's affinity for R* is extinguished before R* reaches its maximum phosphorylation state.
The newly built model was tested alongside the previously published model, in a variety of light stimulation paradigms and under a variety of mutational scenarios. For example, the model reproduced the typical behavior of light adaptation in the presence of a non-saturating background light, consisting of reduced sensitivity and accelerated recovery (Additional file 1: Figure S1). When the signaling effector regulating protein RGS was knocked out, signal recovery was found to be drastically slower than normal, with perhaps more accurate kinetics than what was previously simulated (Additional file 1: Figure S2).
Previously, Ca2+-mediated feedback on RK by Rec was modeled based on a quasi-steady-state assumption using a standard Hill equation to determine the quantity of Rec in the “relaxed”, Ca2+-bound form [14, 15]. In order to more-accurately model the Ca2+-mediated feedback on RK by Rec, we replaced the original kinetics based on the quasi-steady-state assumption (Reaction 30 in ref. ) with a set of two fundamental reactions. First, Rec undergoes a conformational change upon binding Ca2+, from a “tense” form to a “relaxed” form with a covalently bound myristoyl group exposed (Figure 1):
To accommodate the difference in measurement in the model of Rec (number of molecules) and Ca2+ (micromolar concentration), when determining the rate of change of Ca2+ concentration the molecular quantities of RecT and RecRCa are dynamically converted to molar concentration using a volume of 1pL.
Dell’Orco et al. measured k Rec1 to have a value of approximately 0.011 μM-1 s-1 and kRec2 was measured to be 0.05 s-1 using surface plasmon resonance , and these values were found to fit the data well. In the previous version of the model, the conformational change did not explicitly affect the Ca2+ concentration. Instead, the ion’s binding to Rec was included implicitly in its binding to a concentration of anonymous Ca2+ buffers (parameter eT). To account for Rec’s activity as an explicit buffer in the present model, the concentration eT was reduced by an appropriate amount.
Following the conformational change to the relaxed state, Rec is free to bind with RK:
kRec3 and kRec4 are analogous to the synonymous parameters in previous versions of the model. The values of kRec3 and kRec4 were held at the previously used values of 9.69 μM-1 s-1 and 0.61 s-1, respectively (Table 3), which had been attained by parameter estimation techniques. To accommodate the representation in the present model of Rec state in number of molecules rather than the previous representation in micro-molar concentration, the value of kRec3 was converted into units of s-1, using an outer segment volume of 1pL. Because RK and Rec reach a steady state in the dark, the initial quantities of RecT, RecRCa, RecRCa · RK, and RK were set at their steady-state values for all simulated experimental conditions (Table 2). We stress the fact that while minimal parameter tuning is necessary when new mechanisms are introduced in the model, the values of the optimized parameters were not subsequently changed in further simulations. Hence, all the predicted dynamic behaviors have to be considered as robust model validations.
The oligomerization of Arr was modeled according to mass-action kinetics in a simple two-reaction process. First, two Arr monomers bind reversibly to form an Arr dimer. Next, two dimers reversibly bind to form a tetramer. Thus, the first reaction is
while the second reaction is
To minimize the number of introduced parameters, both reactions were assumed to occur at the same rates. Because the rates of these reactions have not been measured in amphibians, we chose to approximate them based on the known dissociation constants in mice, which approximately equal 60 μM for both reactions [32, 44]. To this end, kA4 was estimated by parameter optimization techniques while kA5 was held fixed to maintain a KD of 60 μM given an outer segment volume of 1pL (Table 3).
Since the concentration of Ca2+ changes dynamically upon illumination and the cation triggers regulation mechanisms throughout the cascade, we considered the possibility of a heretofore unknown mechanism of direct Ca2+ feedback on Arr oligomerization. This was implemented by adding an exponential dependence on free Ca2+ concentration as a fraction of the dark Ca2+ concentration for the two forward rates described above. For example, the forward rate of the first reaction would be
This extra mechanism did not impart any significant improvement on the fit of the model to experimental data. It was concluded that the indirect feedback of Ca2+ on Arr via Rec inhibition of RK was sufficient.
Model implementation, parameter estimation and numerical simulations
Because of the novel mechanisms in the model, some signal recovery-related parameters from the original model required re-tuning. Parameter optimization was performed using local and global methods to minimize the change in quantitative response characteristics produced by the previous version of the model to 0.024 s flashes (1 R*/flash, 985 R*/flash, 15,600 R*/flash and 118,000 R*/flash) and 60s steps of light (48 R*s-1, 220 R*s-1, and 450 R*s-1). Local optimization was performed using the Nelder-Mead Simplex Method while global optimization was performed using a particle swarm method . Parameter values were held fixed in all subsequent simulation experiments. Overall, the new model required moderate adjustments (Table 3). The Arr-related parameters required more significant re-tuning, which is not surprising given the change from an exponential dependence on R* phosphorylation state to a linear one and the addition of Arr oligomerization. The change in the relationship of RK affinity for R* and the phosphorylation state of R* from an exponential one to a linear one required significant re-tuning of the basal rate of binding of Gt to R* (kG10). Because the affinity of RK for R* now decays at a linear rate, it continues to bind and phosphorylate R* when R* has more phosphates attached compared to the previous implementation; as a result, R* accumulates phosphates more quickly, resulting in a concomitant rapid decrease in the rate of Gt binding. To overcome this, the basal rate of Gt binding was required to be significantly faster. Note, however, that to the best of our knowledge, this rate has not been measured for amphibian Gt, thus the actual value is unknown. However, a faster interaction between Gt and R* is in line with recent surface-plasmon resonance determinations .
The model was implemented using SBTOOLBOX2 for Matlab (http://www.sbtoolbox2.org) . SBTOOLBOX2 or SBML model files are available upon request. All numerical simulations were carried out in this framework, including parameter estimation. Deterministic simulations were run from automatically generated and compiled C-code models, based on the CVODE integrator from SUNDIALS .
Model validation and experimental simulations
Simulated stimulus intensities were adapted from the published experimental procedures. Where experimental stimuli were measured in photons μm-2, they were converted to R* per flash using a collecting area of 0.4 μm2. Where experiments were performed on mice, stimulus intensities were doubled, to qualitatively account for species differences in sensitivity between the model and the experimental animals.
To facilitate quantitative comparisons of Tsat and τD measurements between simulated experiments, we opted to standardize the methodology rather than simulate exactly the methodology of the published experiments. More specifically, whereas experimental procedures may define saturation time as the time spent above either 80% or 90% of the maximum current suppression, and τD may be calculated as the slope of either the first three or four points of the linear phase in a Pepperberg plot , we opted to exclusively measure currents above 90% of the maximum as saturating and to estimate τD from the first four data points on the curve. Furthermore, stimulus intensities for Pepperberg plots were standardized to half-log steps, with the same set of stimuli used for all such plots, regardless of those originally used in the published experiments.
Guanylate Cyclase Activating Proteins
Ca2+-bound, “relaxed” recoverin
Regulator of G-protein Signaling (protein)
Photoisomerized (activated) rhodopsin
The length of time that a saturating photoresponse remains at greater than 90% of its maximum current suppression
The time constant of a single exponential function fit to the second half of the recovery phase of a non-saturating photoresponse
The dominant time constant of recovery from a saturating response, measured as the slope of Tsat over logarithmically increasing stimulus intensities
This research was funded by grant BFU2010-19443 (subprogram BMC) awarded by Ministerio de Ciencia y Tecnología (Spain) and by the Direccío General de Recerca, Generalitat de Catalunya (Grup de Recerca Consolidat 2009SGR 1101). BMI is supported by FI-DGR and BE-DGR grants from AGAUR, Generalitat de Catalunya (2011 F1 B1 00275). LM acknowledges funding from the Juan de la Cierva Program of the Spanish Ministry of Science and Innovation (MICINN). DDO acknowledges funding from the Hanse-Wissenschaftkolleg Delmenhorst and from the Italian Ministry for Research and Education (Fur2011).
IBE – Institute of Evolutionary Biology (CSIC-Universitat Pompeu Fabra), CEXS-UPF-PRBB
Department of Neurosciences, Biochemistry Group, University of Oldenburg
Department of Life Sciences and Reproduction, Section of Biological Chemistry and Center for BioMedical Computing (CBMC), University of Verona
Pugh EN Jr, Lamb TD: Phototransduction in vertebrate rods and cones: molecular mechanisms of amplification, recovery and light adaptation. In Handbook of Biological Physics (Volume 3). Edited by: Stavenga DG, DeGrip WG, Pugh EN. North-Holland: Elsevier; 2000:183–255. doi:10.1016/S1383-8121(00)80008-1
Ames JB, Ishima R, Tanaka T, Gordon JI, Stryer L, Ikura M: Molecular mechanics of calcium-myristoyl switches.Nature 1997, 389:198–202.PubMedView Article
Lange C, Koch K-W: Calcium-dependent binding of recoverin to membranes monitored by surface plasmon resonance spectroscopy in real time.Biochemistry 1997, 2960:12019–12026.View Article
Kawamura S: Rhodopsin phosphorylation as a mechanism of cyclic GMP phosphodiesterase regulation by S-modulin.Nature 1993, 362:855–857.PubMedView Article
Gorodovikova EN, Senin II, Philippov PP: Calcium-sensitive control of rhodopsin phosphorylation in the reconstituted system consisting of photoreceptor membranes, rhodopsin kinase and recoverin.FEBS Lett 1994, 353:171–172.PubMedView Article
Chen C-K, Inglese J, Lefkowitz RJ, Hurley JB: Ca-dependent interaction of recoverin with rhodopsin kinase.J Biol Chem 1995, 270:18060–18066.PubMedView Article
Klenchin VA, Calvert PD, Bownds MD: Inhibition of rhodopsin kinase by recoverin.J Biol Chem 1995, 270:16147–16152.PubMedView Article
Higgins MK, Oprian DD, Schertler GFX: Recoverin binds exclusively to an amphipathic peptide at the N terminus of rhodopsin kinase, inhibiting rhodopsin phosphorylation without affecting catalytic activity of the kinase.J Biol Chem 2006, 281:19426–19432.PubMedView Article
Komolov KE, Senin II, Kovaleva NA, Christoph MP, Churumova VA, Grigoriev II, Akhtar M, Philippov PP, Koch K-W: Mechanism of rhodopsin kinase regulation by recoverin.J Neurochem 2009, 110:72–79.PubMedView Article
Hirsch JA, Schubert C, Gurevich VV, Sigler PB: The 2.8 A crystal structure of visual arrestin: a model for arrestin’s regulation.Cell 1999, 97:257–269.PubMedView Article
Schubert C, Hirsch JA, Gurevich VV, Engelman DM, Sigler PB, Fleming KG: Visual arrestin activity may be regulated by self-association.J Biol Chem 1999, 274:21186–21190.PubMedView Article
Hanson SM, Van Eps N, Francis DJ, Altenbach C, Vishnivetskiy SA, Arshavsky VY, Klug CS, Hubbell WL, Gurevich VV: Structure and function of the visual arrestin oligomer.EMBO J 2007, 26:1726–1736.PubMedView Article
Hamer RD, Nicholas SC, Tranchina D, Liebman PA, Lamb TD: Multiple steps of phosphorylation of activated rhodopsin can account for the reproducibility of vertebrate rod single-photon responses.J Gen Physiol 2003, 122:419–444.PubMedView Article
Hamer RD, Nicholas SC, Tranchina D, Lamb TD, Jarvinen JLP: Toward a unified model of vertebrate rod phototransduction.Vis Neurosci 2005, 22:417–436.PubMedView Article
Dell’Orco D, Schmidt H, Mariani S, Fanelli F: Network-level analysis of light adaptation in rod cells under normal and altered conditions.Mol Biosyst 2009, 5:1232–1246.PubMedView Article
Kiel C, Vogt A, Campagna A, Chatr-aryamontri A, Swiatek-de Lange M, Beer M, Bolz S, Mack AF, Kinkl N, Cesareni G, Serrano L, Ueffing M: Structural and functional protein network analyses predict novel signaling functions for rhodopsin.Mol Syst Biol 2011, 7:551.PubMedView Article
Caruso G, Bisegna P, Andreucci D, Lenoci L, Gurevich VV, Hamm HE, Dibenedetto E: Identification of key factors that reduce the variability of the single photon response.Proc Natl Acad Sci U S A 2011, 108:7804–7807.PubMedView Article
Holcman D, Korenbrot JI: Longitudinal diffusion in retinal rod and cone outer segment cytoplasm: the consequence of cell structure.Biophys J 2004, 86:2566–2582.PubMedView Article
Dell’Orco D, Koch K-W: A dynamic scaffolding mechanism for rhodopsin and transducin interaction in vertebrate vision.Biochem J 2011, 440:263–271.PubMedView Article
Dell’Orco D, Koch K-W: Systems biochemistry approaches to vertebrate phototransduction: towards a molecular understanding of disease.Biochem Soc Trans 2010, 38:1275–1280.PubMedView Article
Moriondo A, Rispoli G: A step-by-step model of phototransduction cascade shows that Ca2+regulation of guanylate cyclase accounts only for short-term changes of photoresponse.Photochem Photobiol Sci 2003, 2:1292–1298.PubMedView Article
Forti S, Menini A, Rispoli G, Torre V: Kinetics of phototransduction in retinal rods of the newt Triturus cristatus.J Physiol 1989, 419:265–295.PubMed
Krispel CM, Chen D, Melling N, Chen Y-J, Martemyanov KA, Quillinan N, Arshavsky VY, Wensel TG, Chen C-K, Burns ME: RGS expression rate-limits recovery of rod photoresponses.Neuron 2006, 51:409–416.PubMedView Article
Sakurai K, Young JE, Kefalov VJ, Khani SC: Variation in rhodopsin kinase expression alters the dim flash response shut off and the light adaptation in rod photoreceptors.Invest Ophthalmol Vis Sci 2011, 52:6793–6800.PubMedView Article
Pepperberg D, Cornwall M, Kahlert M, Hofmann K, Jin J, Jones G, Ripps H: Light-dependent delay in the falling phase of the retinal rod photoresponse.Vis Neurosci 1992, 8:9–18.PubMedView Article
Gross OP, Burns ME: Control of rhodopsin’s active lifetime by arrestin-1 expression in mammalian rods.J Neurosci 2010, 30:3450–3457.PubMedView Article
Makino CL, Dodd RL, Chen J, Burns ME, Roca A, Simon MI, Baylor DA: Recoverin regulates light-dependent phosphodiesterase activity in retinal rods.J Gen Physiol 2004, 123:729–741.PubMedView Article
Chen C-K, Woodruff ML, Chen FS, Chen D, Fain GL: Background light produces a recoverin-dependent modulation of activated-rhodopsin lifetime in mouse rods.J Neurosci 2010, 30:1213–1220.PubMedView Article
Gray-Keller MP, Polans AS, Palczewski K, Detwiler PB: The effect of recoverin-like calcium-binding proteins on the photoresponse of retinal rods.Neuron 1993, 10:523–531.PubMedView Article
Burns ME, Pugh EN: RGS9 concentration matters in rod phototransduction.Biophys J 2009, 97:1538–1547.PubMedView Article
Gurevich VV, Hanson SM, Song X, Vishnivetskiy SA, Gurevich EV: The functional cycle of visual arrestins in photoreceptor cells.Prog Retin Eye Res 2011, 30:405–430.PubMedView Article
Kim M, Hanson SM, Vishnivetskiy SA, Song X, Cleghorn WM, Hubbell WL, Gurevich VV: Robust self-association is a common feature of mammalian visual arrestin-1.Biochemistry 2011, 50:2235–2242.PubMedView Article
Strissel KJ, Sokolov M, Trieu LH, Arshavsky VY: Arrestin translocation is induced at a critical threshold of visual signaling and is superstoichiometric to bleached rhodopsin.J Neurosci 2006, 26:1146–1153.PubMedView Article
Elias RV, Sezate SS, Cao W, McGinnis JF: Temporal kinetics of the light/dark translocation and compartmentation of arrestin and alpha-transducin in mouse photoreceptor cells.Mol Vis 2004, 10:672–681.PubMed
Slepak VZ, Hurley JB: Mechanism of light-induced translocation of arrestin and transducin in photoreceptors: Interaction-restricted diffusion.IUBMB Life 2007, 60:2–9.View Article
Lobanova ES, Finkelstein S, Song H, Tsang SH, Chen C-K, Sokolov M, Skiba NP, Arshavsky VY: Transducin translocation in rods is triggered by saturation of the GTPase-activating complex.J Neurosci 2007, 27:1151–1160.PubMedView Article
Kerov V, Chen D, Moussaif M, Chen Y-J, Chen C-K, Artemyev NO: Transducin activation state controls its light-dependent translocation in rod photoreceptors.J Biol Chem 2005, 280:41069–41076.PubMedView Article
Sokolov M, Lyubarsky AL, Strissel KJ, Savchenko AB, Govardovskii VI, Pugh EN, Arshavsky VY: Massive light-driven translocation of transducin between the two major compartments of rod cells: a novel mechanism of light adaptation.Neuron 2002, 34:95–106.PubMedView Article
Gibson SK, Parkes JH, Liebman PA: Phosphorylation modulates the affinity of light-activated rhodopsin for G protein and arrestin.Biochemistry 2000, 39:5738–5749.PubMedView Article
Vishnivetskiy SA, Raman D, Wei J, Kennedy MJ, Hurley JB, Gurevich VV: Regulation of arrestin binding by rhodopsin phosphorylation level.J Biol Chem 2007, 282:32075–32083.PubMedView Article
Buczyłko J, Gutmann C, Palczewski K: Regulation of rhodopsin kinase by autophosphorylation.Proc Natl Acad Sci 1991, 88:2568–2572.PubMedView Article
Dell’Orco D, Müller M, Koch K-W: Quantitative detection of conformational transitions in a calcium sensor protein by surface plasmon resonance.Chem Commun 2010, 46:7316–7318.View Article
Song X, Vishnivetskiy S, Seo J, Chen J, Gurevich E, Gurevich V: Arrestin-1 expression level in rods: balancing functional performance and photoreceptor health.Neuroscience 2011, 174:37–49.PubMedView Article
Vaz AIF, Vicente LN: A particle swarm pattern search method for bound constrained global optimization.J Glob Optim 2007, 39:197–219.View Article
Schmidt H, Jirstrand M: Systems Biology Toolbox for MATLAB: a computational platform for research in systems biology.Bioinformatics 2006, 22:514–515.PubMedView Article
Hindmarsh A, Brown P, Grant K: SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers.ACM Trans Math Softw 2005, 31:363–396.View Article
Matthews G: Spread of the light response along the rod outer segment: An estimate from patch-clamp recordings.Vision Res 1986, 26:535–541.PubMedView Article
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.