This Article Abstract Full Text (PDF) Appendix A & B Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Andreucci, D. Articles by DiBenedetto, E. Search for Related Content PubMed PubMed Citation Articles by Andreucci, D. Articles by DiBenedetto, E.

Mathematical Model of the Spatio-Temporal Dynamics of Second Messengers in Visual Transduction

D. Andreucci ^*, P. Bisegna , G. Caruso , H.E. Hamm  and E. DiBenedetto ^¶

^* Dipartimento di Metodi e Modelli Matematici, Università di Roma La Sapienza, 00161 Rome, Italy; Dipartimento di Ingegneria Civile, Università di Roma Tor Vergata, 00133 Rome, Italy; ITC-CNR, Rome, Italy; Department of Pharmacology, Medical Center, Vanderbilt University, Nashville, Tennessee 37232; and ^¶ Biomathematics Study Group, Department of Mathematics, Stevenson Center, Vanderbilt University, Nashville, Tennessee 37240

Correspondence: Address reprint requests to E. DiBenedetto, E-mail: em.diben{at}Vanderbilt.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION BACKGROUND THEORY NUMERICAL SIMULATIONS DISCUSSION CONCLUSIONS SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
A model describing the role of transversal and longitudinal diffusion of cGMP and Ca²⁺ in signaling in the rod outer segment of vertebrates is developed. Utilizing a novel notion of surface-volume reaction and the mathematical theories of homogenization and concentrated capacity, the diffusion of cGMP and Ca²⁺ in the interdiscal spaces is shown to be reducible to a one-parameter family of diffusion processes taking place on a single rod cross section; whereas the diffusion in the outer shell is shown to be reducible to a diffusion on a cylindrical surface. Moreover, the exterior flux of the former serves as a source term for the latter, alleviating the assumption of a well-stirred cytosol. A previous model of visual transduction that assumes a well-stirred rod outer segment cytosol (and thus contains no spatial information) can be recovered from this model by imposing a "bulk" assumption. The model shows that upon activation of a single rhodopsin, cGMP changes are local, and exhibit both a longitudinal and a transversal component. Consequently, membrane current is also highly localized. The spatial spread of the single photon response along the longitudinal axis of the outer segment is predicted to be 3–5 µm, consistent with experimental data. This approach represents a tool to analyze pointwise signaling dynamics without requiring averaging over the entire cell by global Michaelis-Menten kinetics.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION BACKGROUND THEORY NUMERICAL SIMULATIONS DISCUSSION CONCLUSIONS SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Diffusion of the second messengers cGMP (cyclic-guanosine monophosphate) and Ca²⁺ mediates phototransduction in rod outer segments. This process occurs within the thin layers between the membranous discs (transversal diffusion) and the equally thin outer shell along the plasma membrane (longitudinal diffusion). An open issue is to understand the physical role of the transversal diffusion (Dumke et al., 1994; Lamb et al., 1981; Olson and Pugh, 1993) and of the longitudinal diffusion (Gray-Keller et al., 1999; Lamb et al., 1981; Olson and Pugh, 1993). More importantly, no description exists in the literature of how these two mutually perpendicular diffusions communicate and interact.

First, the physics of the diffusion of the second messengers is in need of a deeper understanding; for example the very same notion of "transversal" and "longitudinal" diffusions are not well defined. Diffusion within the interdiscal spaces is important because these are the only physical spaces through which cGMP can be depleted by phosphodiesterase (PDE) localized on the faces of the discs. Diffusion along the outer shell is equally important because there is a spread of depletion of cGMP to a distance of several discs (Gray-Keller et al., 1999; Matthews, 1986) and because this is the region where the channels and transporters reside and Ca²⁺ enters.

Second, to our knowledge, no pointwise, predictive model exists for the cascade. Such a pointwise model would have to describe the current of the cGMP-gated channels as a function of position on the lateral boundary of the rod in terms of the number of photons hitting the rod. Because the current and [cGMP] are directly linked, this would require the calculation of [cGMP] as a function of space and time. Thus the biophysical issue of understanding the diffusion process and the issue of creating a pointwise model are intertwined.

Researchers have long recognized the importance of diffusion of the second messenger in the cytosol and its space-time dependence (for a recent discussion see Leskov et al. (2000), and in particular Fig. 6 A). Some have generated a mathematical model for the radial diffusion within a single interdiscal space, neglecting the longitudinal diffusion (Dumke et al., 1994). Others have taken into account only the diffusion along the axis of the rod, neglecting space variables in the interdiscal space (Gray-Keller et al., 1999). Although others have attempted to account for both (Lamb et al., 1981; Olson and Pugh, 1993), they do not provide a description of the mechanism by which interdiscal and outer shell diffusions interact.

