 Research article
 Open Access
 Open Peer Review
 Published:
Modeling and simulation of phototransduction cascade in vertebrate rod photoreceptors
BMC Ophthalmology volume 19, Article number: 55 (2019)
Abstract
Background
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.
Methods
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.
Results
The developed model structure could fit experimental data with small error.
Conclusions
The result indicated that the developed model structure could show the inactivation process of phototransduction cascades in the rod photoreceptors.
Background
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^{+}/Ca^{2+} ion channels in plasma membrane of rods, and leads to changes of electrical current and voltage [17]. 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. [24]. 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.
Method
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 11cis to the alltrans 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 k_{1} is the activation rate of rhodopsin.
Activation of Gprotein (GGDP) by activated rhodopsin
Through diffusion on the disc membrane, R* and inactivated GGDP interact, which results in a series of reactions including [12, 16]: (1) R* binds to GGDP and R*GGDP is formed, (2) GDP is released and R*G is formed, (3) GTP binds to R*G and R*GGTP is formed, (4) R* is released and GGTP 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 GGDP 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 MichaelisMenten kinetics, the intermediate reactions are simplified and these reactions can be approximately represented as:
where C_{1} is the intermediate complex formed by the binding of R* and G, and k_{2}, k_{3}, and k_{4} 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 [11]:
where k_{5} 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 [23]. According to MichaelisMenten kinetics, this process can be represented as:
where C_{2} is the intermediate complex formed by the binding of E* and cG, and k_{6}, k_{7}, and k_{8} 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 (G_{r}), and inactivation of both E* and G*. G_{r} will finally convert to GGDP (G) [1, 6]. These processes can be represented as:
where k_{9} is the rate of E* hydroxylation and k_{10} is the rate of G_{r} 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 Ca^{2+} influx, and thus cytoplasmic Ca^{2+} concentration is dropped. If Ca is used to denote the cytoplasmic Ca^{2+} concentration, the drop of Ca^{2+} concentration can be represented as:
where cG_{0} is the concentration of cGMP after darkadaption, cG is the current cGMP concentration, k_{11} is a constant, n_{cG} is the Hill coefficient that describes cGMP opening channel in a cooperative manner [1]. In this work, n_{cG} is set as 2 [12].
The reduced Ca^{2+} concentration causes recoverin (a guannylyl cyclase activation protein, GCAP) to release its Ca^{2+} and separate from rhodopsinkinase (RK). The increase of Ca^{2+} concentration is initiated by Ca^{2+} concentration drop from the darkadapted value Ca_{0} and is assumed as proportional to the concentration drop, which implies that the concentration change can be represented as:
where k_{12} 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 Ca^{2+} and originates from the reduction of Ca^{2+} concentration. For simplicity, the inactivation of R* is assumed as proportional to the Ca^{2+} concentration drop from its darkadapted value Ca_{0} and thus can be represented as:
where k_{13} is the deviation rate of R* mediated by the reduction of Ca^{2+} concentration.
The Ca^{2+}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 [3], 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 Ca^{2+}. For simplicity, the synthesis of cGMP is assumed as proportional to the Ca^{2+} concentration drop from its darkadapted value Ca_{0} and the drop of cG concentration from its darkadapted value cG_{0}. It is thus the change of cG concentration can be represented as:
where k_{14} is the synthesis rate of cGMP mediated by the reduction of Ca^{2+} 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 [23]. 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 k_{15} is a constant.
Independent of RK and Arr, R* may be inactivated in an abnormal and prolonged manner [12]. This process is represented as:
where k_{16} is the rate of R* deactivation that is not caused by the activities mediated by Ca^{2+} concentration change as in Eq. (9).
All major chemical reaction transitions are summarized in Fig. 1.
Let x_{l} through x_{8} denote the concentrations of R*, G*, C_{1}, E*, C_{2}, cG, Gr, and Ca^{2+} respectively. The concentration at the balance after darkadaptation for R, G, E, cG, Ca are denoted as R_{0}, G_{0}, E_{0}, cG_{0}, Ca_{0}, 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 (x_{6}) lead to closing or opening of cGMPgated ion channels and changes in the circulating current between the outer and inner segments of photoreceptor (Lam and Pugh, 2004; [16]). The deviation of current from the darkvalue is proportional to the changes of closing and opening of cGMPgated channels. The recorded electrical voltage can be represented as:
where f(t) is the measured ERG signal caused by light stimulation, k_{17} 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.).
Results
ERG data were reproduced from the publication of Lu et al. [16] and Lu [15]. ERG were recorded from three wildtype mice (61dayoldm males), three NOB1 mice (healthy, 61dayold), and three NOB1 mice with drug treatment (NmethyNnitrosourea). The bipolar cells of NOB1 mice were genetically disabled to suppress ERG bwave. Details about the experiments can be found in Lu et al. [16]. The LevenbergMarquardt method was used to estimate the model parameters (k_{1}k_{17}) in this work, which was adopted from the Matlab codes in the appendix of Lu [15]. In the LevenbergMarquardt method, the increment of each model parameter was computed according to the Jacobian matrix of partial derivatives of f(t) with respect to (k_{1}k_{17}) 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 LevenbergMarquardt method can be found in Marquardt [18], Levenberg [14], Constantinides and Mostoufi [4]. The RungeKutta algorithm was used to integrate the differentiation equations. Both the RungeKutta and the LevenbergMarquardt 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 R_{0}. It is thus the ratios of the state variables define the shape of ERG. The ratios for the initial values for R_{0}, G_{0}, E_{0}, cG_{0}, and Ca_{0} were set as 50, 5, 0.5, 4, and 0.22, respectively. The initial values for the state variables of x_{1} through x_{8} were set as [0 0 0 0 0 cG_{0} 0 Ca_{0}] since fully dark adaptation was applied before measurements. In this work, the values of total or dark concentrations (G_{0}, E_{0}, cG _{0}, and Ca_{0}) were also optimized in the algorithm while R_{0} 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 R_{0} 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 [16] were used to reduce ERG bwave, part of the bwave 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, photoreceptordamaged, and wildtype 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 k_{1u} to k_{17} 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, photoreceptordamaged NOB1, and wildtype), respectively. The simulated response from photoreceptors can match the initial part of ERGs and return the darkadapted 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 [21]. The fitting residuals in the three figures can serve as predictions of the ERG components from cells other than rod photoreceptors.
Discussions
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 bwave (generate later than the photoreceptor awave) at some degree makes ERG increase faster and the model structure does not cover the generation of ERG bwave, 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 wildtype subjects after the lowest trough of the ERGs. The fitting residual for wildtype group is stronger than the healthy NOB1 group, which suggests that the bwave 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 firstorder 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.
Conclusions
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, photoreceptor damaged NOB1, and Wildtype) 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.
Abbreviations
 cGMP:

