\documentclass[reqno]{amsart}
\usepackage{graphicx}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
Sixth Mississippi State Conference on Differential Equations and
Computational Simulations,
{\em Electronic Journal of Differential Equations},
Conference 15 (2007),  pp. 1--9.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu  (login: ftp)}
\thanks{\copyright 2007 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2006/Conf/15/\hfil Multiphoton response]
{Multiphoton response of retinal rod photoreceptors}

\author[V. Alexiades, H. Khanal\hfil EJDE/Conf/15 \hfilneg]
{Vasilios Alexiades, Harihar Khanal}  % in alphabetical order

\address{Vasilios Alexiades \newline
Department of Mathematics, University of Tennessee,
 Knoxville, TN 37996, USA \newline
 and  Oak Ridge National Laboratory, Oak Ridge TN 37831, USA}
\email{alexiades@utk.edu}

\address{Harihar Khanal \newline
Department of Mathematics, Embry-Riddle Aeronautical University,
Daytona Beach, FL 32114, USA}
\email{Harihar.Khanal@erau.edu}

\thanks{Published February 28, 2007.}
\subjclass[2000]{92C45, 35K60, 65M99}
\keywords{Phototransduction; signaling; diffusion;
 parabolic system; \hfill\break\indent
finite volume scheme; parallel computation}

