article: The single photon response (SPR) in vertebrate phototransduction is regulated by the dynamics of R* during its lifetime, including the random number of phosphorylations, the catalytic activity and the random sojourn time at each phosphorylation level. Because of this randomness the electrical responses are expected to be inherently variable. However the SPR is highly reproducible. The mechanisms that confer to the SPR such a low variability are not completely understood. The kinetics of rhodopsin deactivation is investigated by a Continuous Time Markov Chain (CTMC) based on the biochemistry of rhodopsin activation and deactivation, interfaced with a spatio-temporal model of phototransduction. The model parameters are extracted from the photoresponse data of both wild type and mutant mice, having variable numbers of phosphorylation sites and, with the same set of parameters, the model reproduces both WT and mutant responses. The sources of variability are dissected into its components, by asking whether a random number of turnoff steps, a random sojourn time between steps, or both, give rise to the known variability. The model shows that only the randomness of the sojourn times in each of the phosphorylated states contributes to the Coefficient of Variation (CV) of the response, whereas the randomness of the number of R* turnoff steps has a negligible effect. These results counter the view that the larger the number of decay steps of R*, the more stable the photoresponse is. Our results indicate that R* shutoff is responsible for the variability of the photoresponse, while the diffusion of the second messengers acts as a variability suppressor.
  In retinal rod photoreceptors, rhodopsin activated by photons of light, denoted by, initiates a signal transduction cascade to produce a suppression of electrical current flowing into rod outer segment (ROS). Following isomerization, a molecule undergoes a random number of phosphorylations by rhodopsin kinase (RK) and finally is inactivated by arrestin (Arr) binding. Activated rhodopsin, moving along its random path, during its random lifetime from isomerization to Arr binding, keeps activating its cognate G-protein (G) transducin, while its catalytic activity declines with increasing level of phosphorylation. The active G-protein () associates with the effector protein phosphodiesterase (E) forming an active - complex, which by hydrolyzing cGMP reduces its concentration, thereby generating a current response on the outer shell of the ROS. The dynamics of during its lifetime, including the random number of phosphorylations, the catalytic activity and the random sojourn time at each phosphorylation level, regulates the production of and therefore the current response. Because of the randomness in the components of the activation/deactivation cascade, the electrical responses are expected to be inherently variable. However, the single photon response (SPR) exhibits a low variability in the sense that the amplitude and shape of the electrical responses, corresponding to a set of activation-deactivation events, are similar. It is reported that the Coefficient of Variation (CV = standard deviation/mean) of the SPR area for mouse is about [1]. However, the mechanisms that confer high reproducibility of the SPR are not completely understood. Several studies [1]–[7] attribute the high reproducibility of the SPR mainly to the mechanisms regulating rhodopsin deactivation. Although the models proposed in these studies account for the low variability of the response, they impose, in one way or another, certain restrictions on the biochemistry of rhodopsin deactivation. For example, if rhodopsin' s integrated activity occurs in independent steps, it is assumed that each step controls an equal fraction of rhodopsin' s integrated catalytic activity [1], [2]. It is then natural to ask what is the statistical mean of the number, as a way of testing both the models and the supporting biochemical assumptions. Mechanistically, one might ask which of the components of the deactivation cascade contribute more importantly to the variability. A major difficulty with these issues is to experimentally separate the various components that contribute to the variability. To our knowledge, the activation/deactivation module of the cascade is not, to date, experimentally separable from the transduction module. We have shown in [8] that diffusion of the second messengers in the cytoplasm acts as a variability suppressor. The separation between the activation cascade on the disks and the diffusion of the second messengers cGMP and in the cytoplasm is realized by a mathematical model [8]–[11]. Likewise several fine properties of the biochemical and biophysical mechanisms regulating the recovery and reproducibility of SPR are not, to our knowledge, experimentally separable. Here we attempted to tease apart the various components of the shutoff mechanism and analyze to what extent each of them contributes to the variability of the SPR. Unlike the transduction part of the cascade, where the intricacy is of geometrical nature [9]–[11], the main difficulty here is stochastic. Rhodopsin inactivation can occur by several mechanisms, including Arr binding and thermal decay to opsin. We only model the former, as the latter occurs on a much longer time scale [12], [13]. Shutoff of by Arr binding can follow, in principle, an infinite number of paths, depending on the random number of phosphorylated states, and the random sojourn times in those states. The biochemistry that regulates rhodopsin deactivation is put into a stochastic framework, which reproduces the SPR both in WT and in mutant mice, and is capable of analyzing the randomness of each phosphorylation state of. This is interfaced with the spatio-temporal model in [8], [9], [11], capable of tracking the diffusion of the second messengers in the cytoplasm and of detecting the effects of geometrical changes of the ROS on the photoresponse. We find that the randomness of the sojourn times of in each of its phosphorylation states acts as the dominant factor contributing to the CV of the response. At the same time the number of available phosphorylation sites or the random number of phosphorylations before shutoff, is shown to contribute little to variability suppression. We also find that, in addition to changed biochemistry, the geometry of the ROS might be important for the light response in mutant mice.
  Deactivation of rhodopsin with one or several mutant phosphorylation sites, can be simulated by suitable choices of the sequences and as indicated in the section of numerical procedures and methods. Mutant mouse rhodopsins bearing fewer than 6 phosphorylation sites generate SPRs of significantly extended durations (Figure 2A). The rate of recovery increases with increasing numbers of phosphorylation sites (Figure 2A), in qualitative and quantitative agreement with the experimental results of [3] (Figure 2B). Inactivation of all rhodopsin phosphorylation sites is realized by either mutation of all six serines and threonines to alanines [1], [3], or rhodopsin kinase knockout [23]. The corresponding SPRs are similar, exhibiting larger amplitude and longer duration than WT (Figure 2B, D for (0P) ). A prolonged SPR in mutant mouse rods lacking arrestin is reported in [24] (Figure 2D). This is realized by setting for all in the model. The activated rhodopsin gets phosphorylated until all six sites are occupied. Its activity is reduced with increased phosphorylations, and kept fixed after the last phosphorylation for the remainder of the process. The remaining activity yields a response with an asymptotic tail at almost half of its peak value. The initial fall of the response is triggered by phosphorylation. The simulations are shown by (–/–) in Figure 2C, and are qualitatively and quantitatively in agreement with the experimental studies of [24] (Figure 2D). In Table 2 we report the simulated characteristics of SPRs from WT rods and those expressing rhodopsin mutants. By increasing the number of phosphorylation sites, the peaks of the current response decrease; the time to peak decreases; and the SPR area decreases significantly. For mutants that exhibit very slow recovery (0P, 1P, 2P) the corresponding is large because the current remains high for an extended period of time. The value of has been computed by integrating the photocurrent over the time of simulation (). When only one phosphorylation site was mutated, the SPR was almost like that of WT but recovery was slightly slower. Consistent with this slower recovery, the SPR area of the response of rhodopsin with five phosphorylation sites (5P) was about larger than those for wild type. Taken together, these results are consistent with the experimental observations of [3] and the notion that normal kinetics of deactivation requires the presence of all six phosphorylation sites. We finally comment on the largest rising curves coded in red in Figure 2B and D. Various experimental studies [3], [23] show that the response amplitude for the case (0P) is roughly twice as large as the response for the case (1P). In [3] the case (0P), is realized by CSM, and in [23] by RK knockout. In both cases all phosphorylation sites are removed or made inoperative, and both cases exhibit the double amplitude response, suggesting a common mechanism. This issue is not discussed in the indicated papers and we are not aware of an explanation or hypothesis for a possible biochemical mechanism. However, Figure 2 of [3] shows that the ROS in CSM mice were about shorter than WT. Geometrical changes due to genetic manipulations are also discussed in [23] (page 3720), and [24], (Figure 3d, page 506). We repeated the simulations with a ROS whose height was reduced by, while all the remaining parameters were kept fixed. In particular, the number of channels was kept fixed, thereby increasing their density. Since the response is localized close to the activation site [10], [11], the augmented channel density yields a larger response. The resulting simulation is reported in Figure 2A for (0P*) as the largest amplitude (red curve). While the agreement with corresponding experimental curve in Figure 2B is striking, at this point we refrain from suggesting that this as the only functional mechanism. The CTMC model permits one to test independently the effects of the random components of the variability on the response. For example one can separate the effects of the randomness of the sojourn time from the randomness of number of shutoff steps. To achieve this, we performed the following sets of simulations: Stochastic simulations are effected for WT and each of the knock-out cases of COOH-terminal truncations [24], [25] and RK knockout [1], [3]. After about 5,000 numerical simulations, up to 3 s, mean, standard deviation and CV are computed for effector and normalized current suppression. Further technical details are in Methods.
  Variability of the photoresponse hinges on a coordinated system behavior of several components. The main two modules are the activation/deactivation part and the transduction part of the cascade. The latter, given its input, is essentially deterministic as it involves the diffusion of the second messengers cGMP and in the cytoplasm and a subsequent current drop through the closure of the cGMP -gated channels. The former is essentially stochastic as it involves the biochemistry of rhodopsin shutoff, which occurs in several random steps. An understanding of the process hinges upon teasing apart all these components, analyzing them separately and blending them together into a system behavior. This point of view began in [8], by separating the role of the transduction from that of the activation/deactivation. This separation was made possible by a mathematical model capable of distinguishing the biochemistry of shutoff, from the functional role of the transduction [8], [9], [11]. A surprising finding was that, while shutoff is responsible for the variability of the photoresponse, the diffusion of the second messengers acts as a variability suppressor. Here we have further separated the various steps of the deactivation cascade by (a) prescribing a probabilistic mechanism (CTMC) by which the system selects its random states, and (b) by interrogating the known biochemistry to trace patterns and parameters. It is not sufficient to determine these parameters unambiguously using WT mice. Experimental information from some mutant and knock-out animals is needed. Specifically, the choice of the catalytic activities by formula Eq: 6, while based on known biochemistry [26], hinges upon the basic parameter, which in turn is determined by the biochemistry of the cascade in mutant mouse (section on parameters in Methods). The same holds true for the transition parameters, given by formula Eq: 7 and depending upon the parameter. Thus, a first remark is that our approach, while mathematical and computational, parallels the biology; that is, information is extracted in a complementary way from the data on genetically modified as well as WT animals. Next the model populated by the indicated parameters is validated against WT and mutant responses as in Figure 2. The model has a deterministic component, and a stochastic component. The first regards the transduction part of the cascade, which is geometry dependent, and deterministic, being based on the diffusion of the second messengers cGMP and in the cytoplasm. Importantly, this model permits one to test the response against geometrical variations of the ROS. The response in mice expressing CSM or RK knock out is rather unusual, exhibiting a double amplitude with respect to WT [3], [23]. An examination of the immunofluorescence micrographs in Figure 2 of [3], suggests that the length of ROS in CSM mice is reduced by about relative to WT. Geometrical modifications presumably due to genetic manipulations are also discussed in [23]. Keeping the same stochastic biochemical scheme and changing the length of the ROS, the model reproduced the double-amplitude phenomenon described in [3], [23] (Figure 2 A, B), suggesting that the modified geometry of mutant ROS, might contribute, along with the changed biochemistry, to this phenomenon. This results, along with a recent study of rod signaling in mice expressing supra-physiological levels of rhodopsin ([27]), emphasize the importance of investigating the ROS geometry in genetically modified mouse lines. Our analysis shows that the changes in ROS length, which were analyzed in very few mouse lines, can have dramatic effects on photoresponse. The stochastic component permits one to single out those parts of the activation/deactivation cascade that dominantly contribute to the variability of the response. The main result is that variability is largely generated by the randomness of the sojourn times of in its phosphorylation states. The prevailing point of view is that the activation cascade is responsible for the variability, although in a non quantified way, and that deactivation of is responsible for variability suppression, and further, the larger the number of decay steps of, the more stable the photoresponse [1], [3]–[7]. This view was expressed in [1], where mice expressing rhodopsin with 0,1, 2,5, and 6 phosphorylation sites were used. The analysis presented in [1] has some inconsistencies. Although the experimental points seem to be best fitted by a straight line (Figure 1 of [1]) the authors describe them by, with being the number of available phosphorylation sites. The lines with 3 and 4 phosphorylation sites, which would have allowed to discriminate between these functions, were not analyzed in [1]. In addition, by comparing the CV of mice with 0,1, and 2 sites, which do not demonstrate rapid recovery ([4]), with those having 5 or 6 (WT) sites that recover with essentially the same fast rate, the authors inappropriately lump together two disparate phenomena. In the latter case, normal two-step rhodopsin inactivation by RK phosphorylation and arrestin binding is fully operative, whereas in the former rhodopsin is inactivated by stochastic thermal decay taking place on a much longer time course. The idea that multiple inactivation steps are necessary to suppress variability was recently expressed in [2], where the authors conclusions were largely based on two assumptions. The first is that activity is nearly equally distributed among the deactivation steps. The second is that in Ames' solution, that yields much greater and longer-lasting SPR than Locke' s ([2], [28]), rhodopsin inactivation is rate-limiting and dominates the recovery kinetics. The biochemical scheme we propose argues against the first assumption, on experimental grounds (see a discussion below and On the Parameters and). The second assumption has been recently questioned in [28], where the authors showed that RGS9 overexpression similarly accelerates the recovery measures in Locke' s and Ames' solutions, indicating that transducin inactivation is rate limiting in both cases. Additional issues with data analysis of [1] were discussed in [28]. Thus, no compelling experimental evidence that the number of inactivation steps reduces variability can be found in the literature. Our results offer a different perspective; demonstrating that variability is generated by the randomness of the sojourn times of in its phosphorylated states, and that increasing the number of these states does not lead to variability suppression. The number of steps to deactivation does not coincide with the number of available phosphorylation sites. The experimental studies of [20] suggest that one phosphorylation is not sufficient for Arr binding, and the probability of quenching becomes large after 3 phosphorylations. Specifically corresponds to by which Eq: 7–Eq: 8 give and hence by Eq: 9. Thus, the system remains indefinitely activated (in reality it is stochastically inactivated by the thermal decay of rhodopsin, which is too slow to be captured by 3 s simulations used here). The case for corresponds to respectively and one computes from Eq: 9 and hence from Eq: 12; the system goes through steps and then remains “indefinitely” active (see above about thermal decay). Thus the CV of and in Table 1 and Table S1 in the supplementary material, are not due to variations caused by shutoff by Arr binding. The first case when and deactivation is possible, is the case corresponding to reported in Eq: 13. Going from 3P to 6P, the and the remain essentially the same. To illustrate the rationale of our main results, consider mutant rhodopsin with only 3 available phosphorylation sites. Since 3 phosphates are needed for Arr binding [20], no randomness is present in the deactivation process, if randomness is only measured in terms of steps to shutoff. This suggest that the source of variability is in other components of the process. Table 1 indeed shows that if the sojourn times of in each of its phosphorylation states are taken to be deterministic (Case 2), then the variability of the photoresponse is negligible. If on the other hand such sojourn times are permitted to be random, then the variability rises to experimentally observed levels (Figure 1), both for WT and genetically modified with 3-5P. This should not be interpreted, however, as though the reproducibility decreases as the number of available phosphorylation sites increases. We stress that increasing does not necessarily mean that the mean number to shutoff increases. The latter depends on the biochemistry of the process via Eq: 10–Eq: 11. Likewise the expected average of the random lifetime of is generated by the biochemistry in Eq: 7–Eq: 12 and; in particular it is different for different genetically modified mice (0P, 1P, etc.). The lifetime of is randomly chosen by the biochemistry in each of its random trials. For WT mouse, and only in this case, the expected lifetime of, as defined by formula Eq: 9–Eq: 12, coincides with the experimentally measured, effective average lifetime. In [29] it is reported as an upper limit, whereas several recent studies [28], [30], [31] suggest that might be as low as 40 ms (see Parameters). Therefore, we performed all simulations for both values, which yielded very similar CVs, both functionally and numerically (Figure 1 and Tables 2–3, and Figure S1 and Tables S1, S3, in the supplementary material). These similarities suggest that reproducibility is independent of the actual value of and depends only on the functional, sequence of the deactivation cascade, as predicted by our biochemical scheme (Eq: 9– Eq: 12). Further remarks on these two parameters and corresponding CVs are in On the Parameters and. We stress that the model includes only the deactivation mechanism due to Arr binding and does not include inactivation due to other causes such as thermal decay to opsin occurring over a time course of ([13]). In [1] the is computed over a time course of over 15 s, which is beyond the time course of inactivation. According to our scheme, based on direct biochemical measurements of arrestin binding to separated rhodopsin species with different numbers of attached phosphates ([20]), cases, and do not permit shutoff by Arr binding and remains active much longer than the 3 s of our simulations. Thus the CV due to deactivation reflects its thermal decay to opsin ([12]). In this case, shutoff is an abrupt 1-step process, implying, by Poisson statistic, CV = 1. This is essentially what is reported in [1]. For the cases and, although the experiments of [1] are carried over a time course of 15 s the is essentially due to shutoff by arrestin binding, which occurs within a time course of 0. 1 s, whereas decay to opsin is much slower. Considering the slow rate of thermal inactivation of rhodopsin, the probability of thermal decay within the first 0. 1 s is negligible relative to the probability of decay due to Arr binding. Accordingly, the reported in [1] for and is similar, as we find. The crucial cases and were not measured in [1]. The CTMC scheme we propose here differs from the Poisson statistics used in [1], [2], where the CV of is claimed to be proportional to. It should be noted that the number of available sites does not coincide with the average number of steps to shutoff and that each step weighs differently in the deactivation process, due to its biochemical history. We return briefly to the explicit, theoretical formula Eq: 3, valid under the assumptions of Eq: 4, and hence for the cases 3-6P. We have already remarked that its theoretical values (for Case 1) are in agreement with our simulations (lines 2 and 3 of Table 1). If one would artificially concoct a biochemistry by which all the products are the same for all, then formula Eq: 3 would give. This occurrence might suggest that the CV of the photoresponse decreases as the reciprocal of the square root of the number of steps to shutoff. A calculation from Eq: 7–Eq: 11, in agreement with known biochemistry ([20]), shows that the products are not constant (Table 3). In addition, even if this were the case, the variability of the photocurrent is very different from that of, as the relation between these functionals is highly non-linear [8], [11]. A further examination of Table 3 for 3-6P, shows that in all cases (WT or mutant), only the first few steps contribute significantly to the total activity; the remaining ones being negligible. In view of the theoretical formula Eq: 3, this is further evidence that increasing the number of steps does not significantly decrease the CV (). In all cases (WT or mutant) we found that the diffusion of the second messengers cGMP and in the cytoplasm acts as the dominant variability suppressor, thereby confirming the results of [8] and extending the analysis to a variety of transgenic models. These results are made possible by separating the activation/deactivation module from the transduction module. In addition, in the activation/deactivation module, one further separates the biochemical effects of each phosphorylation contributing to the responses, thereby allowing an examination of the role of the underlying biochemistry during deactivation. Incorporating the sequence of biochemical steps, described in Methods, allowed us to recapitulate experimental results qualitatively and quantitatively (Figure 2). It is worth noting that with realistic biochemistry, where Arr acquires a high binding affinity after 3 phosphorylation steps [20], the number of inactivation steps actually involved in shutting down individual SPRs varies very little. Therefore the fact that this number contributes virtually nothing to SPR variability, is one of the mechanisms maintaining the reproducibility of SPR.
  The state diagram of the CTMC describing deactivation by Rhodopsin Kinase (RK) phosphorylating the C-terminal serines and threonines in rhodopsin, is shown in Figure 3 with circles and arrows denoting states and transitions respectively. The states are labeled by the indices, and the transitions between connected states are labeled by transition rates and. The catalytic activity in its state is. The number of phosphorylation levels is determined by the number of phosphorylation sites of rhodopsin, which varies in different species. In mouse, rhodopsin has six phosphorylation sites [3]. State 1 is the non-phosphorylated level, representing newly activated rhodopsin with catalytic activity; the state represents fully deactivated rhodopsin with catalytic activity; states 2 to represent different phosphorylation levels, in which rhodopsin holds sites available for phosphorylation, with sites already phosphorylated, and has catalytic activity. The states 1 to are active states and the state is the inactive state. Specifically for WT mouse, there are seven () active states, including state 1 where is active and not phosphorylated. Transitions between active states are governed by the phosphorylation rates. For notation consistency, we let. Transitions between active states and the inactive state are governed by the arrestin binding rates. Arrestin binds with high affinity only to phosphorylated rhodopsin [20], [26], [32], [33], therefore, . A newly isomerized rhodopsin is in state 1. It undergoes a random number of phosphorylations before it transitions to the fully deactivated state. A rhodopsin with available phosphorylation sites could be phosphorylated at most times to state. Generally, in state, rhodopsin either interacts with rhodopsin kinase adding one more phosphate and transitions to the next phosphorylation level with a rate of, or it binds arrestin which quenches its catalytic activity, and transitions to the inactive state with a rate of. This process is a Bernoulli trial with the probability of a further phosphorylation given by and the probability of arrestin binding given by. This statistical scheme permits one to model rhodopsin deactivation also in transgenic animals with different number of rhodopsin phosphorylation sites. For example, if of the phosphorylation sites are mutated, we could setto reflect the effect of the mutation. It should be pointed out that, given mutated sites, the model removes any of the available sites with no discriminating criterion such as their ordering on the C-terminus. Although the phosphorylation of different sites in WT rhodopsin apparently proceeds in some order [34], the overall number of rhodopsin-attached phosphates, rather than their positions on the C-terminus, determines arrestin binding [3], [20], [35]. Accordingly, the model treats them as equal by attributing them the same biochemical function of holding a phosphate. Let denote the probability that a single is in the state. Then the mathematical description of the CTMC model shown in Figure 3 is [14], [15] (5) Note that the integer used to label the state of is one plus the corresponding level of phosphorylation, which is. For example, the phosphorylation level 0 corresponds to state 1, and sites are phosphorylated in state. The sojourn time of in state, is taken as an exponentially distributed random variable with mean. The sequences of the phosphorylation rates by RK, the activities of Arr, and the catalytic activities, depend on the underlying biochemistry, and vary with phosphorylation levels [26], [36]. The catalytic activity of in the state is the production rate of activated G protein by. While is active in each state, including the first unphosphorylated state (), its activity decreases with increasing phosphorylation levels. The catalytic activity of rhodopsin with different numbers of attached phosphates was experimentally measured by Wilden [26]. In this study differentially phosphorylated rhodopsin species were actually separated, so the conclusions were based on direct measurements and did not involve untested assumptions. Although similar conclusions were later reached by Gibson et al. [32], these authors did not separate differentially phosphorylated rhodopsin species, using preparation with different average phosphorylation levels instead. Their calculations are based on the assumption of Poisson distribution of rhodopsin species with different number of phosphates ([32]). This does not appear realistic, considering the distribution determined by rhodopsin fractionation ([26], [33]). Therefore, our model assumes that the binding affinity of phosphorylated for decreases exponentially with each added phosphate. Thus, based on data published by Wilden [26], we assume (6) where is the catalytic activity of in its initial, unphosphorylated state, and is positive. The value of in Eq: 6 has been extracted from the published data after an extensive consistency and sensitivity analysis ([37]). The parameters and are linked and subject to mutual limitations. It has been shown that in arrestin knockout mice, the initial kinetics of single photon response deactivation closely resembles that of WT, whereas the later phase of deactivation is abrogated (Xu et al. , 1997). Initial deactivation is attributable to rhodopsin phosphorylation, which is preserved in these animals. Then deactivation stops at about of the peak current suppression, and remains essentially steady thereafter. This level of current drop reflects the ability of fully phosphorylated mouse rhodopsin to activate transducin, corresponding to the catalytic activity when all 6 sites are phosphorylated. Thus from Eq: 6 one has. Mutual calibration of and is discussed in the section on parameters. Here we stress that they are determined from experimental data for both WT and mutant mice, and not chosen by fitting. While the explicit dependence of -RK binding affinity and the phosphorylation rates on the various biochemical states is not known, there is qualitative biochemical support for the notion that -RK affinity decreases systematically with the phosphorylation level of [38]. It is shown in Buczylko et al. [39], that phosphorylated RK has significantly lower ability to phosphorylate already phosphorylated than unphosphorylated. Moreover, Mendez et al. [3] showed that the rate of deactivation depends not on the identity of the available sites, but on their total number. We used the biochemically realistic assumption that the rate of phosphorylation is proportional to the number of serines and threonines still available for RK on the rhodopsin molecule. Mechanistically this means that the probability that upon binding to light-activated rhodopsin RK dissociates without adding another phosphate increases with the number of phosphates present, reaching 1 when all six sites are already phosphorylated. This assumption is consistent with in vivo observations by Mendez et al [3] that the removal of even one or two rhodopsin phosphorylation sites slows down photoresponse inactivation. Note that is the rate at which RK phosphorylates in its state. It depends on the on-rate of RK binding to in this state, and the rate of phosphate transfer, which were never separated experimentally and were not separated in our model. We set the sequence as linearly decreasing by increasing phosphorylation levels, that is the phosphorylation rate is proportional to the total number of the available sites and is independent of their biochemical identity. Thus (7) where is a rate constant, discussed and calibrated in the section on parameters. Formula Eq: 7 can be arrived at by postulating that RK has a fixed affinity for binding to and that each of the phosphorylation sites could be occupied with an equal rate. Therefore phosphorylation rate depends on the number of phosphorylation sites available for RK. Since in state has available phosphorylation sites it has a phosphorylation rate. Note that this model, similar to previously proposed ones, is based on the assumption that a single site is phosphorylated as a result of each rhodopsin encounter with RK. This assumption has not been experimentally tested. Arrestin binding ensures the timely termination of signaling, and it depends on the -Arr affinity. Several studies [26], [32], [36], suggest that arrestin affinity increases with increasing phosphorylation levels. Note that only the “irreversible” binding that terminates rhodopsin activity is taken into account here; that is the binding which occurs with high enough affinity to make the complex half-life much greater than the time course of the SPR. In a recent study, Vishnivetskiy et al. [20] found that unphosphorylated and mono-phosphorylated show the same low Arr binding levels. In particular, a single receptor-attached phosphate does not facilitate Arr binding; two are necessary to induce higher affinity interaction, and with three phosphates is fully capable of binding Arr with the affinity that makes the interaction essentially irreversible on the time scale of the SPR. Moreover, higher phosphorylation levels do not increase the stability of Arr complex with light-activated rhodopsin [20]. Based on the data in [20], we set the sequence for Arr binding rate by the phosphorylation level as (8) where is the Arr binding rate when Arr affinity reaches its maximum after several phosphorylations. Note that in the model describes arrestin binding that terminates transducin activation. Thus, it reflects the rate of formation of arrestin-rhodopsin complexes that are stable enough to survive significantly longer than the time course of a single photon response analyzed here. Since the stability of arrestin complex with unphosphorylated and mono-phosphorylated rhodopsin is much lower [20], [40]–[42], allowing for arrestin dissociation and consequent rhodopsin reactivation within this time, we set. Since, in states surely transitions to the states respectively. The data in vivo [3] and in vitro [20] also suggest that two rhodopsin-attached phosphates are not sufficient to induce Arr binding with high enough affinity for rapid deactivation. Therefore, we set. The effect of the level of rhodopsin phosphorylation on arrestin binding was explored in two studies. Gibson et al [32] concluded that arrestin affinity linearly increases with the level of phosphorylation in the range of 1–4 phosphates per rhodopsin. The authors used preparations of phosphorylated rhodopsin in native disc membranes that are well known to be highly heterogeneous, containing rhodopsin species carrying from zero to seven phosphates (bovine rhodopsin has seven RK phosphorylation sites) [26], [43]. The authors attempted to solve this problem by using several assumptions (that were not experimentally tested) to compute the fraction of rhodopsin molecules with different phosphorylation levels as a function of average phosphorylation, which was the only parameter actually measured [32]. The authors calculations were based on an additional assumption that unphosphorylated rhodopsin does not bind arrestin, even though specific low affinity binding of wild type arrestin to light-activated unphosphorylated rhodopsin in vitro [40], [41], and its role in inactivation of unphosphorylated rhodopsin in vivo [12], [23], [42] was shown. Arrestin binding in this study was measured using “extra Meta II” assay developed by Schleicher et al in 1989 [44]. This assay is based on the stabilization of active Metarhodopsin II state by bound arrestin. The most significant drawback of this assay is that it does not work above. At physiological temperatures extra Meta II is not detectable, even though it is obvious that arrestin effectively quenches rhodopsin signaling in mammals at –. In another study Vishnivetskiy et al [20] separated rhodopsin species with different levels of phosphorylation by chromatofocusing. Importantly, the authors quantitatively determined the presence of particular phospho-rhodopsin species in each fraction by mass-spectrometry of proteolytically removed rhodopsin C-terminus [34], obviating the need for calculations based on untested assumptions. Moreover, the binding assay in this study was performed at physiological temperature, . Based on their data, Vishnivetskiy et al concluded that arrestin demonstrates the same low affinity for rhodopsin carrying zero and one phosphate. The presence of two phosphates somewhat increases arrestin affinity, whereas arrestin binds rhodopsin with three, four, five, and six phosphates with the same high affinity, forming physiologically relevant long-lived complexes with stability sufficient for reliable quenching without possibility of reactivation on the time scale of the photoresponse [20]. These conclusions are in remarkable agreement with the work of Mendez et al in genetically modified mice expressing rhodopsin with different number of phosphorylation sites [3]. The authors of this study found that in vivo light activation of rhodopsin carrying zero, one, or two phosphorylation sites yields responses that last for many seconds, whereas rhodopsin carrying three or more phosphorylation sites is inactivated by wild type arrestin with sub-second kinetics [3]. Therefore, we based the modeling on the conclusions of these two studies [3], [20]. The parameters and appearing in Eq: 7–Eq: 8 are calibrated by WT and mutant experimental data, and not by fitting (see section on parameters). In the state, maintains its catalytic activity for a random time, until further phosphorylation or Arr binding. The sojourn times, for are exponentially distributed random variables with mean. The average lifetime of being deactivated after random biochemical states visited by before quenching, and is the sum of the up to. Thus (9) Hence and are determined by the biochemistry of the process through the sequences and. The number of steps after which binds to Arr is itself a random variable. The probability of shutoff in steps, or equivalently the probability of undergoing steps of phosphorylation and a final step for Arr binding, is (10) The mean steps of shutoff, or equivalently the expected value of the first moment of, is denoted by and is given by (11) Thus is the mean of the discrete valued random variable. The lifetime of is itself a random variable with expected value (12) These remarks permit one to detect the pattern of the means of the random sojourn times. First, the expected lifetime of, as a function of the number of available phosphorylation sites, decreases with increasing; second, the sojourn times, while increasing in number, each have a smaller mean. For example for 3P (), from Eq: 7–Eq: 9, and Eq: 10–Eq: 11, one computes (13) We stress that the number of available phosphorylation sites does not coincide with the mean number of steps to quenching. The number is fixed by the structure of rhodopsin, whereas depends on the biochemistry, through the probabilities. The dynamics of during deactivation is analyzed by the CTMC model in Eq: 5, which is numerically integrated in the Matlab platform. Its output is integrated into the spatio-temporal model and its Matlab code introduced in [8], [11]. This produces pointwise values of the effector, and on the ROS and permit one to compute the current response as a function of time and thus the functionals of effector and current in Eq: 1–Eq: 2. Simulations were performed for each of the 3 Test Cases. Random numbers are generated according to the exponential distribution of the random sojourn times, for each of wild type, transgenic and Arr knock-out cases indicated above. The corresponding dynamics is computed and the functionals in Eq: 1–Eq: 2 are evaluated. For each case, after about 5,000 numerical simulations, we compute the mean, the standard deviation and the coefficient of variation (CV) of these functionals. For WT mice the sequences and are chosen as in Eq: 6–Eq: 7. For mice lacking phosphorylation whether by COOH-terminal truncations ([24], [25]) or RK knockout ([1], [3]), in the CTMC one sets in Eq: 5 and, so that would remain in state 1 with catalytic activity for the whole process. If of the six phosphorylation sites are mutated out, we would have a CTMC model with. For mice lacking Arr ([24]), we let, for in the CTMC model in Eq: 5. In [37], we have generated a complete, self-consistent set of parameters for the mouse rod phototransduction, calibrated by least square fitting of the model in [8], [9], [11]. to a set of experimental data kindly provided Dr. C. Makino, leading to Table S2 (Text S1). Note that these parameters describe SPR recorded in Locke' s solution used in [3], [7], [24], [25], [28], [29], [42], [45]–[49], rather than in Ames' solution used in [1], [2]. For reasons that remain to be elucidated, the latter has greater amplitude and duration, although the recovery in both conditions is rate-limited by transducin inactivation ([28]). Figure 4 compares the simulated SPR by our model with the parameters of Table S2, and the experimental SPR kindly provided by C. Makino. The new parameters involved in the present investigation are the biochemical sequences, and. Below we indicate in detail how they have been determined. Their estimated values are reported in Table 4. Given the catalytic rate the sequence is determined from Eq: 6 whence the rate is known.
summary: Reception and transmission of biological stimuli such as vision, olfaction, taste, and hormone and neurotransmitter signal transduction, contain inherently variable components. Yet, biological functions are stable and reliable. For each signaling process, it is of interest to investigate the causes of variability and the mechanisms by which variability is mitigated to yield responses that reliably reflect the strength of the stimulus. We have investigated the variability of the single photon response in rod photoreceptors. A photon of light is captured by a receptor rhodopsin, and it goes through a series of biochemical states ending with a random shutoff. We have created a mathematical model of such a process, based on the recent biochemical findings on activation/deactivation, capable of reproducing the peculiar experimental features of visual trasduction both in wild type and genetically modified mice. We have found that the randomness of the time that rhodopsin sojourns in each of these biochemical states is the dominant cause of variability, whereas diffusion of molecules carrying the signal within the cell acts as variability mitigators.
keywords: physiology/sensory systems
  mathematics
  biophysics/theory and simulation
  physiology/cell signaling
  computational biology/systems biology
title: Kinetics of Rhodopsin Deactivation and Its Role in Regulating Recovery and Reproducibility of Rod Photoresponse