View larger version (7K): [in this window] [in a new window]   FIGURE 6  Activation by a point source and the diffusion process (Eq. 24). Activation site (point source) at and (A) Plots of (B) Plots of 1 - Jrel(z, 0, t); spr(0, 0.6 s) = 2.86 µm; spr(0, 1.2 s) = 4.17 µm. (C) Plots of values of

	BACKGROUND

TOP ABSTRACT INTRODUCTION BACKGROUND THEORY NUMERICAL SIMULATIONS DISCUSSION CONCLUSIONS SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
The phototransduction cascade
The outer segment of a rod photoreceptor in vertebrates is a right cylinder of height H and radius R + _o, housing a longitudinal stack of n_o equispaced parallel cylinders, called discs, each of radius R, and width _o, and mutually separated by a distance _o (Figs. 1 and 2). Each disc is made up of two functionally independent layers of lipidic membrane where proteins are embedded, such as rhodopsin (Rh), the light receptor, G-protein (G), also called transducin, and cGMP phosphodiesterase (PDE), the effector. These membrane-associated proteins can diffuse on the face of the disc where they are located. The plasma membrane of the rod contains cGMP-gated channels. In absence of light, these channels are open and allow a positive influx of sodium and calcium (Ca²⁺) ions. The space within the rod, and not occupied by the discs, is filled with fluid cytosol, in which cGMP and Ca²⁺ diffuse.

View larger version (20K): [in this window] [in a new window]   FIGURE 1  Geometrical description of the ROS and its outer shell.

View larger version (17K): [in this window] [in a new window]   FIGURE 2  Geometry of the ROS and its discs.

Active PDE^* hydrolyzes cGMP in the cytosol, thereby lowering its concentration. The decrease of concentration of the cGMP causes closure of some of the cGMP-gated channels of the plasma membrane, resulting in a lowering of the influx of positive ions, and thus a lowering of the local current across the outer membrane. Because of the Na⁺/K⁺/Ca²⁺ exchanger that continues to remove Ca²⁺ from the cytosol, there is a decrease in the calcium concentration, which in turn results in an increase in cGMP production by stimulation of Ca²⁺-inhibited guanylyl cyclase, and thus a consequent reopening of the channels. The same decreased Ca²⁺ closes the cycle by causing disactivation of rhodopsin through stimulation of rhodopsin kinase. Rhodopsin ceases activating new G-protein. Thus PDE^* decays to basal, ending depletion of cGMP.

This cascade is well known and it is supported by a sizable amount of published experimental data (e.g., Wald, 1968; Liebman et al., 1987; Schnapf and Baylor, 1987; Pugh and Lamb, 2000; Burns and Baylor, 2001). Its formal mathematical description, however, is less developed. Here we have used a novel mathematical formulation of the excitation phase, which allows us to take into account the complex geometry of the rod outer segment and the diffusion of both second messengers through it.

The geometry of the rod outer segment
Let denote a parameter in the range (0, _o] and let n be a positive integer larger or equal to n_o. Denote by a disc, and of radius R + , where R, , and H are positive numbers. Let and be the cylinders, of height H and cross sections the disc and respectively. The cylinder is included in , is coaxial with it, and it is formally obtained from by setting = 0. The outer shell S is the gap between these cylinders. Coordinates and are introduced as in Fig. 1.

The cylinder houses a longitudinal stack of n parallel, equispaced cylinders C_j, j = 1, 2, ..., n, coaxial with and with cross section a disc of radius R. They are thin in the sense that their height

The C_j are equally spaced, i.e., the upper face of C_j has distance from the lower face of C_j+1, where is a positive number. The first C₁ has distance from the lower face of the rod and the last C_n has distance from the upper face of the rod. The indicated geometry implies that,

Also the volume fraction of the union of the C_j with respect to is

(1)

The upper and lower faces of the cylinders C_j are denoted by whereas L_j denotes their lateral surface (Fig. 2). The spaces between two contiguous cylinders C_j and C_j+1 and within are the interdiscal spaces. We label them by I_j, j = 0, 1, 2, ..., n by defining I_o as the space between the lower face {z = 0} of and the lower face of C₁ and I_n as the space between the upper face {z = H} of and the upper face of C_n. The disc hit by the photon on one of its faces is called the activated disc and is denoted by for some 1 j_o n. For definiteness we assume that the photon hits on its lower face The interdiscal space adjacent to the lower face of is called the activated interdiscal space.

This geometrical description is done in terms of the parameter ranging in the interval (0, _o]. The actual physical width of the discs C_j is _o and their actual number is n_o. Thus the geometrical description of the actual physical rod outer segment is obtained from this by taking for = _o and for n = n_o. Such a description is motivated by the idea of regarding _o as a parameter that will be let go to zero.

	THEORY

