Sixth Mississippi State Conference on Differential Equations and Computational Simulations, Electronic Journal of Differential Equations, Conference 15 (2007), pp. 1–9. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu ftp ejde.math.txstate.edu (login: ftp)

MULTIPHOTON RESPONSE OF RETINAL ROD PHOTORECEPTORS

VASILIOS ALEXIADES, HARIHAR KHANAL

Abstract. Phototransduction is the process by which light is converted into an electrical response in retinal photoreceptors. Rod photoreceptors contain a stack of (about 1000) disc membranes packed with photopigment rhodopsin molecules, which absorb the photons. We present computational experiments which show the profound effect on the response of the distances (how many discs apart) photons happen to be absorbed at. This photon-distribution effect alone can account for much of the observed variability in response.

1. Introduction

Phototransduction is the process by which photons of light generate an electrical response in retinal rod and cone photoreceptors, thereby initiating vision. Phototransduction in rod photoreceptors is among the best understood biological signaling processes, with the underlying biochemistry, geometry, and physiology of the rod outer segment (ROS) known in fair detail. Thus, it is ripe for detailed quantitative modeling. Great strides were made over the past decade in developing basic models for the time–evolution of the response [11, 13, 12, 14, 15], by viewing the ROS as a single, homogeneous (“well-stirred”) compartment. More recently, spatio–temporal models accounting for diffusion of second messengers in the cytoplasm have been developed, [1, 2, 8, 9], focusing mainly on the single photon response (SPR) of the photoreceptor.

Here we consider multiple photon activation. We investigate how the response depends on the distribution of activation sites and identify those that yield maximal response. We find that the distance between ROS discs capturing photons can greatly contribute to response variability in dim light. Considerable nonlinear summation of SPR’s arises even with two photons.

After a brief description of the phototransduction process in §2, a spatio - temporal model (see [4]) is outlined, in dimensionless form, in §3. Simulations and their significance are described in §4, and conclusions in §5.

2000 Mathematics Subject Classification. 92C45, 35K60, 65M99.

Key words and phrases. Phototransduction; signaling; diffusion; parabolic system;

finite volume scheme; parallel computation.

\begin{matrix}c \\ ◯ \end{matrix}2007 Texas State University - San Marcos.

Published February 28, 2007.

2. Phototransduction

The first stage of vision occurs in photoreceptor cells in the back of the retina, which capture light and produce an electrical response. Rod photoreceptors contain a stack of (about 1000) “disc” membranes with embedded rhodopsin molecules, which absorb the photons and trigger a complex biochemical cascade ([7]) resulting in the production of activated phosphodiesterase (PDE^{∗}). This depletes cyclic guanosine monophosphate (cGMP) in the cytosol of the rod outer segment (ROS) of the photoreceptor. The plasma (lateral) membrane of the ROS contains cGMP– gated channels, which are open in darkness, permitting influx of Na^{+} and Ca^{2+} ions; a steady dark current is maintained by the Na^{+}/K^{+}/Ca^{2+} exchanger mechanism. Depletion of cGMP causes local closing of channels, thus lowering the local current across the plasma membrane. This is the signal that eventually reaches the brain enabling vision. A Ca^{2+}–mediated feedback mechanism deactivates rhodopsin and increases cGMP production, thus reopening the channels and restoring the dark current. We refer to [15] for a detailed description of phototransduction.

The essential players in the process are the effector, PDE^{∗}, produced on the discs, and the second messengers, cGMP and Ca^{2+}, which by diffusion in the cytosol carry the signal to and from the plasma membrane. The mathematical model assumes a specified number of activated PDE subunits in the entire ROS, E(t) (see Eq.(3.7) below) and accounts for the diffusion of cGMP and Ca^{2+} and their interactions on disc surfaces and on the plasma membrane. The current, J(t), across the plasma membrane can be found directly in terms of the boundary values of the concentrations [cG] and [Ca] at time t (§3). We are interested in the cellular response, defined as J_{dark}−J(t), or rather in the normalized response 1−J(t)/J_{dark}, with J_{dark} the dark steady-state current.

