Multi-scale, numerical modeling of spatio-temporal signaling in cone phototransduction

Mammals have two types of photoreceptors, rods and cones. While rods are exceptionally sensitive and mediate vision at very low illumination levels, cones operate in daylight and are responsible for the bulk of visual perception in most diurnal animals, including humans. Yet the mechanisms of phototransduction in cones is understudied, largely due to unavailability of pure cone outer segment (COS) preparations. Here we present a novel mathematical model of cone phototransduction that explicitly takes into account complex cone geometry and its multiple physical scales, faithfully reproduces features of the cone response, and is orders of magnitude more efficient than the standard 3D diffusion model. This is accomplished through the mathematical techniques of homogenization and concentrated capacity. The homogenized model is then computationally implemented by finite element method. This homogenized model permits one to analyze the effects of COS geometry on visual transduction and lends itself to performing large numbers of numerical trials, as required for parameter analysis and the stochasticity of rod and cone signal transduction. Agreement between the nonhomogenized, (i.e., standard 3D), and homogenized diffusion models is reported along with their simulation times and memory costs. Virtual expression of rod biochemistry on cone morphology is also presented for understanding some of the characteristic differences between rods and cones. These simulations evidence that 3D cone morphology and ion channel localization contribute to biphasic flash response, i.e undershoot. The 3D nonhomogenized and homogenized models are contrasted with more traditional and coarser well-stirred and 1D longitudinal diffusion models. The latter are single-scale and do not explicitly account for the multi-scale geometry of the COS, unlike the 3D homogenized model. We show that simpler models exaggerate the magnitude of the current suppression, yield accelerated time to peak, and do not predict the local concentration of cGMP at the ionic channels.


The Homogenized Model in Weak Form
The homogenized weak form for cGMP diffusion is given by Here, where [E * ] is generated by the G-protein activation cascade, The homogenized weak form for Ca 2+ diffusion is given by Here S ∇ S ϕ · ∇ Sŵ dS = [λ(z)] 2 r 2 ϕ θŵθ + cos 2 (γ)ϕ zŵz λ(z)rsec(γ)dθdz (4) while also Note that J sat ex , J max cG are defined in (7), that from (4) the principal part in the sliver is a Laplace-Beltrami driven scaled by the normal width of the sliver, and that the test function class consists of all ϕ which are C ∞ smooth in the space-time variables.

GC Activity.
The synthesis rate of cGMP follows a spatially localized Hill-type law owing to the binding of Ca 2+ to GCAPS: The quantities α min and α max are respectively the least and greatest rate of cGMP synthesis by GC. m cyc is the Hill coefficient, and K cyc is the concentration for the half-maximal rate. To estimate α max , the concentration of guanylate cyclase in carp cones and carp rods respectively has been experimentally measured in (1) to be 72µM and 4.2µM respectively. Their ratio is then used to scale the (2) reported mouse rod value of α max = 76.5µM s −1 . The value of α max /α min is taken from mouse in (2).

Dark Hydrolysis.
The synthesis rate of cGMP in the dark is proportional with the dark concentration of cGMP through the proportionality constant β dark . The value β dark = 67s −1 was chosen and used in models by appealing to mass balance principles. It achieved equilibrium balance between PDE hydrolysis and cyclase synthesis for the given choice of other model parameters:

Buffers:
For B cG , in (3) the total quantities of PDE in carp cone and rods were reported as similar, so that their ratio is taken as 1. The buffering power for cGMP has been estimated as that ratio, 1, times the mouse rod value B cG = 1 reported in (2). For B Ca , (4) argues that the calcium buffer mechanism is more complicated in cones and varies as a function of calcium. Text equation (1.9) and Table 3 there propose a model and parameterize this dependence. At dark concentrations, [Ca 2+ ] dark = 0.4µM , which leads to an estimate B Ca = 20.01. This is consistent with (2)'s value of 20 in mouse rod.

Coupling Coefficient
The value c T E is the ratio ν RE /ν RG . The mouse rod value, 1, in (2) is taken.

Dark Steady State.
The concentration of cGMP in the dark was estimated in Carp (1) while measured by a model from samples of experimental dark current values in Striped Bass (4,5). Calcium concentration in the dark was reported for striped bass in (4)

Diffusion Coefficients
The values D cG − D R are taken as they are for rods in (2).

Volume to Surface Ratio
η is the asymptotic conversion ratio to pass from the volumic density defined in a chamber and the surface area densities defined on discal faces (2):

Fraction of Current Carried by Ca 2+
(8) has found the fraction of current carried by calcium in cones to be .33 and larger than the .06 value reported for rods (2).

Circulating Dark Current
In the homogenized model, dark current is not a free parameter but is determined by other parameter choices. Present simulations report a value J dark = 14.95pA. (See the discussion of K ex ).