TOP ABSTRACT INTRODUCTION BACKGROUND THEORY NUMERICAL SIMULATIONS DISCUSSION CONCLUSIONS SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Diffusion of cGMP and Ca²⁺ within the cytosol
Ca²⁺ and cGMP diffuse in the cytosol, within the rod outer segment (ROS). This is the domain obtained from the cylinder from which the internal discs C_j have been removed. In Fig. 2 it corresponds to the white area that is left in the cylinder when the discs are removed. This domain consists of the outer shell which is a thin cylindrical layer, and the union of the parallel, transversal thin layers of the interdiscal spaces. We denote it by . Because within there are no volume sources for either cGMP or Ca²⁺,

(2)

where t is time, ² is the Laplacian operator in the space variables and D_cG, D_Ca are the respective free diffusivity constants in the cytosol.

The system (Eq. 2) is complemented with source terms supported on the lateral boundary of the rod and the faces of the discs C_j. These sources, positive or negative, are due to volume-to-surface first-order reactions. Postponing the discussion of these sources, here we indicate how we intend to interpret the diffusion phenomenon in the cytosol. One might study the system (Eq. 2) in the structured layered geometry of by regarding such a domain as macroscopic. However consists of layers, transversal or cylindrical, whose thickness is of three orders of magnitude less with respect to the dimensions of the rod outer segment. In situations such as these, the homogenization theory seeks to extract physical information from the system by letting the thickness of the layers go to zero, without altering the total relative volume available for diffusion. In practical terms the parameter _o is replaced with a parameter (0, _o]. Such a parameter is then let to go to zero. Because is the thickness of the new fictitious discs, their number n must increase to infinity to keep constant the volume available for diffusion. From Eq. 1 and the geometry of the rod outer segment, the limit is carried out so that,

(3)

In such a limiting process we impose that although the activated disc changes its width as 0, its lower face, which is the face where the photon hits, remains fixed. Diffusion starts at time t = 0 from a steady state, i.e., with [cGMP](0) and [Ca²⁺](0) constant in the space variables and given by their dark values [cGMP]_dark and [Ca²⁺]_dark.

Boundary source terms for [cGMP]
We will model the activation phase of an idealized experiment by which a single photon hits a disc on the lower face, at coordinate z_o along the axis of the rod.

Production or depletion of molecules of cGMP occurs through binding phenomena on the lower and upper faces of each of the cylinders C_j. Precisely, cGMP is depleted as it binds to dark-activated phosphodiesterase, at a rate,

Here [PDE] is defined as the surface concentration of PDE, uniformly distributed on the total area of the faces of the discs C_j. Precisely denoting by n_o the number of discs and by N_AV Avogadro's number,

where [PDE] is the volumic concentration of PDE regarded as uniformly distributed in the rod. This is the quantity actually being experimentally measured under the assumption of well stirred. Thus,

(4)

Guanylyl cyclase (GC), which is bound to the faces of the discs C_j, synthesizes cGMP. Molecules of guanosine triphosphate bind guanylyl cyclase to generate cGMP. Such activity is modulated by Ca²⁺, which is bound to guanylyl cyclase-activating protein (GCAP). As the concentration of Ca²⁺ decreases, GCAP is released and is free to bind to guanylyl cyclase and to activate it. Diffusion of GCAP is assumed to be negligible, so that molecules of GCAP are essentially still within . Thus only those near the faces of the cylinders C_j and in contact with the GC affect the process. The rate of conversion of guanosine triphosphate into cGMP in terms of [Ca²⁺] is given by an experimental Hill-type relation,

where m is a positive parameter, k_GC is the catalytic rate of guanylyl cyclase, [GC] is the surface density of GC, uniformly distributed on the total area of the faces of the discs C_j, and ß is the Ca²⁺ concentration that achieves half of the maximum rate. Proceeding as before [GC] = _o[GC]/2, where [GC] is the measured volumic concentration of GC regarded as uniformly distributed in the rod. Setting = k_GC[GC], such a rate of production takes the form,

(5)

where is the maximum rate of production of cGMP, corresponding to absence of Ca²⁺. Let be the disc hit by the photon on one of its faces, say for example the lower one, and let be the resulting, pointwise surface density of activated phosphodiesterase. Then assuming full activation of PDE,

(6)

where k^* is the catalytic rate of the light-activated PDE.

Combining these various contributions (Eqs. 4–6), the physical flux of [cGMP] on the faces of the discs C_j, takes the form

(7)

where

[cGMP] does not penetrate the lateral part L_j of the boundary of the cylinders C_j, and does not outflow the boundary of the rod. Therefore [cGMP] has zero physical flux on each of the L_j, on the lateral boundary and on the top (z = H) and bottom (z = 0) of the rod outer segment.

Boundary source terms for [Ca²⁺]
Calcium does not penetrate the discs C_j carrying the rhodopsin, so that its flux across the boundary of each C_j is zero. Calcium does not outflow the bottom (z = 0) and the top (z = H) of the outer segment.

Calcium ions are lost through the lateral boundary of the rod, by electrogenic exchange and are gained by their influx, through the cGMP-activated channels.

The pointwise current density J_ex across the boundary of the rod (charge flux), due to electrogenic exchange, is modeled by the Michaelis-Menten type relation,