\begin{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.
\end{abstract}

\maketitle
\numberwithin{equation}{section}

\section{Introduction}\label{S:1}

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
\cite{LP92, NEP98, Leskov00, NLP, PL},
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,
\cite{Andreucci03a, Andreucci03b, Khanal03, Khanal04}, 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
\S \ref{S:2},
a spatio - temporal model (see \cite{caruso}) is outlined, in
dimensionless form, in \S\ref{S:3}.
Simulations and their significance are described in \S\ref{S:4},
and conclusions in \S\ref{S:5}.


\section{Phototransduction}\label{S:2}

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 (\cite{Hamer03})
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 {\it 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 \cite{PL} for a detailed description of phototransduction.

The essential players in the process are the effector, \text{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.(\ref{eq:6}) 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 $\text{[cG]}$ and $\text{[Ca]}$
at time $t$ (\S\ref{S:3}).
We  are interested in the cellular {\it response},
defined as $J_\text{dark} - J(t)$,
or rather in the {\it normalized response}
$1 - J(t)/J_\text{dark}$,
with $\text{J}_{\text{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 \cite{GrayKeller99} by
assuming radially uniform concentrations, and to
the lumped model (\cite{NEP98, NLP, PL})
under the assumption of uniform (bulk) concentrations.


\section{ Mathematical Model}\label{S:3}

\subsection*{Geometry}
The ROS of a rod photoreceptor in vertebrates can be considered
as a right circular cylinder of height $H$ and radius $R_{\rm rod}$,
housing a vertical stack of $N$ equispaced
parallel discs $\mathcal{D}_j$, $j=1,2,\dots N$, coaxial with the cylinder,
each of radius $R_{\rm dsc}$, and thickness $\varepsilon$.
The distance between discs, and the gap $R_{\rm rod} - R_{\rm dsc}$
are also small, $\sim \varepsilon$.
%
Values for salamander ROS are given in \S\ref{S:4}.

The region inside the ROS not occupied by the discs is filled with
cytosol.  This is the region $\Omega$
where diffusion of second messengers takes place.
We denote by  $F_j^{\pm}$ the upper/lower disc faces, and by
$\partial_o\Omega$ 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 $\eta$.



\subsection*{Dimensionless Form of the Model}
Detailed formulation of the mathematical model, in physical variables,
can be found in \cite{caruso}. 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, \theta, 0) = u_o$,
$w(r, z, \theta, 0) = w_o$ with $u_o>0$, $w_o>0$, the initial uniform
steady-state for the dark adapted system,
find $u(r, z, \theta, t)$,  $w(r, z, \theta, t)$ for
$t>0$, such that
%
\begin{equation}\label{eq:1}
\frac{\partial u}{\partial t}  -  \nabla\cdot (D_u\, \nabla u) = 0,
\quad
\frac{\partial w}{\partial t} - \nabla\cdot (D_w\, \nabla w) = 0,
\quad \text{in } \Omega, \quad \text{for } t>0,
\end{equation}
%
where $D_u$ and $D_w$ are the
respective (dimensionless) diffusion coefficients.

Consider a beam of photons hitting a disc $\mathcal{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}^{\pm}$
of each disc $\mathcal{D}_j$.
Calcium enters or leaves the cytosol only through the plasma membrane
$\partial_o\Omega$, (via the cGMP-gated channels
and the electrogenic exchanger).
Thus the two parabolic partial
differential equations in (\ref{eq:1}) are coupled weakly via the
following nonlinear boundary conditions:
%
\begin{gather}
-  D_{u}\frac{\partial u}{\partial z}
= \eta \left[ \pm  (C_\text{min} + C_1  f_1(w)) \mp
C_2  u  \right] + \delta_{j} C_0 P u
\quad \mbox{on  } F_j^{\pm}, \quad j=1,\dots,N, \; t>0,\label{eq:2a} \\
-  D_{w}\frac{\partial w}{\partial r}
= g_1(w)  -  C  g_2(u)
\quad \mbox{on  } \partial_o\Omega,  \; t>0, \label{eq:2b}
\end{gather}
where
$\delta_{j} =1$ if $j=j_{*}$ (activated face), and zero otherwise.
$C, \, C_\text{min},\, C_1, \, C_2$ and $C_0$
are positive constants
expressing various interaction rates
(synthesis or hydrolysis of cGMP, etc,
see \cite{caruso}).
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 $\Omega$ are zero.
The functions $f_1$, $g_1$ and $g_2$ are given by
%
\begin{equation}\label{eq:3}
f_1(w)=\frac{1}{1+(\gamma w)^{m_{\text{Ca}}}} \, ,
\quad
g_1(w)=\frac{w}{1+w} \, ,
\quad
g_2(u)= \frac{u^{m_{\text{cG}}}}{1+u^{m_{\text{cG}}}} \, ,
\end{equation}
%
where $\gamma$ is a constant related to the channel opening
and the cyclace rate, and
$m_{\text{Ca}}$, $m_{\text{cG}}$ are Hill constants
($6.5 \le \gamma\le 16$,
$m_{\text{Ca}} \approx 2$,
$m_{\text{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_{\text{cG}}$, and exchanger, $J_{\text{ex}}$,
circulating currents (\cite{NLP, PL})
%
\begin{equation}\label{eq:4}
J = J_{\text{cG}}+J_{\text{ex}} \, ,
\quad \mbox{with } \quad J_\text{cG}=
j_{\text{cG}}^{\max} g_2(u)
\quad \mbox{  and  } \quad
J_\text{ex} =
j_{\text{ex}}^{\text{sat}} g_1(w) \, ,
\end{equation}
%
where
$j_{\text{cG}}^{\text{max}}$ is the maximal cG-gated current,
and
$j_{\text{ex}}^{\text{sat}}$ is the saturation exchanger current.

\subsection*{Activation Mechanism}
Light activation is embodied
in the term $\delta_{j}\,C_0\,P\,u \,$,
appearing in Eq. (\ref{eq:2b}).
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 \cite{PL, Koutalos95, NEP98, NLP, Hamer00, Andreucci03a, Khanal03,
Khanal04, Hamer03}.
%
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_\text{activ}$
of activated discs,
%
\begin{equation}\label{eq:5}
P(t) =
{\tfrac12}  \text{E}^*(t)/ A_\text{activ}
\end{equation}
%
where $\text{E}^*(t)$ is the number of $\gamma$-subunits of PDE
at time $t$ in the entire rod (PDE is considered activated when
both of its $\gamma$-subunits have been removed).
%
Following  \cite{NEP98, PL},
the quantity E${}^*(t)$ is approximated in terms of
two first-order rate constants
$k_\text{R}$, $k_\text{E}$, representing decay rates of activated
rhodopsin R${}^*$ and $\text{PDE}^*$, as
%
\begin{equation}\label{eq:6}
\text{E}^* (t) = \Phi \cdot
\left(\frac{\nu_\text{RE}}{k_\text{R}-k_\text{E}}\right)
\left( e^{-k_\text{E}t} -
e^{-k_\text{R}t} \right), \quad t > 0 \, ,
\end{equation}
%
where $\Phi$ is the number of photoisomerisations per rod per flash,
and $\nu_\text{RE}$ is the effective rate with which
a single $\text{R}^*$ triggers activation of $\text{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).

\section{Numerical Simulations}\label{S:4}


\subsection*{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 $\times$ 4*800 axial nodes,
for a total of  108,800 nodes) takes about
%  $8$ hours on $21$ processors of cheetah (IBM pSeries System)
% $13$ hours on $4$ networked dual Intel Xeon 3 GHz workstations.
$4.5$ hours on $11$ processors of a linux cluster
(AMD Opteron, with EKO pathscale f90 compiler).


\subsection*{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 (\cite{NLP, PL, Koutalos95, NEP98, LP92, Leskov00}),
some by matching terms in our model and the bulk model of
\cite{NLP, PL}, 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 \cite{caruso}.

Simulations were performed for a typical {\it salamander}
photoreceptor ROS with
$N = 800$ discs,
$R_{\text{dsc}} = 5.5\,\mathrm{\mu m}$,
$\mathrm{\varepsilon} = 14\,\mathrm{nm}$,
$R_{\text{rod}} = 5.515\,\mathrm{\mu m}$ and
$H = 22.4\,\mathrm{\mu m}$.
Typical parameters for the whole cell electrical properties
at the dark resting state are taken to be
$j_{\text{cG}}^\text{max} = 7000\, \mathrm{pA}$ and
$j_{\text{ex}}^\text{sat} = 0.17\, \mathrm{pA}$.
The initial state is the dark steady-state with concentrations
$u_0 \equiv u_\text{dark} = 3\, \mathrm{\mu M}$,
$w_0 \equiv w_\text{dark} = 0.66\, \mathrm{\mu M}$,
found by solving the system \eqref{eq:2a}-\eqref{eq:2b} after setting
the fluxes to zero.
The resulting dark current is $\text{J}_{\text{dark}} = 66$ pA.
For the PDE-activation step,
the lumped method described in \S\ref{S:3} was employed
with rate constants $\nu_\text{RE} = 195 \,\, \mathrm{s^{-1}}$,
$k_\text{E} = 0.6 \,\, \mathrm{s^{-1}}$ and
$k_\text{R} = 2.6\,\, \mathrm{s^{-1}}\,$,
and various values of $\Phi$ 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 $\Omega$. The local circulating current $J(z,t)$ at
height $z$ at time $t$ is obtained  from (\ref{eq:4}) using
the boundary values (at $r=R_{\text{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 $\partial_o\Omega$.
Results are presented in terms of the
{\it normalized response}
$1-J/J_\text{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.



\subsection*{Single Photon Response}
To put the multiphoton responses in perspective, we begin with
the case of activation by a single photon
($\Phi = 1$ in Eq.(\ref{eq:6}))  delivered
at the $400^\text{th}$ disc of a salamander rod photoreceptor,
Fig. \ref{fig:p1}.

%*********************************  Fig1 a,b 	

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{figures/fig1a}
\includegraphics[width=0.7\textwidth]{figures/fig1b}
\end{center}
\caption{\label{fig:p1} 
Simulations of response to a {\bf single photon} 
flash delivered at the 
$400^\text{th}$ disc.
{\bf (A)} Total response $1-J(t)/J_\text{dark}$ versus time.
{\bf (B)} Local response $1-J(z,t)/J_\text{dark}$ versus 
location $z$ at times $t =$ 100, 200, and 1000 ms.}
\end{figure}

The peak response is $0.82\%$   %% $0.54$ pA
of the $66$ pA dark current, occuring at $860$ ms,
which agrees with experimental data of F. Rieke (\cite{caruso}).
Fig. \ref{fig:p1}(B) shows the longitudinal profile of the local
response ($1-J(z,t)/\text{J}_{\text{dark}}$); it peaks at $\sim 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 $\sim 25$
discs from the ends of the rod outer segment produces identical behavior
(except, of course, with $z$-profiles shifted at the activation site).



\subsection*{Multi-photon Response}
Simulations with a higher intensity stimulus
of $\Phi=7$ photoisomerizations are presented in Fig.\ref{fig:p7}
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 {\bf 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.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fig:p7    Fig 2

\begin{figure}[ht]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.47\textwidth]{figures/fig2a}
&\includegraphics[width=0.47\textwidth]{figures/fig2b}
\end{tabular}
\end{center}
\caption{\label{fig:p7}
Multiphoton simulations with $\Phi=7$ photoisomerisations
for three different arrangements of activation sites. \,
(a) only one disc is activated (disc \# 400) by 7 photons
({\it green});  
(b) seven adjacent discs are activated (discs \# 397 - 403), 
each by one photon ({\it blue}); 
(c) seven discs are activated, 100 discs apart from each 
other (disc \# 100, 200, 300, 400, 500, 600, 700) ({\it red}).
{\bf (A):} Total response versus time. 
{\bf (B):} z-profiles of local response. }
\end{figure}
 

Indeed, for higher $\Phi$ the effect can be much more pronounced than
seen in Fig. \ref{fig:p7}. For example, as shown in
Fig. \ref{fig:p700}, 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 $\Phi=700$).
Minimal response is produced when all the photons are
on one disc (adaptation takes over reducing the response),
and maximal when they are maximally separated,
with as few as possible on each disc (to minimize adaptation).


%%%%%%%%%%%%%%%%%%%------------ fig:p700        fig 3
\begin{figure}[ht] 
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.47\textwidth]{figures/fig3a}
&\includegraphics[width=0.47\textwidth]{figures/fig3b}
\end{tabular}
\end{center}
\caption{\label{fig:p700}
Multiphoton simulations with $\Phi=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, \dots, 785) (blue);
(c) one photon on each of 700 consecutive discs
(discs \#50--749) (red). 
%
{\bf (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.
{\bf (B):}  z-profiles of local response at time 1000 ms 
for each of the three cases in (A).}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\subsection*{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.\ref{fig:p2},
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 $\sim 150$ discs apart,
as their spreads begin to overlap,
and reduces down to $\sim 90\%$ of maximal when the activated discs are
adjacent (or when both photons are on 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.

%%%%%%%%%%%%%%%%%%%------------           fig 4
\begin{figure}[ht] 
\begin{center}
\resizebox*{0.6\textwidth}{!}
{{\includegraphics{figures/fig4}}}
\caption{\label{fig:p2}
Two-photon simulations ($\Phi=2$) exhibiting nonlinear summation
of SPR's as the distance between activated discs decreases.
Left scale shows peak total response as \% of $\text{J}_{\text{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.}
\end{center}
\end{figure}

\section{Conclusions}\label{S:5}
%
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 $\sim 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.

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

\begin{thebibliography}{00}

\bibitem{Andreucci03a}
D. Andreucci, P. Bisegna, G. Caruso, H.E. Hamm and E. DiBenedetto,
\emph{Mathematical Model of
the Spatio-Temporal Dynamics of Second Messengers in Visual
Transduction}, { Biophysical J.} {\bf 85}: 1358--1376, 2003.
%
\bibitem{Andreucci03b}
D. Andreucci, P. Bisegna, and E DiBenedetto,
\emph{Homogenization and Concentrated Capacity for a Problem in Visual
Transduction}, { Annali di Mat. Pura et Appl.} {\bf 182}(4):
375--407, 2003.
%
\bibitem{Baylor79}
D. A. Baylor, T.D. Lamb and K.-W. Yau,
\emph{The Membrane Current of Single Rod Outer Segments},
{ J. Physiol.}, {\bf 288}: 589--611, 1979.
%
\bibitem{caruso}
G. Caruso, H. Khanal, V. Alexiades, F. Rieke, H.E. Hamm, E. DiBenedetto,
\emph{Mathematical and Numerical Modeling of Spatio-Temporal Signaling
in Rod Phototransduction},
{IEE Proc. Systems Biology}, {\bf 152}(3): 119-137, 2005.
%
\bibitem{GrayKeller99}
M. Gray-Keller, W. Denk, B. Shraim and P.B. Detwiler,
\emph{Longitudinal Spread of Second
Messenger Signals in Isolated Rod Outer Segments of Lizards},
{J. Physiol.}, {\bf 519}: 679--692, 1999.
%
\bibitem{Hamer00}
R.D. Hamer,
\emph{Analysis of Ca${}^{++}$-dependent gain changes in PDE activation in
vertebrate rod phototransduction},
{Molecular Vision} {\bf 6}: 265--286, 2000.
%
\bibitem{Hamer03}
R.D .Hamer, S.C. Nicholas, D. Trachina, P.A. Liebman, and T.D. Lamb,
\emph{Multiple steps of phosphorylation of activated rhodopsin can account
for the reproducibility of vertebrate rod single-photon response},
{J. Gen. Physiol.} {\bf 122}: 419--444, 2003.
%
\bibitem{Khanal03} %AIPproc03
H. Khanal, V. Alexiades, E. DiBenedetto and H. Hamm,
\emph{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.
%
\bibitem{Khanal04}
H. Khanal, V. Alexiades, E. DiBenedetto,
\emph{Response of Dark-adapted Retinal Rod Photoreceptors},
pp. 138-145 in { Dynamic Systems and Applications 4},
editor M. Sambandham, Dynamic Publishers, 2004.
%
\bibitem{Koutalos95}
Y Koutalos, K Nakatani, and K-W Yau,
\emph{Cyclic GMP Diffusion Coefficients in Rod Photoreceptors Outer
Segments}, {Biophysical J.}, {\bf 68}: 373--382, 1995.
%
\bibitem{LP92}
T.D. Lamb and E.N. Pugh Jr.,
\emph{A quantitative account of the activation steps involved in
phototransduction in amphibian photoreceptors},
{J. Physiol.} {\bf 449}: 719--758, 1992.
%
\bibitem{Leskov00}
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,
\emph{The Gain of Rod Phototransduction: Reconciliation of Biochemical and
Electrophysical Measurements},
{Neuron}, {\bf 27}: 525--537, 2000.
%
\bibitem{NEP98}
S. Nikonov, N. Engheta and E.N. Pugh, Jr
\emph{Kinetics of Recovery of the Dark-adapted Salamander Rod
Photoresponse}, {J. Gen. Physiol}, {\bf 111}: 7--37, 1998.
%
\bibitem{NLP}
S .Nikonov, T.D. Lamb  and E.N. Pugh, Jr.
\emph{The Role of Steady Phosphodiestarase Activity in the Kinetics
and Sensitivity of the Light-Adapted Salamander Rod Photoresponse},
{J. Gen. Physiol.}, {\bf 116}: 795--824, 2000.
%
\bibitem{PL}
E.N. Pugh, Jr. and T.D. Lamb,
\emph{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.

\end{thebibliography}

\end{document}