The model incorporates the main mechanisms presently known to operate in phototransduction, and reduces to simpler models proposed by physiologists. Namely, it reduces to the (one-dimensional) longitudinal model of [5] by assuming radially uniform concentrations, and to the lumped model ([13, 14, 15]) under the assumption of uniform (bulk) concentrations.

3. Mathematical Model

Geometry. The ROS of a rod photoreceptor in vertebrates can be considered as a right circular cylinder of height H and radius R_{rod}, housing a vertical stack of N equispaced parallel discs D_{j}, j = 1, 2, … N, coaxial with the cylinder, each of radius R_{dsc}, and thickness ε. The distance between discs, and the gap R_{rod} − R_{dsc} are also small, ∼ ε. Values for salamander ROS are given in §4.

The region inside the ROS not occupied by the discs is filled with cytosol. This is the region Ω where diffusion of second messengers takes place. We denote by F_{j}^{±} the upper/lower disc faces, and by ∂_{o}Ω the ROS lateral outer boundary (plasma membrane). The ratio of the cytosolic volume to the surface area of all disc faces will be denoted by η.

Dimensionless Form of the Model. Detailed formulation of the mathematical model, in physical variables, can be found in [4]. For brevity and clarity, here we present the model in dimensionless variables.

Let u and w denote dimensionless concentrations of the second messengers cGMP and Ca^{2+} respectively. We rescale lengths and time so that the geometric parameters and independent variables r, z, t, are all dimensionless. Employing cylindrical coordinates, the mathematical model for the diffusion of u and w in cytosol is expressed, in dimensionless form, as follows.

Given u(r, z, θ, 0) = u_{o}, w(r, z, θ, 0) = w_{o} with u_{o} > 0, w_{o} > 0, the initial uniform steady-state for the dark adapted system, find u(r, z, θ, t), w(r, z, θ, t) for t > 0, such that

\frac{∂u}{∂t} − ∇ · (D_{u} ∇u) = 0, \frac{∂w}{∂t} − ∇ · (D_{w} ∇w) = 0, in Ω, for t > 0, (3.1)

where D_{u} and D_{w} are the respective (dimensionless) diffusion coefficients.

Consider a beam of photons hitting a disc D_{j_{∗}} on one of its faces, say for example the lower one, F_{j_{∗}}^{−}, at coordinate z_{∗} along the axis of the rod. Generation and removal of free cGMP in the cytoplasm occurs through binding phenomena on the upper and lower faces F_{j}^{±} of each disc D_{j}. Calcium enters or leaves the cytosol only through the plasma membrane ∂_{o}Ω, (via the cGMP-gated channels and the electrogenic exchanger). Thus the two parabolic partial differential equations in (3.1) are coupled weakly via the following nonlinear boundary conditions:

−D_{u} \frac{∂u}{∂z} = η [±(C_{min} + C_{1}f_{1}(w)) ∓ C_{2}u] + δ_{j}C_{0}P u on F_{j}^{±}, j = 1, …, N, t > 0, (3.2)

−D_{w} \frac{∂w}{∂r} = g_{1}(w) − Cg_{2}(u) on ∂_{o}Ω, t > 0,  (3.3)

where δ_{j} = 1 if j = j_{∗} (activated face), and zero otherwise. C, C_{min}, C_{1}, C_{2} and C_{0} are positive constants expressing various interaction rates (synthesis or hydrolysis of cGMP, etc, see [4]). The quantity P represents the strength of PDE* - cGMP interaction, and thus the effect of activation by light (see Activation Mechanism below). The fluxes on the remaining parts of the boundary of Ω are zero. The functions f_{1}, g_{1} and g_{2} are given by

f_{1}(w) = \frac{1}{1 + (γw)^{m_{Ca}}} , g_{1}(w) = \frac{w}{1 + w} , g_{2}(u) = \frac{u^{m_{cG}}}{1 + u^{m_{cG}}} ,  (3.4)

where γ is a constant related to the channel opening and the cyclace rate, and m_{Ca}, m_{cG} are Hill constants (6.5 ≤ γ ≤ 16, m_{Ca} ≈ 2, m_{cG} = 2).