(8)

Here is the surface area of the lateral boundary of the rod, j_ex;sat is the maximal, or saturation current as [Ca²⁺] , and K_ex is the half-maximal constant. The pointwise current density J_cG carried by the cGMP-activated channels, across the boundary of the rod, is given by the Hill's type law,

(9)

where j_max is the maximal current as [cGMP] , K_cGMP is the half-maximal constant and is a positive parameter. Let n denote the unit normal to the lateral boundary of the rod, pointing outside the rod. Then, the total pointwise flux of Ca²⁺ across such a surface is given by,

(10)

where is a positive parameter, and f_Ca is a dimensionless number in (0, 1) (see Table 1 footnote). In this formula the product f_CaJ_cG is the portion of the flux of current J_cG carried by Ca²⁺. The fluxes in Eqs. 4–10 contain no recovery mechanisms; accordingly the model is currently valid for the activation phase only.

The diffusivity and capacity coefficients in Eq. 2, in the outer shell and in the activated interdiscal space are, roughly speaking, multiplied by _o/ to compensate for a shrinkage of the domain of the same order. This is the role of concentrated capacity (Andreucci, 1990; Ciarlet and Lods, 1996; Magenes, 1998). The precise mathematical implementation of this idea is in Andreucci et al. (2002, 2003) (see also Supplementary Material, Appendix A, A2 and A3).

Here we discuss what the limiting [cGMP] and [Ca²⁺] look like and the equations they satisfy.

As 0 the layered domain tends formally to the cylinder and the activated interdiscal spaces tend to the disc Likewise the outer shell S tends to the surface The functions [cGMP](x, t) and [Ca²⁺](x, t) generate three pairs of limiting functions, each representing [cGMP] and [Ca²⁺] in different parts of the rod outer segment. Precisely:

The interior limits [cGMP] and [Ca²⁺] are defined on a volumic domain and their physical dimensions remain unchanged. Although the last two limits are defined on surfaces, they keep their physical dimensions in µM. To make this point precise, consider for example the limit [cGMP]_s in the limiting outer shell S. Describe the approximating outer shell S in terms of cylindrical coordinates,

and express [cGMP] in terms of these coordinates. It can be shown that [cGMP]_s is a function of (, z, t) defined as the limit, as 0, of the radial integral average of [cGMP] in the approximating outer shell S, i.e., (see Supplementary Material, Appendix A, A3),

This formula implies that [cGMP]_s, although defined on the surface (0, 2] x (0, H), keeps its physical dimensions in µM, because it is the limit of integral averages of volume densities. The factor 1/ in this limiting formula, arises form the rescaling of the capacity and diffusivity coefficients in the outer shell, as part of the process of concentrated capacity. Similar considerations apply to (see Supplementary Material, Appendix A A3).

We next give the equations satisfied by these limiting quantities, each in its own geometric portion of the rod outer segment. More importantly we elucidate how these seemingly separate diffusion processes interact with each other. To gain in simplicity we do this for the equations satisfied by the limiting [cGMP]. The analogs for [Ca²⁺] are in Supplementary Material, Appendix A, A3, where justifications and proofs are provided.

The limiting equations will contain in various forms the forcing terms generated by Eq. 7. To simplify the symbolism we will set,

(11)

In these expressions [cGMP] and [Ca²⁺] are the interior limits of [cGMP] and [Ca²⁺] and [cGMP]_o and [Ca²⁺]_o are the limits at the activated level z_o.

Form of the interior limit of [cGMP]
The interior homogenized limit is computed by a local average in each of the interdiscal spaces I_j. Such an average takes into account the boundary conditions in Eq. 7. The net result is that the interior limiting [cGMP] satisfies the equation,

(12)

Here is the Laplace diffusion operator acting only on the transversal variables i.e., formally,

Because is a function of the transversal variables and the longitudinal variable z, these can be regarded as diffusion processes, parameterized with z (0, H), taking place on the disc Thus the volumic diffusion in Eq. 2 in the layered structure of the rod, is transformed into a family of two-dimensional diffusions. Also, the homogenized limit transforms the boundary fluxes in Eq. 7 into volumic source terms holding in .

Form of the limiting [cGMP]_o at the special level z_o
The limiting [cGMP]_o on the activated level z_o is also computed by averaging Eq. 2 over the interdiscal space adjacent to the activated disc and by letting its thickness go to zero. The limiting [cGMP]_o satisfies the equation (Andreucci et al., 2003; also see Supplementary Material, Appendix A, A3),

(13)

on the activated limiting disc _R x {z_o}. Thus, also at the activated level z_o, the volumic diffusion in Eq. 2 is transformed into a two-dimensional diffusion on the layer _R x {z_o} and the fluxes in Eq. 7 are transformed into sources defined in the interior of the same disc and keeping the same form. Note that in this case the limit equation contains also the term F^* due to activated PDE.

