Modeling and simulation of phototransduction cascade in vertebrate rod photoreceptors
BMC Ophthalmology volume 19, Article number: 55 (2019)
The activation of phototransduction cascade in rod photoreceptors has been well studied in literature, but there is a lack of a mature kinetic model structure covering both the activation and inactivation processes.
In this work, a kinetic model structure is developed to describe the major activation and inactivation processes in vertebrate rod photoreceptors with the electroretinogram (ERG) as output. Simulation was performed to validate developed model structure.
The developed model structure could fit experimental data with small error.
The result indicated that the developed model structure could show the inactivation process of phototransduction cascades in the rod photoreceptors.
When the molecules of rod photoreceptors, rhodopsin, absorb photons, a cascade of reactions are triggered [10, 13]. This eventually results in the decline of intracellular cGMP concentration and closure of cGMP gated Na+/Ca2+ ion channels in plasma membrane of rods, and leads to changes of electrical current and voltage . The electrical current or voltage from the photoreceptors is an important part of electroretinogram (ERG), which is the summation of a few components of extracellular currents from different cells [13, 25, 28]. ERG has become an important diagnosis tool as the measurement standard protocols are developed and commercial instruments are available [8, 12, 19, 20, 22]. The response in the photoreceptor is very important to vision because it is the early stage of ERG and triggers latter vision activities [2, 5, 7, 16].
In literature, there are many efforts on modeling and simulation of phototransduction cascade in rods and ERG generation [6, 9, 11, 12, 16, 17, 21, 26, 27]. A recent review on modeling the retina in health, development and disease was well conducted by Roberts et al. . Phenomenological analytical functions or stochastic events are often used to simulate the processes. These efforts have significantly contributed to the understanding of phototransduction cascade and vision system. However, there are limitations if they are applied for the estimation of reaction rates from measured ERG signal because the reaction rates are not explicit variables in phenomenological model structure and there are many difficulties (e.g. computation speed) in estimating model parameters from stochastic event models. The models in literature usually only focus on the activation processes. Estimation of reaction rates from measured ERG signal can be very useful for diagnosis; therefore, it is meaningful to develop a kinetic model structure with reaction rates as parameters for both the activation and inactivation processes. A kinetic model structure for this purpose should be fundamentally based on major reaction kinetics but not include too many reaction details; otherwise, the complexity from the details will incur difficulty in parameter estimation algorithm convergence and reduce computation speed. Currently, there is a lack of a mature kinetic model structure covering both the activation and inactivation processes. In this work, a kinetic model structure is developed, which included the major reactions phototransduction cascade. The capability of the model structure to represent the inactivation process has been demonstrated.
Based on published work on the molecular mechanisms for phototransduction cascade in vertebrate rod photoreceptors, the involved major activation and inactivation biochemical reactions are summarized and the model structure is developed.
Activation of rhodopsin
A captured photon may set an inactivated rhodopsin (R) to its activated state (R*) by isomerizing its chromophore from 11-cis to the all-trans form [1, 12, 13]. The production rate of R* is proportional to the concentration of R and light intensity u, which can be represented by:
where k1 is the activation rate of rhodopsin.
Activation of G-protein (G-GDP) by activated rhodopsin
Through diffusion on the disc membrane, R* and inactivated G-GDP interact, which results in a series of reactions including [12, 16]: (1) R* binds to G-GDP and R*-G-GDP is formed, (2) GDP is released and R*-G is formed, (3) GTP binds to R*-G and R*-G-GTP is formed, (4) R* is released and G-GTP is formed, which then separates into submits Gβγ and Gα*-GTP. Gα*-GTP is activated and can trigger further reactions. R* is not altered in the processes and serves as an enzyme. For simplification, G is used to denote G-GDP and G* is used to denote the activated Gα*-GTP. It is not practical to include all the reaction details; otherwise, the complexity of the model structure will be dramatically increased. According to Michaelis-Menten kinetics, the intermediate reactions are simplified and these reactions can be approximately represented as:
where C1 is the intermediate complex formed by the binding of R* and G, and k2, k3, and k4 are the reaction rate constants.
Activation of phosphodiesterase (PDE) by activated G*
Two G* units bind to the two inhibitory subunits of inactivated PDE (denoted as E) and activate the α and β catalytic subunits of PDE and make PDE in active state [1, 12]; however, the two PDE subunits seem to work independently. If E* is used to denote a single activated subunit of PDE, the activation process of PDE can thus be represented as :
where k5 is the activation rate of E.
Change of cyclic GMP (cGMP) concentration
E* catalyzes the hydrolysis of cGMP (denoted as cG for conciseness) and produces GMP . According to Michaelis-Menten kinetics, this process can be represented as:
where C2 is the intermediate complex formed by the binding of E* and cG, and k6, k7, and k8 are the reaction rate constants.
Inactivation of E* and G*
The GTP bound to G* in the complex of E* is hydrolyzed to GDP, which results in separation of the complex, production of PDE and Gα-GDP (Gr), and inactivation of both E* and G*. Gr will finally convert to G-GDP (G) [1, 6]. These processes can be represented as:
where k9 is the rate of E* hydroxylation and k10 is the rate of Gr converts to G.
Channel activity, R* inactivation, and cGMP restoration
The decrease of cGMP concentration as represented in Eq. (4) leads to closure of cGMP gated channels and reduction of Ca2+ influx, and thus cytoplasmic Ca2+ concentration is dropped. If Ca is used to denote the cytoplasmic Ca2+ concentration, the drop of Ca2+ concentration can be represented as:
where cG0 is the concentration of cGMP after dark-adaption, cG is the current cGMP concentration, k11 is a constant, ncG is the Hill coefficient that describes cGMP opening channel in a cooperative manner . In this work, ncG is set as 2 .
The reduced Ca2+ concentration causes recoverin (a guannylyl cyclase activation protein, GCAP) to release its Ca2+ and separate from rhodopsinkinase (RK). The increase of Ca2+ concentration is initiated by Ca2+ concentration drop from the dark-adapted value Ca0 and is assumed as proportional to the concentration drop, which implies that the concentration change can be represented as:
where k12 is the rate that Ca is increased.
The free form of RK phosphorylizes R* and allows arrestin (Arr) to bind and inactivate R*. The inactivation of R* in this way is mediated by Ca2+ and originates from the reduction of Ca2+ concentration. For simplicity, the inactivation of R* is assumed as proportional to the Ca2+ concentration drop from its dark-adapted value Ca0 and thus can be represented as:
where k13 is the deviation rate of R* mediated by the reduction of Ca2+ concentration.
The Ca2+-free form of GCAP will bind to gyanylyl cyclase (GC) and turn on the enzymatic activity of GC, which catalyzes the synthesis of cGMP from GTP [16, 23]. According to Biernbaum and Bownds , GTP concentration change only occurs when light intensity is stronger than the saturation level and much later than the GC catalyzed cGMP production. In this work, the concentration of GTP is thus considered as a constant. Again, this activity is initiated and mediated by Ca2+. For simplicity, the synthesis of cGMP is assumed as proportional to the Ca2+ concentration drop from its dark-adapted value Ca0 and the drop of cG concentration from its dark-adapted value cG0. It is thus the change of cG concentration can be represented as:
where k14 is the synthesis rate of cGMP mediated by the reduction of Ca2+ concentration.
Reduction of cGMP concentration will also cause cGMP buffer in the cytoplasm releasing free form of cGMP and the diffusion of cGMP under concentration gradient . If these activities are assumed to proportional to the drop of cGMP concentration, the increase of free cGMP concentration in cytoplasm can be represented as:
where k15 is a constant.
Independent of RK and Arr, R* may be inactivated in an abnormal and prolonged manner . This process is represented as:
where k16 is the rate of R* deactivation that is not caused by the activities mediated by Ca2+ concentration change as in Eq. (9).
All major chemical reaction transitions are summarized in Fig. 1.
Let xl through x8 denote the concentrations of R*, G*, C1, E*, C2, cG, Gr, and Ca2+ respectively. The concentration at the balance after dark-adaptation for R, G, E, cG, Ca are denoted as R0, G0, E0, cG0, Ca0, respectively. According to chemical reaction kinetics, reaction speed is usually proportional to species concentration or the probability that reactants meet each other, species concentration changes rate with time represented by the phototransduction cascade in Eqns. (1) through (12) can be written as:
Changes of cGMP concentration (x6) lead to closing or opening of cGMP-gated ion channels and changes in the circulating current between the outer and inner segments of photoreceptor (Lam and Pugh, 2004; ). The deviation of current from the dark-value is proportional to the changes of closing and opening of cGMP-gated channels. The recorded electrical voltage can be represented as:
where f(t) is the measured ERG signal caused by light stimulation, k17 is a gain factor accounting for the number of photoreceptors and instrumentation gain (Because the produced current is proportional to the number of photoreceptors and instrumentation gain, the two factors will appear as product, which is merged into one term and cannot be separated.).
ERG data were reproduced from the publication of Lu et al.  and Lu . ERG were recorded from three wild-type mice (61-day-oldm males), three NOB1 mice (healthy, 61-day-old), and three NOB1 mice with drug treatment (N-methy-N-nitrosourea). The bipolar cells of NOB1 mice were genetically disabled to suppress ERG b-wave. Details about the experiments can be found in Lu et al. . The Levenberg-Marquardt method was used to estimate the model parameters (k1-k17) in this work, which was adopted from the Matlab codes in the appendix of Lu . In the Levenberg-Marquardt method, the increment of each model parameter was computed according to the Jacobian matrix of partial derivatives of f(t) with respect to (k1-k17) evaluated at all data points, the error between the model prediction and experimental data, and a damping factor for improving algorithm convergence. The algorithm of the Levenberg-Marquardt method can be found in Marquardt , Levenberg , Constantinides and Mostoufi . The Runge-Kutta algorithm was used to integrate the differentiation equations. Both the Runge-Kutta and the Levenberg-Marquardt were programmed in Matlab. The model parameters and the total or dark concentrations in Eqns. (13) through (21) can always be normalized by redefining the state variables and maintaining a constant R0. It is thus the ratios of the state variables define the shape of ERG. The ratios for the initial values for R0, G0, E0, cG0, and Ca0 were set as 50, 5, 0.5, 4, and 0.22, respectively. The initial values for the state variables of x1 through x8 were set as [0 0 0 0 0 cG0 0 Ca0] since fully dark adaptation was applied before measurements. In this work, the values of total or dark concentrations (G0, E0, cG 0, and Ca0) were also optimized in the algorithm while R0 was maintained as a constant for all the groups. By doing this, the number of unknown parameters was reduced by 1, which was significant to improving computation speed and algorithm convergence. Accordingly, the estimated model parameters were effective rates. The constant R0 value served as a reference and made the estimated effective rates comparable and could be used for classification as demonstrated in this work.
Although NOB1mouse and a chemical technique  were used to reduce ERG b-wave, part of the b-wave and/or other ERG components from cells other than rod photoreceptors still existed. Because the response of photoreceptors is the early part of ERG, only the initial data segment from the beginning of light stimulus to a moment slightly after ERG trough was used to estimate model parameters. In this way, the effect of ERG components from other cells was reduced since the proposed model structure was only for the phototransduction cascade in rod photoreceptors.
Figures 2, 3, and 4 show the comparisons of the averaged value of measured ERG data with model predictions and fitting residual for healthy NOB1, photoreceptor-damaged, and wild-type mice till 0.28 s, respectively. As shown in these figures, the model prediction could largely fit the initial stage of ERG, which is mainly from the rod photoreceptors. The relative fitting error for each of the three individual groups is less than 10% and the averaged value of relative errors for all the three groups of mice is 6.8%.
The mean value of k1u to k17 for the three groups of mice are shown in Table 1.
In order to verify the capability of the model structure to represent the inactivation process of rod phototransduction cascade, the parameters estimated from the initial part of the ERG data were used to simulate the photoreceptor response for a much longer time. Figures 5, 6, and 7 show the simulation results, measured ERG data, and the difference between the simulated results and ERG data (fitting residual) for the three types of mice (healthy NOB1, photoreceptor-damaged NOB1, and wild-type), respectively. The simulated response from photoreceptors can match the initial part of ERGs and return the dark-adapted value after the stimulation light is turned off for a while, which agree with the long term behavior of rods and ERGs without the responses from other cells like bipolar cell . The fitting residuals in the three figures can serve as predictions of the ERG components from cells other than rod photoreceptors.
Observations of Figs. 2, 3 and 4 reveal that the fitting errors mainly come from two sources. The first one is the high frequency noise on top of the signal, the other one is the discrepancy between experimental data and model prediction after the lowest trough of ERG. Because the existence of b-wave (generate later than the photoreceptor a-wave) at some degree makes ERG increase faster and the model structure does not cover the generation of ERG b-wave, the fitting error is larger at the later stage (after 0.2 s). The experimental data are bigger than model prediction for healthy NOB1 and wild-type subjects after the lowest trough of the ERGs. The fitting residual for wild-type group is stronger than the healthy NOB1 group, which suggests that the b-wave in ERG of the wild type group has more weight.
It is worth mentioning that the estimated reaction rates are effective rates because the initial concentrations were normalized. They are not comparable to the commonly reported pseudo first-order rates, but this does not affect the capability of the model structure to fit experimental data. In future research, further experiments including more subjects under other experimental conditions are needed to validate the developed model structure and compared model parameter changes with retina diseases.
In this work, a kinetic model structure is developed to describe the major reactions of activation and inactivation of phototrasduction cascade in rods. Measured ERGs from three groups of mice (NOB1, photo-receptor damaged NOB1, and Wild-type) are used to validate the model structure. The developed model structure could fit experimental data with small error. The developed model structure can represent both the activation and inactivation activities in rod phototransductions. In future research, more effort is needed to validate and establish the developed model structure.
cyclic guanosine monophosphate
Arshavsky VY, Lamb TD, Pugh EN. G Proteins and Phototransduction. Annual Review of Physiology. 2002;64,153–87.
Armington JC. The electroretinogram (pp. 275–278). New York: Academic Press; 1974.
Biernbaum MS, Bownds MD. Light-induced Changes in GTP and ATP in Frog Rod Photorecepters. Journal of General Physiology. 1985;85:107–121.
Constantinides A, Mostoufi N. Numerical methods for chemical engineers with MATLAB applications, vol. 443. Upper Saddle River, NJ: Prentice Hall PTR; 1999.
Einthoven W, Jolly WA. The form and magnitude of the electrical response of the eye to stimulation by light at various intensities. Exp Physiol. 1908;1(4):373–416.
Felber S, Breuer HP, Petruccione F, Honerkamp J, Hofmann KP. Stochastic simulation of the transducinGTPase cycle. Biophys J. 1996;71(6):3051–63.
Granit R. The components of the retinal action potential in mammals and their relation to the discharge in the optic nerve. J Physiol. 1933;77(3):207–39.
Karaśkiewicz J, Drobek-Słowik M, Lubiński W. Pattern electroretinogram (PERG) in the early diagnosis of normal-tension preperimetric glaucoma: a case report. Doc Ophthalmol. 2014;128(1):53–8.
Keener, J. & Sneyd, J. Mathematical Physiology II: Systems Physiology. New York: Springer; 2009.
Laitko U, Hofmann KP. A model for the recovery kinetics of rod phototransduction, based on the enzymatic deactivation of rhodopsin. Biophys J. 1998;74(2):803–15.
Lamb TD. Gain and kinetics of activation in the G-protein cascade of phototransduction. Proc Natl Acad Sci. 1996;93(2):566–70.
Lamb TD, Pugh EN. A quantitative account of the activation steps involved in phototransduction in amphibian photoreceptors. J Physiol. 1992;449(1):719–58.
Lamb TD, Pugh EN Jr. Dark adaptation and the retinoid cycle of vision. Prog Retin Eye Res. 2004;23:307–80.
Levenberg K. A method for the solution of certain problems in least squares. Q Appl Math. 1944;2:164–8.
Lu L. Modeling of phototransduction in vision systems (master thesis, University of Missouri--Columbia); 2007.
Lu L, Lei B, Guo Y, Tan J. Modeling of the phototransduction cascade in vertebrate rod photoreceptor. In Biomedical Engineering and Informatics (BMEI), 2010 3rd International Conference on (Vol. 3, pp. 1221–1228). Yantai: IEEE; 2010.
Mahroo OA, Ban VS, Bussmann BM, Copley HC, Hammond CJ, Lamb TD. Modelling the initial phase of the human rod photoreceptor response to the onset of steady illumination. Doc Ophthalmol. 2012;124(2):125–31.
Marquardt DW. An algorithm for least-squares estimation of nonlinear parameters. J. Soc. Ind. Appl. Math. 1963;11(2):431–41.
Marr BP, Singh AD. Retinoblastoma: evaluation and diagnosis. In: Clinical ophthalmic oncology: Springer Berlin Heidelberg; 2015. p. 1–11.
Nair SS, Paul Joseph K. Wavelet based electroretinographic signal analysis for diagnosis. Biomed Signal Process Control. 2014;9:37–44.
Nikonov S, Engheta N, Pugh EN. Kinetics of recovery of the dark-adapted salamander rod photoresponse. J. Gen. Physiol. 1998;111(1):7–37.
Praidou A, Hagan R, Newman W, Chandna A. Early diagnosis of Stargardt disease with multifocal electroretinogram in children. Int Ophthalmol. 2014;34(3):613–21.
Pugh EN Jr., Lamb TD. Phototransduction and Vertebrate Rods and Cones: Molecular Mechanisms of Amplification, Recovery and Light Adaptation. Molecular Mechanisms in Visual Transduction. ed. Stavenga WJ, de Grip WJ, Pugh EN Jr., Amsterdam; New York; Elsevier. 2000;183–255.
Roberts PA, Gaffney EA, Luthert PJ, Foss AJE, Byrne HM. Mathematical and computational models of the retina in health, development and disease. Prog Retin Eye Res. 2016;53:48–69.
Sneyd J, Tranchina D. Phototransduction in cones: an inverse problem in enzyme kinetics. Bull Math Biol. 1989;51:749–84.
Song Z, Coca D, Billings S, Postma M, Hardie RC. Biophysical Modeling of a Drosophila Photoreceptor, Neural Information Processing. Eds. Juusola M, Leung CS, Lee M, Chan JH. Berlin: Springer. 2009;5863:57–71.
Song Z, Postma M, Billings S, Coca D, Hardie R, Juusola M. Stochastic, adaptive sampling of information by microvilli in Fly photoreceptors. Curr Biol. 2012;22:1371–80.
Tranchina D, Sneyd J, Cadenas ID. Light adaptation in turtle cones. Testing and analysis of a model for phototransduction. Biophys J. 1991;60:217–37.
This project is partially supported by the National Natural Science Foundation of China (No: 31771680, used in data analysis), the Fundamental Research Funds for the Central Universities of China (No: JUSRP51730A, used in data interpretation), the 111 Project (B12018, used in data interpretation), and the Research Funds for New Faculty in Jiangnan University (used in manuscript writing).
Availability of data and materials
All the data supporting our findings are available through email request from the corresponding author.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Pan, G., Tan, J. & Guo, Y. Modeling and simulation of phototransduction cascade in vertebrate rod photoreceptors. BMC Ophthalmol 19, 55 (2019). https://doi.org/10.1186/s12886-019-1048-7