The local current J at a point of the plasma membrane (with local concentrations u, w at that point) is the sum of the cG-gated, J_{cG}, and exchanger, J_{ex}, circulating currents ([14, 15])

J = J_{cG} + J_{ex} , with J_{cG} = j_{cG}^{max}g_{2}(u) and J_{ex} = j_{ex}^{sat}g_{1}(w) ,  (3.5)

where j_{cG}^{max} is the maximal cG-gated current, and j_{ex}^{sat} is the saturation exchanger current.

Activation Mechanism. Light activation is embodied in the term δ_{j} C_{0} P u, appearing in Eq. (3.3). A satisfactory full modeling of the function P (x, t), for x ranging over a face F_{j_{∗}}^{−} hit by a photon, is an open problem under investigation. The literature contains various attempts to describe such a quantity (see [15, 10, 13, 14, 6, 1, 8, 9, 7]. Here we consider a simple activation mechanism with a lumped model, by taking the surface density of activated PDE molecules as the total PDE^{∗} in the rod divided by the area A_{activ} of activated discs,

P (t) = \frac{1}{2}E^{∗}(t)/A_{activ}  (3.6)

where E^{∗}(t) is the number of γ-subunits of PDE at time t in the entire rod (PDE is considered activated when both of its γ-subunits have been removed). Following [13, 15], the quantity E^{∗}(t) is approximated in terms of two first-order rate constants k_{R}, k_{E}, representing decay rates of activated rhodopsin R^{∗} and PDE^{∗}, as

E^{∗}(t) = Φ · (\frac{ν_{RE}}{k_{R} − k_{E}} ) (e^{−k_{E}t} − e^{−k_{R}t}) , t > 0 ,  (3.7)

where Φ is the number of photoisomerisations per rod per flash, and ν_{RE} is the effective rate with which a single R^{∗} triggers activation of PDE^{∗}.

Since this activation method applies the source P(t) uniformly on each activated disc, the process is axially symmetric, reducing the computation to two-dimensional (in r, z coordinates).

4. Numerical Simulations

Discretization and Parallelization. We employ Finite Volume discretization and explicit-implicit time-stepping, implemented in Fortran, with time-steps sufficiently small to ensure numerical stability of the scheme.

Due to the intricate geometry of the cytosol, the problem involves very intensive computations demanding high performance computing. This was addressed by parallelization, via domain decomposition, for clusters of distributed memory multiprocessors. The idea is to decompose the spatial domain into sections and assign a section to each processor. Here a section consists of a group of disc units. The parallel implementation employs the MPI (Message Passing Interface) library, following the master/slaves paradigm generated in SIMD (Single Instruction Multiple Data) mode, where one processor acts as a master and the rest as slaves. The master loads I/O, distributes tasks to the slaves, controls and synchronizes the slaves, whereas the slaves all solve the same problem but on their own segment of the mesh, exchange boundary values with their neighbors, and send their output to the master. A typical 1 sec simulation of the axisymmetric case, for a rod with 800 discs, with a fairly fine grid (34 radial nodes × 4*800 axial nodes, for a total of 108,800 nodes) takes about 4.5 hours on 11 processors of a linux cluster (AMD Opteron, with EKO pathscale f90 compiler).

Simulation Setup. Finding a consistent set of parameters for the model has been a major issue. We obtained some of the parameters from the literature ([14, 15, 10, 13, 11, 12]), some by matching terms in our model and the bulk model of [14, 15], and some by testing many combinations of parameters attempting to match the peak response, and the time at which it occurs, with experimental data of Rieke on single photon response in salamander rods. Details and parameter values appear in [4].

Simulations were performed for a typical salamander photoreceptor ROS with N = 800 discs, R_{dsc} = 5.5 μm, ε = 14 nm, R_{rod} = 5.515 μm and H = 22.4 μm. Typical parameters for the whole cell electrical properties at the dark resting state are taken to be j_{cG}^{max} = 7000 pA and j_{ex}^{sat} = 0.17 pA. The initial state is the dark steady-state with concentrations u_{0} ≡ u_{dark} = 3 μM, w_{0} ≡ w_{dark} = 0.66 μM, found by solving the system (3.2)-(3.3) after setting the fluxes to zero. The resulting dark current is J_{dark} = 66 pA. For the PDE-activation step, the lumped method described in §3 was employed with rate constants ν_{RE} = 195 s^{−1}, k_{E} = 0.6 s^{−1} and k_{R} = 2.6 s^{−1} , and various values of Φ as shown in the simulation plots below.

The numerical solution of the mathematical model gives the evolution of the spatial distribution of u(r, z, t) and w(r, z, t) in the cytosol Ω. The local circulating current J(z, t) at height z at time t is obtained from (3.5) using the boundary values (at r = R_{rod}) of u and w at that z and t. The current J(t) across the entire plasma membrane at time t is the mean value of J(z, t) over ∂_{o}Ω. Results are presented in terms of the normalized response 1−J/J_{dark} for either local or total J. To quantify the longitudinal spread, we count the number of discs for which response is more than 0.1% of the dark value.

Single Photon Response. To put the multiphoton responses in perspective, we begin with the case of activation by a single photon (Φ = 1 in Eq.(3.7)) delivered at the 400^{th} disc of a salamander rod photoreceptor, Fig. 1.

The peak response is 0.82% of the 66 pA dark current, occuring at 860 ms, which agrees with experimental data of F. Rieke ([4]). Fig. 1(B) shows the longitudinal profile of the local response (1−J(z, t)/J_{dark}); it peaks at ∼ 14%, and the maximal longitudinal spread is 311 disc units, occuring at 1300 ms.

We have verified that activating any single disc farther away than ∼ 25 discs from the ends of the rod outer segment produces identical behavior (except, of course, with z-profiles shifted at the activation site).

Multi-photon Response. Simulations with a higher intensity stimulus of Φ = 7 photoisomerizations are presented in Fig.2 for three different arrangements of activation sites: (a) all 7 photons on one disc (green curves), (b) one photon on each of seven adjacent discs (blue curves), (c) one photon every 100 discs (red curves). The peak response is only 2.1% in case (a), and 3.7% in case (b), much less than the sum of seven individual SPR’s; thus there is “nonlinear summation” of SPR’s in these cases. On the contrary, in case (c) the peak total response is 5.7%, and the local response at each of the activated discs peaks at 14%, the same as a single SPR; the discs, being sufficiently far apart, hardly interact and their contributions to the total response are additive, resulting in “linear summation” of SPR’s.

Proximity reduces longitudinal gradients, hence the total response, and thus maximal separation produces maximal total response. This phenomenon creates variability of responses to non-saturating light, apart from any other factors, merely due to where (how far apart) photons happen to be absorbed.

Indeed, for higher Φ the effect can be much more pronounced than seen in Fig. 2. For example, as shown in Fig. 3, whereas 700 photons all on one disc produce 4% peak response, they produce 38% peak when distributed on 12 discs lying 70 discs apart from each other, and 86% peak when spread one each on 700 consecutive discs. Considering the enormous number of possible ways that 700 photons could be distributed among 800 discs (of the order of 10^{25}), it is clear this photon-distribution effect alone can generate tremendous variability in response, essentially any response in the range between minimal and maximal (e.g. between 4% and 86% when Φ = 700). Minimal response is produced when all the photons are on one disc (adaptation takes over reducing the response), and maximal when Figure 1. Simulations of response to a single photon flash delivered at the 400^{th} disc. (A) Total response 1 − J(t)/J_{dark} versus time. (B) Local response 1 − J(z, t)/J_{dark} versus location z at times t = 100, 200, and 1000 ms.

they are maximally separated, with as few as possible on each disc (to minimize adaptation).

Nonlinear summation of SPR’s. In experiments, nonlinear summation of single photon responses becomes noticeable (about 10% deviation from linearity) for flashes that suppress 20-25% of the dark current. On the other hand, we already saw nonlinear summation arising for 7 photons above. To explore this further, we simulate the case of 2 photons placed at various distances (0, 1, 10, 20, 40, 50, 100, 200, 400 discs apart), centered about the middle (400^{th} disc) of the ROS. As seen in Fig.4, the response is maximal, at twice that of a single photon (linear summation), when the two activated discs are at least 200 discs apart. It begins to diminish at ∼ 150 discs apart, as their spreads begin to overlap, and reduces down to ∼ 90% of maximal when the activated discs are adjacent (or when both photons are on Figure 2. Multiphoton simulations with Φ = 7 photoisomerisations for three different arrangements of activation sites. (a) only one disc is activated (disc # 400) by 7 photons (green); (b) seven adjacent discs are activated (discs # 397 - 403), each by one photon (blue); (c) seven discs are activated, 100 discs apart from each other (disc # 100, 200, 300, 400, 500, 600, 700) (red). (A): Total response versus time. (B): z-profiles of local response.

Figure 3. Multiphoton simulations with Φ = 700 photoisomerisations for three different arrangements of activation sites, to illustrate the dramatic effect of photon distribution. (a) all 700 photons on one disc (disc # 400) (green); (b) on 12 discs located 70 discs apart (discs #15, 85, …, 785) (blue); (c) one photon on each of 700 consecutive discs (discs #50–749) (red). (A): Total response versus time. Peak response is only 4% in (a), but 38% in (b), and 86% in (c). Other arrangements can generate essentially any intermediate response. (B): z-profiles of local response at time 1000 ms for each of the three cases in (A).

the same disc); at that point, the 9.7% deviation from linear summation would become “noticeable”. Thus, our simulations show that nonlinear summation occurs at much lower light, even for 2 photons, and it is strictly due to how far apart the photons act upon. At present, such a setting cannot be achieved experimentally, since individual photons cannot be directed to impinge at specific locations.

Figure 4. Two-photon simulations (Φ = 2) exhibiting nonlinear summation of SPR’s as the distance between activated discs decreases. Left scale shows peak total response as % of J_{dark}, while the right scale shows deviation from linear summation as % of maximal relative response, vs distance in disc units. When the two photons are sufficiently far apart (at least 200 discs), the (total) response is maximal (1.64%) at twice the SPR (0.82%), thus their effect is additive. The deviation becomes 9.7% when the two photons act on adjacent discs or on the same disc, where the response is minimal (1.48%). The peak times range from 830 to 860 ms.

5. Conclusions

Employing a spatio-temporal model for rod phototransduction, that incorporates all the mechanisms presently known to operate in generating the photoreceptor response, we examined the effect of the spacing of activating sites on the photoreceptor response.

Activating discs far apart from each other produces considerably higher response than when the photons act on one or a few nearby discs. Our numerical simulations show that the effect can be much more significant than has hitherto been appreciated, in generating response variability at low light, apart from any other factor. This variability arises already with 2 photons, and can be very pronounced for a higher number of photons of non-saturating light. E.g. in the case of 700 photons, the response can be anything between 4% and 86% of the dark current, by merely varying the arrangement of activated discs.

Another way of viewing the phenomenon is that of additivity of individual responses. Proximity reduces longitudinal concentration gradients, resulting in lower total response than the sum of individual single photon responses. For 2 photons, deviation from linear summation begins when the photons act ∼ 150 discs apart and it reaches 10% when the photons are on the same or on adjacent discs.

The results presented here point out the usefulness of the detailed spatio - temporal model of the process in conducting (virtual) experiments that cannot at present be achieved experimentally.

Acknowledgments. This work was partially supported by NIH grant NIH-1-RO1-GM 68953-01, and by the University of Tennessee.

References

[1] D. Andreucci, P. Bisegna, G. Caruso, H.E. Hamm and E. DiBenedetto, Mathematical Model of the Spatio-Temporal Dynamics of Second Messengers in Visual Transduction, Biophysical J. 85: 1358–1376, 2003.

[2] D. Andreucci, P. Bisegna, and E DiBenedetto, Homogenization and Concentrated Capacity for a Problem in Visual Transduction, Annali di Mat. Pura et Appl. 182(4): 375–407, 2003.

[3] D. A. Baylor, T.D. Lamb and K.-W. Yau, The Membrane Current of Single Rod Outer Segments, J. Physiol., 288: 589–611, 1979.

[4] G. Caruso, H. Khanal, V. Alexiades, F. Rieke, H.E. Hamm, E. DiBenedetto, Mathematical and Numerical Modeling of Spatio-Temporal Signaling in Rod Phototransduction, IEE Proc. Systems Biology, 152(3): 119-137, 2005.

[5] M. Gray-Keller, W. Denk, B. Shraim and P.B. Detwiler, Longitudinal Spread of Second Messenger Signals in Isolated Rod Outer Segments of Lizards, J. Physiol., 519: 679–692, 1999.

[6] R.D. Hamer, Analysis of Ca^{++}-dependent gain changes in PDE activation in vertebrate rod phototransduction, Molecular Vision 6: 265–286, 2000.

[7] R.D .Hamer, S.C. Nicholas, D. Trachina, P.A. Liebman, and T.D. Lamb, Multiple steps of phosphorylation of activated rhodopsin can account for the reproducibility of vertebrate rod single-photon response, J. Gen. Physiol. 122: 419–444, 2003.

[8] H. Khanal, V. Alexiades, E. DiBenedetto and H. Hamm, Numerical Simulation of Diffusion of Second Messengers cGMP and Ca^{2}+ in Rod Photoreceptor Outer Segment of Vertebrates, pp. 165-172 in Unsolved Problems of Noise and Fluctuations in Physics, Biology and High Technology, editor Sergey Bezrukov, American Institute of Physics, 2003.

[9] H. Khanal, V. Alexiades, E. DiBenedetto, Response of Dark-adapted Retinal Rod Photoreceptors, pp. 138-145 in Dynamic Systems and Applications 4, editor M. Sambandham, Dynamic Publishers, 2004.

[10] Y Koutalos, K Nakatani, and K-W Yau, Cyclic GMP Diffusion Coefficients in Rod Photoreceptors Outer Segments, Biophysical J., 68: 373–382, 1995.

[11] T.D. Lamb and E.N. Pugh Jr., A quantitative account of the activation steps involved in phototransduction in amphibian photoreceptors, J. Physiol. 449: 719–758, 1992.

[12] I.B. Leskov, V.A. Klenchin, J.W. Handy, G.G. Whitelock, V.I. Govardovskii, M.D. Bownds, T.D. Lamb, E.N. Pugh, & V.Y. Arshavsky, The Gain of Rod Phototransduction: Reconciliation of Biochemical and Electrophysical Measurements, Neuron, 27: 525–537, 2000.

[13] S. Nikonov, N. Engheta and E.N. Pugh, Jr Kinetics of Recovery of the Dark-adapted Salamander Rod Photoresponse, J. Gen. Physiol, 111: 7–37, 1998.

[14] S .Nikonov, T.D. Lamb and E.N. Pugh, Jr. The Role of Steady Phosphodiestarase Activity in the Kinetics and Sensitivity of the Light-Adapted Salamander Rod Photoresponse, J. Gen. Physiol., 116: 795–824, 2000.

[15] E.N. Pugh, Jr. and T.D. Lamb, Phototransduction in Vertebrate Rods and Cones: Molecular Mechanisms of Amplification, Recovery and Light adaptation, pp.183–255 in Molecular Mechanism in Visual Transduction, edited by DG Stavenga, WJ Degrip & EN Pugh Jr, Elsevier, Amsterdam, 2000.

Vasilios Alexiades Department of Mathematics, University of Tennessee, Knoxville, TN 37996, USA and Oak Ridge National Laboratory, Oak Ridge TN 37831, USA

E-mail address: alexiades@utk.edu

Harihar Khanal Department of Mathematics, Embry-Riddle Aeronautical University, Daytona Beach, FL 32114, USA

E-mail address: Harihar.Khanal@erau.edu