Form of the limiting [cGMP]_s in the outer shell
The limiting [cGMP]_s(, z, t) on S is a function of the angular variable [0, 2), of the longitudinal variable z (0, H) and of time. Outside the activated level z_o it must equal the interior limit when this is computed on S. For consistency, on the activated level z_o it must equal the limiting [cGMP]_o when this is computed on S. Therefore,

(14)

Moreover, the interior limit [cGMP] and the limit [cGMP]_o on the activated level z_o are linked to the limit [cGMP]_s in a more essential way. Describe the limiting cylinder in cylindrical coordinates (, , z). Then the fluxes of [cGMP] and [cGMP]_o on S are given by

Denote by the Laplace-Beltrami diffusion operator on S, i.e., formally

Then these fluxes and the limiting [cGMP]_s(, z, t) on the outer shell, satisfy the surface-diffusion equation,

(15)

in S. Here _zo is the Dirac delta function on S with mass on the level z_o. Thus, the diffusion of [cGMP]_s on the limiting outer shell S, is forced by the exterior fluxes on S, of the interior limit [cGMP] and the limit [cGMP]_o on the activated z_o.

This is the biophysical law by which the homogenized-concentrated limiting diffusions interact with each other. Although it is somewhat intuitive that [cGMP] coming from the transversal interstices should provide the driving force for the movement of the [cGMP] on the longitudinal surface S, Eqs. 14 and 15 provide a precise law by which this occurs. In particular they contain a precise combination of the original geometric parameters _o, _o, , . This combination of geometric parameters expresses the balance of mass between [cGMP]_s on S and the outflow through S of the interior [cGMP]. Multiplying Eq. 15 by _o, the left-hand side represents the pointwise space-time variation of _o[cGMP]_s. The latter quantity can be regarded as a surface density of cGMP concentrated on S, starting from the original shell S_o of thickness _o. The factor (1 - _o) on the right-hand side signifies that only a fraction of (1 - _o) of S is exposed to the outflow of the homogenized interior limit of [cGMP]. This is the same fraction of surface exposed to inflow/outflow of cGMP from the interdiscal spaces into the outer shell S_o in the original, nonhomogenized configuration of the rod outer segment.

An integral version of Eqs. 11–15
The form of Eqs. 14 and 15 precisely describes the interaction between the "interior" and the "boundary" diffusion of the second messengers cGMP and Ca²⁺. However, from a mathematical point of view, Eqs. 11–15 must be interpreted in a suitable weak sense (Andreucci et al., 2003 and Suppplementary Material, Appendix A, A4). Here we give an integrated form of Eqs. 11–15, which is a particular case of such a weak formulation.

Far from being an artificial construct, such a rigorous mathematical interpretation has dense physical consequences that will be brought to light in the next section. Its main feature is that it combines the geometrical properties of the various compartments. This permits one to specialize it under various simplifying assumptions such as transverse or global well-stirred cytosol.

One of the outcomes of such an integral form is that Eqs. 11–15 contain, as a particular case, the known well-stirred theories (either in the transversal variables (x, y) or in all space variables). In addition, even in the well-stirred assumption, they represent a significant improvement with respect to the existing theories in that they distinguish the diffusion of the second messengers outside the activation site z_o from the diffusion on the activated level z_o. This is the content of the next mathematical derivations from Eqs. 11–15.

Integrate Eq. 12 over the disc _R x {z} at the generic level z. Applying the Gauss-Green theorem,

where d is the line measure on the circle An entirely similar operation on Eq. 13 yields,

Next integrate Eq. 15 on _R x {z} for a generic level z (0, H), to obtain

In performing such an integration we have taken into account that the angular part of the Laplace-Beltrami operator does not give any contribution because [cGMP](, z, t) is a periodic function of . The last two integrals are substituted from the previous two formulae, thereby eliminating the explicit calculation of the flux of [cGMP] across S. Regrouping the resulting terms we arrive at

(16)

This is a particular case of the notion of "weak formulation" for the problem (Eqs. 11–15). A more general weak formulation is in Andreucci et al. (2002, 2003) (see also Supplementary Material, Appendix A, A4).

Cytosol well stirred in the transversal variables (x, y)
Assume the cytosol is well stirred in the transversal variables (x, y). Thus, the rod outer segment is ideally lumped on its axis and transversal diffusion effects are immaterial. Such an assumption is suggested by the idea that the system diffuses with infinite speed on each transversal cross section and thereby responds with an instantaneous transversal equilibration. Although not rigorous on physical and mathematical grounds, such an assumption here is made in the sense that errors originating from it are neglected.

The analysis below will permit one to compare our model to the existing ones based on the assumption of well stirred.

If [cGMP] and [Ca²⁺] are regarded as lumped on the axis of the rod, they depend only on z and t, and are independent of (x, y). Because there is no dependence on the (x, y) variables, by Eq. 14

We insert this information into the formula (Eq. 16) and compute the resulting integrals, to obtain,