cyclic guanosine monophosphate
 ERG:

Electroretinogram
 GGDP:

Gprotein
 PDE:

Phosphodiesterase
 R:

Rhodopsin
 R*:

activated rhodopsin
References
 1.
Arshavsky VY, Lamb TD, Pugh EN. G Proteins and Phototransduction. Annual Review of Physiology. 2002;64,153–87.
 2.
Armington JC. The electroretinogram (pp. 275–278). New York: Academic Press; 1974.
 3.
Biernbaum MS, Bownds MD. Lightinduced Changes in GTP and ATP in Frog Rod Photorecepters. Journal of General Physiology. 1985;85:107–121.
 4.
Constantinides A, Mostoufi N. Numerical methods for chemical engineers with MATLAB applications, vol. 443. Upper Saddle River, NJ: Prentice Hall PTR; 1999.
 5.
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.
 6.
Felber S, Breuer HP, Petruccione F, Honerkamp J, Hofmann KP. Stochastic simulation of the transducinGTPase cycle. Biophys J. 1996;71(6):3051–63.
 7.
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.
 8.
Karaśkiewicz J, DrobekSłowik M, Lubiński W. Pattern electroretinogram (PERG) in the early diagnosis of normaltension preperimetric glaucoma: a case report. Doc Ophthalmol. 2014;128(1):53–8.
 9.
Keener, J. & Sneyd, J. Mathematical Physiology II: Systems Physiology. New York: Springer; 2009.
 10.
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.
 11.
Lamb TD. Gain and kinetics of activation in the Gprotein cascade of phototransduction. Proc Natl Acad Sci. 1996;93(2):566–70.
 12.
Lamb TD, Pugh EN. A quantitative account of the activation steps involved in phototransduction in amphibian photoreceptors. J Physiol. 1992;449(1):719–58.
 13.
Lamb TD, Pugh EN Jr. Dark adaptation and the retinoid cycle of vision. Prog Retin Eye Res. 2004;23:307–80.
 14.
Levenberg K. A method for the solution of certain problems in least squares. Q Appl Math. 1944;2:164–8.
 15.
Lu L. Modeling of phototransduction in vision systems (master thesis, University of MissouriColumbia); 2007.
 16.
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.
 17.
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.
 18.
Marquardt DW. An algorithm for leastsquares estimation of nonlinear parameters. J. Soc. Ind. Appl. Math. 1963;11(2):431–41.
 19.
Marr BP, Singh AD. Retinoblastoma: evaluation and diagnosis. In: Clinical ophthalmic oncology: Springer Berlin Heidelberg; 2015. p. 1–11.
 20.
Nair SS, Paul Joseph K. Wavelet based electroretinographic signal analysis for diagnosis. Biomed Signal Process Control. 2014;9:37–44.
 21.
Nikonov S, Engheta N, Pugh EN. Kinetics of recovery of the darkadapted salamander rod photoresponse. J. Gen. Physiol. 1998;111(1):7–37.
 22.
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.
 23.
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.
 24.
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.
 25.
Sneyd J, Tranchina D. Phototransduction in cones: an inverse problem in enzyme kinetics. Bull Math Biol. 1989;51:749–84.
 26.
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.
 27.
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.
 28.
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.
Acknowledgements
Not applicable.
Funding
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.
Author information
Affiliations
Contributions
YG and JT conceived the idea. YG and PG did the modeling and data analysis. All the authors participated in the manuscript writing, and have read and approved the manuscript.
Corresponding author
Correspondence to Ya Guo.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Phototransduction
 ERG a wave
 Rod
 Photoreceptor
 Kinetic model structure