Maximum cGMP-gated Current and Exchanger Current
The functional form of the currents is given by local Michaelis-Menten and Hill Laws (9). The exchanger current density and cGMP-gated channel current density are given respectively by The current values J sat ex , J max cGMP are the maximum currents measured across the whole COS, respectively for either the exchanger as [Ca 2+ ] becomes saturating or the cGMP-gated current as [cG] becomes saturating. Σ cone is the surface area of the cone at the sliver. This normalization assumes that the channels are distributed uniformly on the sliver.
In (5) and (10), the striped bass measurement for J max cG is reported. A range for J sat ex is also reported in (5). The upper value of this range is used in simulation.

Hydrolytic Efficiency of Activated PDE Dimer
In (11) Table 1, k cat /K M was reported to be (5 * 10 3 molecules/s/(10µM )). This value is within the range reported in (2) for mouse rod.

Surface Hydrolysis Rates of cGMP
The dark surface hydrolysis rate of cGMP, k σ,hyd , may be computed from Eq. 23 in (2): The light activated surface hydrolysis rate of cGMP, k * σ,hyd , may be computed from the expression just after Eq. 24 in (2): .02 * 10 23 /mol = .83µm 3 s −1

Decay of Rhodopsin
Rhodopsin activity is mediated through phosphorylation by rhodopsin kinase and arrestin binding. In (12), a continuous time markov chain (CTMC) framework is developed to account for the stochasticity of rhodopsin shutoff in rods. In the case of carp cones, it has been found that total kinase activity is much higher than in rods (13). In principle same CTMC framework may be used for both rods and cones while the parameters decribing cone opsin shut-off differ. The rate at which an opsin with i − 1 phosphorylations acquires an i th is λ i = (n step − i)λ 0 for i = 1, . . . n step .
The λ 0 value was estimated by scaling the mouse rod value in (12) by the ratio of carp cone GRK to carp rod GRK found in (14). The parameter n step counts the number of phosphorylations an opsin can undergo beginning at the first step with 0 phosphorylations. Arrestin binds an opsin in the i th phosphorylation state with rate µ i . For present simulations a single step to arrestin shutoff was chosen and µ 0 was taken as the mouse rod decay rate k R in (2). This was done due to a shortage of experimental measurements of arrestin binding phosphorylated cone opsin. For these choices the CTMC framework approximates an exponential decay model, which is the reasoning for substituting k R for µ 0 .

Catalytic Activity of Phosphorylated Opsin
Following (12) the activation rate of G-protein from opsin R * is assumed to decrease exponentially with incremental phosphorylation. This is described by the relation ν i = ν RG e −kv(i−1) for i = 1, . . . , n step .
The value k v for rods reported in (12) was used in simulations. The rate ν RG was taken from (11).

Michaelis-Menten and Hill Constants
K cyc is the Michaelis-Menten constant for cyclase appearing in equation (6), and the value reported in (5) is used. The value of m cGMP is also taken from (5). For K cG and m cyc , the mouse rod values reported in (2) were taken.
K ex is used as an adjusted quantity for fit in Table 3 of (5). There a value of 19nm was obtained. In our simulations K ex 's value of 0.69µM was chosen to fit the dark current of numerical simulations to the value 14.95pA ∼ 15pA. This K ex value is much closer to the mouse rod range of 0.9 − 1.6µM reported in (2) than the adjusted-to-fit value value reported in (5).

Geometric Constants
From the measurements of and ν = ν * ν , it follows that ν is taken as unity. The number of chambers may be computed from H, ν and through the relation n = H ν + = 15µm/(2 * .015µm) = 500

Effector Surface Density
The highest reported mouse rod value of (2) was taken in simulations.
To simulate the effects of reduced PDE inhibition, as in certain retinal rod disorders (15)(16)(17), the dark hydrolysis parameter β dark has been relatively increased and the effect on the cone's ten photon response shown. The numerical findings support the conclusion that disorders which increase the basal activity of PDE diminish the photoreceptor's sensitivity to light.

Biochemistry Drives Differences Between Rods and Cones in their HOM vs N-HOM's Drop-Relative-Errors
One observes that the relative errors in drop for the cone homogenized model are bigger than the relative errors in drop for the rod homogenized model. A back of the envelope calculation can explain this empirical observation.
The absolute error in drop between the homogenized and nonhomogenized models is given by the difference of the peak drops between HOM, drop h , and N-HOM, drop nh : e a = |drop h − drop nh |. Differences in rod and cone biochemistry ensure that the rod drop is much larger than the cone drop, and so, for example, even if the homogenized and nonhomogenized model attain the same absolute error on both rods and cones, still the cone relative error will be substantially larger than the rod relative error.

Verifying Text Eq (2)
The following Maple worksheet is to verify the proportion of vol ( C j ) in the cone. One minus this number is the value reported in Text Eq (2) The lowest height of the chamber C j is 1/2ν + (j − 1)(ν + ) The upper height of the chamber C j is 1/2ν + (j − 1)(ν + ) + .

1/3
The nonleading term is seen to be O( 2 ). Now simplify the leading term.