Set,

These two parameters have a geometric and physical significance. Specifically, up to higher-order corrections, f_A is the fraction of the cross-sectional area of the outer segment that is available for longitudinal diffusion, and f_V is the fraction of the total outer segment volume occupied by the cytosol (Lamb et al., 1981; Olson and Pugh, 1993). Then, dividing by f_V the previous equation takes the more concise form,

(17)

This is the law of diffusion of [cGMP] under the assumption that the cytosol is well stirred in the transversal variables. A key feature is that it distinguishes between diffusion outside the activated level z_o and diffusion at z_o by the action of the Dirac delta function If z is different than the activated level z_o, Eq. 17 implies,

(18)

Equation 18 is formally similar to a model proposed by Gray-Keller et al. (1999). In that work, however, the term F_* due to activation is distributed along the longitudinal variable z. To elucidate the effect of the activation site z_o, one has to compute Eq. 17 for z = z_o. In view of the Dirac delta function computation of Eq. 17 for z = z_o can be done only in the sense of distributions (DiBenedetto, 2002, Chap. VII). For example, we might integrate in dz over a small interval (z_o - a, z_o + a) about z_o, where Letting a 0, we obtain a relation expressing the conservation of mass of [cGMP] across the activated level z_o.

Globally well-stirred cytosol
Regard now the rod as a homogeneous bag of cytosol, and [cGMP], [Ca²⁺], [PDE], [PDE^*], as lumped quantities depending only on time. Thus, in particular [cGMP] = [cGMP]_o = [cGMP]_s and similarly for [Ca²⁺]. Rewrite Eq. 17 in the integrated form,

(17)'

Then we may set to zero the term involving the z-derivative and compute the remaining integrals to get,

Now f_V is of the order of one and _o is of three orders of magnitude smaller than H and R. Therefore,

From the expression of f_V and the form (Eq. 11) of F_*,

where V_cyto is the volume of the outer rod segment available for diffusion. Therefore, the assumption of well stirred in all the space variables yields, up to higher-order corrections, the dynamic equation,

(19)

where F is defined in the first of Eq. 11. A similar analysis for calcium gives,

(20)

where the various parameters are the ones occurring in the flux condition (Eqs. 8–10) and discussed there. These formulae coincide with Nikonov et al. (2000) (A3 and A4; p. 39), upon identifying the various parameters. This is a validating point of our model as it fits the experimental data at least as well as Nikonov's model does.

Flexibility of the model
In Eqs. 12–19 and throughout the development of the theory, the functions F, F_o, F_* are those defined in Eq. 11. Although this has been done for notational simplicity, nowhere in the arguments does the specific form in Eq. 11 enter, in the calculation of the homogenized-concentrated limit (Eqs. 12–19). Thus, such a limit is independent of the form (Eq. 11) of F, F_o, F_*, provided these are bounded smooth functions of [cGMP] and [Ca²⁺].

These functions originated from modeling the production and depletion mechanisms of [cGMP] on the faces of the discs. Variants or refinements of these mechanisms might be incorporated into these functions and would produce the very same homogenized-concentrated limit (Eqs. 12–19) with the newly redefined forms of F, F_o, F_*. This affords considerable flexibility to the model.

As an example, consider the rate of production of [cGMP] due to membrane-bound guanylyl cyclase GC, leading to Eq. 5. The mechanism we have adopted is that proposed by Forti et al. (1989) and Gray-Keller et al. (1999). A refinement of such a mechanism is in Nikonov et al. (2000), although in a volumic well-stirred form. When interpreted as a boundary flux it reads,

where m_cyc is a Hill's exponent, K_cyc is a positive parameter, and [GC] is the surface density of GC, regarded as uniformly distributed on the total area of the faces of the discs C_j. Because GC is inhibited by Ca²⁺ its maximum catalytic rate k_GC,max occurs for [Ca²⁺] 0 and its minimum catalytic rate k_GC,max occurs theoretically as [Ca²⁺] . Surface bound GC is converted into volumic [GC], by the same surface-volume mechanism leading to Eq. 5. Setting,

the cyclase mediated rate of production of cGMP takes the form,

This contributes to the total flux of [cGMP] in the faces of the discs. The factor accounts for the surface-volume interpretation. A similar form of the rate of production of cGMP due to cyclase is in Nikonov et al. (2000), where, however, such a term is volumic and offers no spatial resolution. In that work, the value of _max is reported as 50 µM s^-1 whereas (_min/_max) = 0.02. Thus _min is about two orders of magnitude smaller than _max. If _min is neglected by setting it to be zero, the previous rate of production reduces to Eq. 5 with _max = and K_cyc = ß. Starting with this new rate of production, the function F and F_o in Eq. 11 are modified into

(21)

These reduce to Eq. 11 if _min = 0. The function F_* in Eq. 11 remains unchanged.

Multiple photon activation
To gain in simplicity, the theory has been developed to model a single-photon response. Multiple-photon responses are easily treated as follows. Assume discs are activated at the levels z_o,1, z_o,2, ..., z_o,N. The very same theory applies and the resulting equations remain the same in nature with the following variants. The function remains defined as in Eq. 11. Each of the activated discs now generates its own forcing terms F_o, and F_*, defined as in Eq. 11 each at the level z_o,. Precisely,

In these expressions [cGMP]_o, and [Ca²⁺]_o, are the limits of the approximating [cGMP] and [Ca²⁺] on the activated level z_o,. Also [PDE^*]_, is the surface density (in µmol/µm²) of [PDE^*] on the disc _R x {z_o,}.

Equation 12 remains unaltered. Equation 13 is replaced by N equations of exactly the same form as Eq. 13, each representing the limiting [cGMP] on its own activated disc. Precisely,

(13)'

each on its own activated limiting disc _R x {z_o,}. All the terms of Eq. 15 remain the same except the last, which now has to account for the presence of N activated "levels" and takes the form,

24

(15)'

in S. Here _zo, are the Dirac delta functions on S with masses on the levels z_o,.

If several photoisomerizations occur on the same disc at level z_o, this is accounted for in the form of the function F_*, defined in Eq. 11' for its own level z_o,.

If the rate of generation of [cGMP] due to cyclase is modeled by Eq. 21, and several isomerizations occur on more than one disc, then Eq. 11' is modified accordingly. Thus, the model encompasses a large spectrum of experimental settings.

	NUMERICAL SIMULATIONS

TOP ABSTRACT INTRODUCTION BACKGROUND THEORY NUMERICAL SIMULATIONS DISCUSSION CONCLUSIONS SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
The strength of the homogenized model (Eqs. 11–15) is in its spatio-temporal resolution. The numerical simulations presented below build on this feature. The main results are a suppression of total dark current, in close agreement with the experimental data (0.5–1.5%; F. Rieke, unpublished) and the phenomenon of local spread of excitation along the axis of the rod outer segment (Baylor et al., 1979 a,b; Lamb et al., 1981; Gray-Keller et al., 1999; Matthews, 1986). Let J_cG and J_ex be defined as in Eqs. 8 and 9 and set

(22a)

As z ranges over (0, H) and ranges over [0, 2), the variables (z, ) range over the lateral surface S of the rod outer segment. When computed at t = 0 both [cGMP] and [Ca²⁺] are constant and corresponding to their dark values. Consequently J_dark is also a constant. The plots below show relative currents and their deviation from the dark state, i.e.,

(22b)

where dS is the surface measure on S. A set of simulation parameters for the salamander have been collected from a large cross section of the literature and numerically tested for consistency in Khanal et al. (2003) and are reproduced in the Table 1 parameters.

Numerical simulations take also into account the equations for Ca²⁺ (in their weak form, Andreucci et al., 2002, 2003; see also Supplementary Material, Appendix A, A4 Although simulations could be done for a number of modeling combinations, as indicated previously, we assume at this point that a single rhodopsin is activated on a disc at level z_o. In Eqs. 11–15, we take the forcing terms F and F_o as the F^new and defined in Eq. 21. This way all the terms in Eqs. 11– 15 are well identified except the form of the function [PDE^*] on the activated level z_o, which enters in the forcing term F_* of Eq. 11. In the simulations, such a function has been taken two different ways. We call the first "diffused activation" and the second "pointwise activation" (see below).

Diffused activation on a single disc by a lumped/bulk model
The function [PDE^*] depends only on time, and if activation occurs on a single disc, [PDE^*] = E^*/N_AVR², where E^* is the number of activated effectors. Initially E^*(0) = 0. Consider the case when a single rhodopsin is activated, and remains active along the numerical simulation. Then following Nikonov et al. (2000), (see formula (A2) in this reference), activation of the effector starts and continues (at least for the activation phase) by the differential equation,

(23a)

whose solution is

(23b)

Here _RE is the rate of activation of the effector for a single rhodopsin and k_E is the rate of inactivation of E^*. The values of these parameters are reported in the previous table. Such an activation mechanism as proposed in Nikonov et al. (2000), assumes a well-stirred environment. We are assuming that PDE^* is uniformly distributed on the activated disc at level z_o. However, only the input [PDE^*], on the activated disc _R x {z_o}, is taken to be well stirred in the transversal variable of the disc. Starting from this input, the evolution process involves all the spatio-temporal variables and delivers pointwise information on [cGMP], [Ca²⁺] and the resulting current. Starting from t = 0 the number E^* of effectors grows as in Eq. 23b but remain confined on a fixed disc, which, in the simulation, is taken at the middle level of the rod outer segment. Depletion of [cGMP] occurs on the rod as a function of position and time. Current is generated on the boundary of the rod as a function of position and time. Because the process is radially symmetric the current depends only on the variable z along the longitudinal axis of the rod. All simulations are run for 1.2 s with a time-step integration of 10 ms.

As a way of comparing the space-resolved, homogenized model with existing well-stirred ones, we have also generated numerical simulations for: a), the lumped/bulk model as arising in Eqs. 19 and 20; b), the model well stirred in the transversal variables (x, y) as appearing in Eq. 17 with the companion equation for calcium. Activation occurs at the level z_o.

In all cases, the PDE activation mechanism is the one described in Eq. 23b. A result of the numerical simulations is that current suppression is less the more space resolved the model is (Fig. 3). This is due to the damping effect of the diffusion mechanism. In a well-stirred model all the molecules of cGMP in the rod are regarded as contributing instantaneously to the closing of the channels, thereby generating a larger current suppression.

View larger version (22K): [in this window] [in a new window]   FIGURE 3  Activation by a single R* with [PDE*] given by the lumped/bulk model (Eqs. 23a and 23b). History of the relative integrated current Jint. After 800 ms, dark current suppression is 1.54% for a well-stirred ROS, 1.34% for ROS well stirred in the transversal variables, and 0.94% for a fully space-resolved model. Thus, the more space resolved is the model, the less is the current suppression. The three lowest curves represent current suppression for the homogenized model activated by a punctual source as in Eq. 24. The three activation sites are in Fig. 4. Current suppression is dramatically dependent on the activation site.

View larger version (9K): [in this window] [in a new window]   FIGURE 4  Transversal cross section of the rod outer segment at the level zo of the disc where activation occurs. Activations are simulated at 0, at and and on the boundary of the rod at R and

View larger version (7K): [in this window] [in a new window]   FIGURE 5  Plots of 1 - Jrel(z, t) at times 0.2, 0.4, 0.6, 0.8, 1.2 s. (A, B) Activation by Eqs. 23a and 23b. (A) Model with ROS well stirred in the transversal variables; spr(0.6 s) = 3.28 µm; spr(1.2 s) = 4.78 µm. (B) Homogenized model with disc R x {zo} activated; spr(0.6 s) = 3.21 µm; spr(1.2 s) = 4.58 µm. (C) Activation by the diffusion process (Eq. 24). Homogenized model with activation site at the center of the disc R x {zo}; spr(0.6 s) = 2.66 µm; spr(1.2 s) = 4.01 µm. Common to these panels is that they exhibit radially symmetric solutions and therefore there is no dependence on the angular variable . Both dark current suppression and spread decrease for higher space resolution of the model.

1. E^* is zero at time t = 0 when activation starts.
   
2. Activation of the effector E starts at time t = 0 by a single R^* and it continues for all the duration of the simulation.
   
3. The activating rhodopsin is localized at a fixed point (x_o, y_o) on the disc _R x {z_o}. Therefore its action is that of a Dirac point-mass
   
4. The resulting molecules of E^* diffuse on the disc _R x {z_o} and are depleted by a decay term of the type k_EE^*. Such a depletion term has been inserted to keep the model consistent with Eq. 23a above.
   

Thus, the surface density of molecules P^* = [PDE^*]N_AV is a solution of,

(24)

Integrating this over _R x {z_o} gives precisely Eq. 23a for the variable

In this sense Eq. 24 can be regarded as a space-resolved version of Eq. 23a. By varying the position of (x_o, y_o) on _R x {z_o} one can trace numerically the spatio-temporal dependence of the response.

Simulations have been run for three different activation sites (locations where the photon hits). The first is at the center of the disc; the second is half-way between the center and the rim (in polar coordinates (, ) such a site is and see Fig. 4); the third is exactly at the rim of the disc ( = R and ). In either case the current, as a function of the longitudinal variable z and time t, has been recorded at three angular locations on the boundary of the rod. As indicated in Fig. 4 at (the farthest point on the disc, from the activation site), (the closest point on the surface of the rod to the activation site), and = 0 (at some intermediate distance from activation).

Fig. 5 C reports current suppression for the fully space-time resolved ROS with homogenized model, with activation mechanism given by Eq. 24. The activating point source is placed at the center of the disc _R x {z_o} so that the response is radial. Relative current suppression (1 - J_rel(z, t)) is plotted along the axis of the rod for the times t = 0.2, 0.4, 0.6, 0.8, 1.2 s.

In Figs. 6 and 7, the activation mechanism is that of Eq. 24, for the fully space-time resolved ROS with homogenized model. In Fig. 6 the activating point source is placed half-way between the center and the rim of the disc _R x {z_o}. In each of the panels the relative current suppression (1 - J_rel(z, , t)) is plotted, as a function of z at three different sites on the lateral boundary of the ROS. Precisely at = /2 in panel A; at = 0 in panel B and at = - /2 in panel C. In each of the panels, the responses are plotted at the same times. In Fig. 7 the activating point source is placed exactly at the rim of the disc _R x {z_o}. In each of the panels the relative current suppression (1 - J_rel(z,