Predictive Mathematical Models for the Spread and Treatment of Hyperoxia-induced Photoreceptor Degeneration in Retinitis Pigmentosa

Abstract

Purpose: To determine whether the oxygen toxicity hypothesis can explain the distinctive spatio-temporal patterns of retinal degeneration associated with human retinitis pigmentosa (RP) and to predict the effects of antioxidant and trophic factor treatments under this hypothesis.

Methods: Three mathematical models were derived to describe the evolution of the retinal oxygen concentration and photoreceptor density over time. The first model considers only hyperoxia-induced degeneration, while the second and third models include mutation-induced rod and cone loss respectively. The models were formulated as systems of partial differential equations, defined on a two-dimensional domain spanning the region between the foveal center and the ora serrata, and were solved numerically using the finite element method.

Results: The mathematical models recapitulate patterns of retinal degeneration which involve preferential loss of photoreceptors in the parafoveal/perifoveal and far-peripheral retina, while those which involve a preferential loss of midperipheral photoreceptors cannot be reproduced. Treatment with antioxidants or trophic factors is predicted to delay, halt, or partially reverse retinal degeneration, depending upon the strength and timing of treatment and disease severity.

Conclusions: The model simulations indicate that while the oxygen toxicity hypothesis is sufficient to explain some of the patterns of retinal degeneration observed in human RP, additional mechanisms are necessary to explain the full range of behaviors. The models further suggest that antioxidant and trophic factor treatments have the potential to reduce hyperoxia-induced disease severity and that, where possible, these treatments should be targeted at retinal regions with low photoreceptor density to maximize their efficacy.

The group of retinal diseases known collectively as retinitis pigmentosa (RP) is the world's most common inherited retinal degeneration.¹ Affecting approximately 1 in 4000 individuals worldwide,^2,3 RP causes a progressive loss of visual function, culminating in complete blindness by approximately the age of 60.³ While a number of hypotheses have been proposed as to the disease etiology, at present no treatments are clinically available to halt or reverse disease progression.⁴ In this paper, we use mathematical models to investigate one such hypothesis, namely the oxygen toxicity hypothesis.^5–7 

Following the development of night-blindness, caused by a loss of rod function, RP progresses via the degeneration of photoreceptors in small, disc-shaped patches, located mainly within the retinal midperiphery (see Fig. 1 for a description of the retinal regions).^8–12 These patches then expand and coalesce, leading to a number of distinctive spatio-temporal patterns of retinal degeneration, correlating with visual field loss patterns. Grover et al.¹³ classified visual field loss into three classes (Patterns 1, 2, and 3), which may be subdivided into six subclasses: Patterns 1A, 1B, 2A, 2B, 2C, and 3 (see Fig. 2). Pattern 1A involves a constriction of the peripheral visual field, while 1B also includes a parafoveal/perifoveal ring scotoma (blind spot). Pattern 2 begins with either a loss of nasal or superior nasal vision (2A), out from which an arcuate (bow shaped) scotoma winds inferiorly through the midperiphery to the temporal side, or with a superior temporal (2C) or complete superior (2B) restriction, out from which an arcuate scotoma winds inferiorly through the midperiphery to the nasal side. In each case, only the central and inferior visual field is preserved. Lastly, in Pattern 3, a complete or partial midperipheral ring scotoma develops initially. The scotoma then expands into either the superior or inferior periphery, leaving a central island of vision, together with a U- or n-shaped peripheral visual field, the arms of the ‘U' or ‘n' subsequently retracting. For all patterns, central vision is generally the best preserved, being lost in the final stage of the disease. Similar patterns have been observed in numerous studies.^14–23 

Figure 1

 

Schematic showing the different regions of the retina (in the left eye, as viewed from the front).

Figure 2

 

Patterns of visual field loss as described by Grover et al.¹³ Vision loss can be classified into three cases and six subcases. Large gray arrows indicate transitions between stages of vision loss and small red arrows indicate the direction of scotoma propagation. See text for details.

While the initial loss of photoreceptors may be attributed to genetic mutations,^2,3 the cause of the expansion of degenerate patches remains a matter of speculation. Five main hypotheses have been proposed to explain this phenomenon: the oxygen toxicity,^5–7 rod trophic factor,^24–30 toxic substance,³¹ microglia,³² and metabolic dysregulation (mTOR)^33,34 hypotheses. The oxygen toxicity hypothesis suggests that the mutation-induced loss of photoreceptors causes an increase in the retinal oxygen tension, such that conditions become toxic (hyperoxic) for the remaining photoreceptors. This results in a positive feedback loop, as subsequent photoreceptor loss further increases oxygen levels, which are maintained due to the inability of the choroid (the retina's chief oxygen supply) to autoregulate.^5,35,36 

Hyperoxia-induced photoreceptor degeneration may be combated using a combination of antioxidant and trophic factor treatments. Antioxidants neutralize the excess reactive oxygen species generated under hyperoxia,³⁷ while trophic factors are thought to increase photoreceptor resistance to apoptosis under hyperoxic conditions.^38,39 Retinal antioxidant and trophic factor concentrations can be increased through dietary intake (antioxidants only) and using ocular gene therapy or encapsulated cell technology.^4,40–44 Both treatments have been shown to reduce retinal degeneration in RP.^{35,38,39,43,45–47} 

Several mathematical models have been developed to describe the rod trophic factor,^48–54 mTOR,⁵⁴ and toxic substance^55–57 hypotheses (see Roberts et al.⁵⁸ for a detailed review of the current state-of-the-art in retinal modeling, including models of RP). In an earlier study we developed a one-dimensional (1D) model (fovea to ora serrata) to describe retinal degeneration under the oxygen toxicity hypothesis.⁵⁹ To the best of our knowledge, this was the first model to consider this hypothesis, or to account for the spatial distribution of rods and cones. In this paper, we extend this model to two dimensions (2D), incorporating the azimuthal dimension, so that we can consider scenarios where radial symmetry is broken. Comparing our modeling predictions with the patterns of retinal degeneration that can be inferred from Grover et al.,¹³ we judge the strengths and weaknesses of the oxygen toxicity hypothesis by our models' ability to recapitulate the characteristic patterns associated with RP. We also use our models to predict the effects of treatment with antioxidants and trophic factors under the oxygen toxicity hypothesis. 

Methods

In what follows we provide a succinct description of our three mathematical models. The first model considers only hyperoxia-induced photoreceptor degeneration, while the second and third include mutation-induced rod and cone degeneration respectively. For further details on the assumptions behind the models see Roberts et al.⁵⁹ The equations were solved using the finite element method, using the Portable, Extensible Toolkit for Scientific Computation libraries (available in the public domain, www.mcs.anl.gov/petsc/) together with piecewise linear basis functions (see Roberts⁶⁰ for further details). 

Our models focus on the outer retina, which we define as the region between Bruch's membrane (BM)—where the retina abuts the inner capillary layer of the choroid known as the choriocapillaris (CC)—and the photoreceptor inner segments (ISs) inclusive. This region also contains the retinal pigment epithelium and the photoreceptor outer segments (OSs). RP typically occurs as a rod–cone dystrophy, in which rod photoreceptors are affected earlier and more severely than cone photoreceptors; however, cone–rod varieties also exist.^2,3 Our models account for the heterogeneous distribution of photoreceptors across the retina, distinguishing between rod and cone photoreceptors when considering the effects of mutation-induced rod or cone loss in the rod–cone and cone–rod forms of RP. 

Adopting a spherical polar coordinate system, we model the eye as a spherical cap (see Fig. 3), orientating the coordinate axes such that the center of the fovea lies at 0° eccentricity (θ = 0 [rad], noting that all angles are given in radians). We model the retina in 2D: polar, θ ∈ [0, Θ] (rad), and azimuthal, ϕ ∈ [0, 2π) (rad), neglecting the radial dimension, r (m), because the retinal thickness (80–320 [μm], see Webvision, available in the pubic domain, http://webvision.med.utah.edu/) is two orders of magnitude smaller than its radius of curvature (≈ 1.2 [cm]⁶¹). The ora serrata is assumed to lie at θ = Θ ≈ 1.33 (rad), based on data supplied by Curcio et al.⁶² (noting that, unlike the retina depicted in Fig. 3a, the retinas measured by Curcio et al.⁶² do not extend beyond the equator at θ = π/2 [rad] in the sector of the eye from which these measurements are taken). We assume the portion of the retina we are modeling comprises an effectively depth-averaged section between BM and the outer tips of the photoreceptor ISs. 

Figure 3

 

Diagrams to show the model geometry and 2D visualization. (a) Diagram of the human eye, viewed in the transverse plane. Model equations are solved on a 2D domain that extends from the center of the fovea (θ = 0 [rad]) to the ora serrata (θ = Θ [rad]). Diagram reproduced, with modifications, from http://www.nei.nih.gov/health/coloboma/coloboma, courtesy: National Eye Institute, National Institutes of Health (NEI/NIH). (b) Equations are solved on a spherical cap (the healthy photoreceptor distribution is shown here), ignoring the optic nerve head and projected onto the xy-plane for visualization. This has the effect of preserving area near the center of the retina and compressing the area in the polar direction toward the periphery.

Oxygen diffuses freely through the retina, is consumed by retinal tissue and exchanged with the CC, while photoreceptors maintain their healthy density or regenerate lost biomass under normoxia (healthy oxygen levels) and degenerate under hyperoxic conditions. We define partial differential equations for the retinal oxygen concentration, c(θ,ϕ,t) (mol m^–3), and the photoreceptor density, p(θ,ϕ,t) (photoreceptors m^–2), over time, t (s), as follows  

\(\def\upalpha{\unicode[Times]{x3B1}}\)\(\def\upbeta{\unicode[Times]{x3B2}}\)\(\def\upgamma{\unicode[Times]{x3B3}}\)\(\def\updelta{\unicode[Times]{x3B4}}\)\(\def\upvarepsilon{\unicode[Times]{x3B5}}\)\(\def\upzeta{\unicode[Times]{x3B6}}\)\(\def\upeta{\unicode[Times]{x3B7}}\)\(\def\uptheta{\unicode[Times]{x3B8}}\)\(\def\upiota{\unicode[Times]{x3B9}}\)\(\def\upkappa{\unicode[Times]{x3BA}}\)\(\def\uplambda{\unicode[Times]{x3BB}}\)\(\def\upmu{\unicode[Times]{x3BC}}\)\(\def\upnu{\unicode[Times]{x3BD}}\)\(\def\upxi{\unicode[Times]{x3BE}}\)\(\def\upomicron{\unicode[Times]{x3BF}}\)\(\def\uppi{\unicode[Times]{x3C0}}\)\(\def\uprho{\unicode[Times]{x3C1}}\)\(\def\upsigma{\unicode[Times]{x3C3}}\)\(\def\uptau{\unicode[Times]{x3C4}}\)\(\def\upupsilon{\unicode[Times]{x3C5}}\)\(\def\upphi{\unicode[Times]{x3C6}}\)\(\def\upchi{\unicode[Times]{x3C7}}\)\(\def\uppsy{\unicode[Times]{x3C8}}\)\(\def\upomega{\unicode[Times]{x3C9}}\)\(\def\bialpha{\boldsymbol{\alpha}}\)\(\def\bibeta{\boldsymbol{\beta}}\)\(\def\bigamma{\boldsymbol{\gamma}}\)\(\def\bidelta{\boldsymbol{\delta}}\)\(\def\bivarepsilon{\boldsymbol{\varepsilon}}\)\(\def\bizeta{\boldsymbol{\zeta}}\)\(\def\bieta{\boldsymbol{\eta}}\)\(\def\bitheta{\boldsymbol{\theta}}\)\(\def\biiota{\boldsymbol{\iota}}\)\(\def\bikappa{\boldsymbol{\kappa}}\)\(\def\bilambda{\boldsymbol{\lambda}}\)\(\def\bimu{\boldsymbol{\mu}}\)\(\def\binu{\boldsymbol{\nu}}\)\(\def\bixi{\boldsymbol{\xi}}\)\(\def\biomicron{\boldsymbol{\micron}}\)\(\def\bipi{\boldsymbol{\pi}}\)\(\def\birho{\boldsymbol{\rho}}\)\(\def\bisigma{\boldsymbol{\sigma}}\)\(\def\bitau{\boldsymbol{\tau}}\)\(\def\biupsilon{\boldsymbol{\upsilon}}\)\(\def\biphi{\boldsymbol{\phi}}\)\(\def\bichi{\boldsymbol{\chi}}\)\(\def\bipsy{\boldsymbol{\psy}}\)\(\def\biomega{\boldsymbol{\omega}}\)\(\def\bupalpha{\unicode[Times]{x1D6C2}}\)\(\def\bupbeta{\unicode[Times]{x1D6C3}}\)\(\def\bupgamma{\unicode[Times]{x1D6C4}}\)\(\def\bupdelta{\unicode[Times]{x1D6C5}}\)\(\def\bupepsilon{\unicode[Times]{x1D6C6}}\)\(\def\bupvarepsilon{\unicode[Times]{x1D6DC}}\)\(\def\bupzeta{\unicode[Times]{x1D6C7}}\)\(\def\bupeta{\unicode[Times]{x1D6C8}}\)\(\def\buptheta{\unicode[Times]{x1D6C9}}\)\(\def\bupiota{\unicode[Times]{x1D6CA}}\)\(\def\bupkappa{\unicode[Times]{x1D6CB}}\)\(\def\buplambda{\unicode[Times]{x1D6CC}}\)\(\def\bupmu{\unicode[Times]{x1D6CD}}\)\(\def\bupnu{\unicode[Times]{x1D6CE}}\)\(\def\bupxi{\unicode[Times]{x1D6CF}}\)\(\def\bupomicron{\unicode[Times]{x1D6D0}}\)\(\def\buppi{\unicode[Times]{x1D6D1}}\)\(\def\buprho{\unicode[Times]{x1D6D2}}\)\(\def\bupsigma{\unicode[Times]{x1D6D4}}\)\(\def\buptau{\unicode[Times]{x1D6D5}}\)\(\def\bupupsilon{\unicode[Times]{x1D6D6}}\)\(\def\bupphi{\unicode[Times]{x1D6D7}}\)\(\def\bupchi{\unicode[Times]{x1D6D8}}\)\(\def\buppsy{\unicode[Times]{x1D6D9}}\)\(\def\bupomega{\unicode[Times]{x1D6DA}}\)\(\def\bupvartheta{\unicode[Times]{x1D6DD}}\)\(\def\bGamma{\bf{\Gamma}}\)\(\def\bDelta{\bf{\Delta}}\)\(\def\bTheta{\bf{\Theta}}\)\(\def\bLambda{\bf{\Lambda}}\)\(\def\bXi{\bf{\Xi}}\)\(\def\bPi{\bf{\Pi}}\)\(\def\bSigma{\bf{\Sigma}}\)\(\def\bUpsilon{\bf{\Upsilon}}\)\(\def\bPhi{\bf{\Phi}}\)\(\def\bPsi{\bf{\Psi}}\)\(\def\bOmega{\bf{\Omega}}\)\(\def\iGamma{\unicode[Times]{x1D6E4}}\)\(\def\iDelta{\unicode[Times]{x1D6E5}}\)\(\def\iTheta{\unicode[Times]{x1D6E9}}\)\(\def\iLambda{\unicode[Times]{x1D6EC}}\)\(\def\iXi{\unicode[Times]{x1D6EF}}\)\(\def\iPi{\unicode[Times]{x1D6F1}}\)\(\def\iSigma{\unicode[Times]{x1D6F4}}\)\(\def\iUpsilon{\unicode[Times]{x1D6F6}}\)\(\def\iPhi{\unicode[Times]{x1D6F7}}\)\(\def\iPsi{\unicode[Times]{x1D6F9}}\)\(\def\iOmega{\unicode[Times]{x1D6FA}}\)\(\def\biGamma{\unicode[Times]{x1D71E}}\)\(\def\biDelta{\unicode[Times]{x1D71F}}\)\(\def\biTheta{\unicode[Times]{x1D723}}\)\(\def\biLambda{\unicode[Times]{x1D726}}\)\(\def\biXi{\unicode[Times]{x1D729}}\)\(\def\biPi{\unicode[Times]{x1D72B}}\)\(\def\biSigma{\unicode[Times]{x1D72E}}\)\(\def\biUpsilon{\unicode[Times]{x1D730}}\)\(\def\biPhi{\unicode[Times]{x1D731}}\)\(\def\biPsi{\unicode[Times]{x1D733}}\)\(\def\biOmega{\unicode[Times]{x1D734}}\)\begin{equation}\tag{1}\frac{\partial c}{\partial t}=\underbrace{\frac{D}{{{R}^{2}}\sin (\theta )}\left( \frac{\partial }{\partial \theta }\left( \sin (\theta )\frac{\partial c}{\partial \theta } \right)+\frac{1}{\sin (\theta )}\frac{{{\partial }^{2}}c}{\partial {{\phi }^{2}}} \right)}_{\rm{diffusion}}-\underbrace{\frac{\alpha Qpc}{\gamma +c}}_{\rm{uptake}}+\underbrace{\beta h({{c}_{ch}}-c)}_{ \rm{{exchange}\ {with}} \atop \rm{{ choriocapillaris}} }{\rm{,}}\end{equation}

 

\begin{equation}\tag{2}\frac{\partial p}{\partial t}=\underbrace{\mu p\left( 1-\frac{p}{\tilde{p}(\theta )} \right){{{\hat{\lambda }}}_{1}}(c)}_{\rm{regrowth\ (normoxia)}}-\underbrace{\delta p{{{\hat{\lambda }}}_{2}}(c)}_{ \rm{{degeneration}} \atop \rm{{(hyperoxia)}}} {\rm{,}}\end{equation}

where ∂c/∂t and ∂p/∂t are the rate of change of the oxygen concentration and photoreceptor density respectively over time. The variable p(θ,ϕ,t) can also be thought of as representing OS biomass density (multiplied by the appropriate scaling factor to maintain dimensional consistency), this interpretation being more appropriate in those cases where regrowth occurs since new photoreceptors cannot be regrown in the adult retina (see Table 1 for a summary of the model variables).  

Table 1

 

Summary of the Physical Variables Used in the Mathematical Models (see Equations 1–14)

 

We note that D (m²s^–1) is the (constant⁶³) diffusivity of oxygen and R (m) is the radial position of the retina. Oxygen consumption is assumed to have Michaelis-Menten kinetics (an increasing, saturating function of the retinal oxygen concentration), where Q (mol s^–1 [tissue unit]^–1) is the maximum rate of oxygen consumption under light-adapted conditions (when consumption is lowest and the risk of hyperoxic damage is highest),⁶⁴ assuming that rods and cones have the same oxygen demand,⁶⁵ and γ (mol m^–3) is the oxygen concentration at which oxygen uptake is half maximal. The parameter α (m^–1) is the ratio of unit surface area to unit volume (ensuring dimensional consistency), c_ch (mol m^–3) is the oxygen concentration in the CC, β (m s^–1) is the effective permeability of the CC vessels and BM to oxygen, h (m^–1) is the capillary surface area per unit volume of tissue and δ (s^–1) is the rate of hyperoxia-induced photoreceptor degeneration. We remark that the rates of oxygen supply and uptake are large enough to maintain a heterogeneous oxygen profile despite the smoothing effect of diffusion. It has been necessary to reduce the values of both Q and β by a factor of 10 from those used in Roberts et al.⁵⁹ to render the simulations computationally feasible. This results in a faster rate of propagation of hyperoxic degeneration, but does not alter the spatial pattern that develops, which is our focus here (see Roberts⁶⁰ for details). OS biomass regrowth is assumed to occur logistically with intrinsic growth rate μ (s^–1) and carrying capacity Display Formula\(\tilde p(\theta )\) (photoreceptors m^–2), where  

\begin{equation}\tag{3}\tilde p(\theta ) = \underbrace {{B_1}{e^{ - {b_1}\theta }} + {B_2}{e^{ - {b_2}\theta }}}_{{\rm{cones}}} + \underbrace {{B_3}\theta {e^{ - {b_3}\theta }}}_{{\rm{rods}}}{.}\end{equation}

This functional form was chosen because of its qualitative similarity to the healthy human photoreceptor distribution, where the first two terms capture the cone profile and the last term captures the rod profile. The values for the parameters B₁ (photoreceptors m^–2), B₂ (photoreceptors m^–2), B₃ (photoreceptors m^–2 rad^–1), b₁ (rad^–1), b₂ (rad^–1), and b₃ (rad^–1) were obtained by using the Matlab curve fitting toolbox (MathWorks, Natick, MA, USA) to fit Equation 3 to the mean of eight experimentally measured healthy adult photoreceptor distributions across the temporal horizontal meridian, using data provided by Curcio et al.⁶² (see Table 2). We note that, for simplicity, we assume that the photoreceptor distribution depends only upon θ and hence is axisymmetric about the z-axis (see Fig. 3).  

Table 2

 

Parameter Values for the Three Mathematical Models (see Equations 1–7; Equations 8–12; and Equations 13,14), Given to a Maximum of Three Significant Figures

 

The functions Display Formula\({\hat \lambda _1}(c)\) (dimensionless) and Display Formula\({\hat \lambda _2}(c)\) (dimensionless) that appear in Equation 2 are defined as  

\begin{equation}\tag{4}\displaylines{ {{\hat \lambda }_1}(c) = 1 - {1 \over 2}(\tanh ({S_c}(c - {c_{crit}})) + 1){\rm{,}} \cr {{\hat \lambda }_2}(c) = 1 - {{\hat \lambda }_1}(c) = {1 \over 2}(\tanh ({S_c}(c - {c_{crit}})) + 1) \approx \left\{ {\matrix{ 0&{{\rm{if\ }}c \lt {c_{crit}}{\rm{,}}} \cr 1&{{\rm{if\ }}c \ge {c_{crit}}{\rm{,}}} \cr } } \right. \cr} \end{equation}

wherein c_crit (mol m^–3) represents the ‘hyperoxic threshold', above which photoreceptors degenerate and below which they remain healthy or regenerate biomass. The parameter S_c (rad m³mol^–1) determines the sharpness of the transition between normoxia and hyperoxia. We use hyperbolic (tanh) functions here, rather than the Heaviside step functions used in our earlier work,⁵⁹ in order to make the governing equations smooth, improving the convergence of the numerical scheme. This remains a biologically realistic choice for these terms.  

We close the system by imposing the following boundary and initial conditions:  

\begin{equation}\tag{5}{{\partial c} \over {\partial \theta }}(\Theta ,\phi ,t) = 0{\rm{,}}\end{equation}

 

\begin{equation}\tag{6}c(\theta ,\phi ,0) = {c_{init}}(\theta ,\phi ), \quad p(\theta ,\phi ,0) = {p_{init}}(\theta ,\phi ) = F(\theta ,\phi )\tilde p(\theta ){.}\end{equation}

The boundary conditions are zero-flux, such that there is no net gain or loss of oxygen to or from the retina along its boundary at the ora serrata. The initial oxygen distribution, c_init(θ,ϕ) (mol m^–3), is the steady-state oxygen concentration (the distribution to which oxygen settles over time) corresponding to the initial photoreceptor profile (that is, the steady-state solution to Equations 1 and 5 with p(θ,ϕ) = p_init(θ,ϕ) [photoreceptors m^–2]). The function 0 ≤ F(θ,ϕ) ≤ 1 (dimensionless) is used to remove a patch of photoreceptors and takes one of the following three forms:  

\begin{equation}\tag{7}\matrix{ {F(\theta ,\phi ) = \left\{ {\matrix{ {} \hfill&{1{\rm{,\ or,}}} \hfill \cr {} \hfill&{{1 \over 2}(\tanh (S({\theta _1} - \theta )) + \tanh (S(\theta - {\theta _2})) + 2){\rm{,\ or,}}} \hfill \cr {} \hfill&{{1 \over 2}(\tanh (\hat S({{(\theta - {\theta _c})}^2} + {{(\phi - {{\hat \phi }_c})}^2}{{\sin }^2}(\theta ) - {\psi ^2})) + 1){\rm{,}}} \hfill \cr } } \right.} \hfill \cr } \end{equation}

which correspond to a healthy retina, the removal of an annulus (ring) of photoreceptors and the removal of a disc of photoreceptors respectively (see Supplementary Fig. S1 for graphs depicting each of these initial conditions). Degenerate discs represent the patchy loss observed in the early stages of RP, while annuli represent a later disease stage. The parameters S (dimensionless) and Display Formula\(\hat S\) (rad^–1) determine the sharpness of the edge of an annulus or disc respectively, while θ₁ (rad) and θ₂ (rad) are the eccentricities of the inner and outer boundaries of the annulus. The disc is centered at Display Formula\((\theta ,\phi ) = ({\theta _c},{\hat \phi _c})\) (rad) (where we write the azimuthal position with a ‘hat' to distinguish it from the rate of mutation-induced cone degeneration below) with radius governed by ψ (rad). The sin²(θ) term is required in order to preserve arc length in the azimuthal direction. We choose hyperbolic (tanh) functions, rather than step functions (as used in our earlier work)⁵⁹ in order to smooth the transition between the healthy and degenerate retina, so that a coarser finite element mesh can be used to resolve the initial wavefront of degeneration. See Table 2 for parameter values.  

Mutation-induced Rod Degeneration

Modifying our model to account for mutation-induced rod degeneration (corresponding to the rod–cone form of RP), we denote rod density by p_r(θ,ϕ,t) (photoreceptors m^–2) and cone density by p_c(θ,ϕ,t) (photoreceptors m^–2), to obtain the following system:  

\begin{equation}\tag{8}{{\partial c} \over {\partial t}} = {D \over {{R^2}\sin (\theta )}}\left( {{\partial \over {\partial \theta }}\left( {\sin (\theta ){{\partial c} \over {\partial \theta }}} \right) + {1 \over {\sin (\theta )}}{{{\partial ^2}c} \over {\partial {\phi ^2}}}} \right) - {{\alpha Q({p_r} + {p_c})c} \over {\gamma + c}} + \beta h({c_{ch}} - c){\rm{,}}\end{equation}

 

\begin{equation}\tag{9}{{\partial {p_r}} \over {\partial t}} = - ({\phi _r} + \delta {\hat \lambda _2}(c)){p_r}{\rm{,}}\end{equation}

 

\begin{equation}\tag{10}{{\partial {p_c}} \over {\partial t}} = \mu {p_c}\left( {1 - {{p_c} \over {{{\tilde p}_c}(\theta )}}} \right){\hat \lambda _1}(c) - \delta {p_c}{\hat \lambda _2}(c){\rm{,}}\end{equation}

where ϕ_r (s^–1) is the rate of mutation-induced rod degeneration and rods are assumed to be incapable of recovering biomass. We close Equations 8 to 10 by imposing zero-flux oxygen boundary conditions at θ = Θ (rad) (given by Equation 5) and the following initial conditions:  

\begin{equation}\tag{11}\eqalign{ c(\theta ,\phi ,0) &= {c_{init}}(\theta ,\phi ){\rm{,}} \cr {p_r}(\theta ,\phi ,0) &= {p_{{r_{init}}}}(\theta ,\phi ) = F(\theta ,\phi ){{\tilde p}_r}(\theta ){\rm{,}} \cr {p_c}(\theta ,\phi ,0) &= {p_{{c_{init}}}}(\theta ,\phi ) = F(\theta ,\phi ){{\tilde p}_c}(\theta ){\rm{,}} \cr} \end{equation}

where F(θ,ϕ) is defined in Equation 7 and the healthy rod and cone distributions, Display Formula\({\tilde p_r}(\theta )\)(photoreceptors m^–2) and Display Formula\({\tilde p_c}(\theta )\) (photoreceptors m^–2), are given by  

\begin{equation}\tag{12}{\tilde p_r}(\theta ) = {B_3}\theta {e^{ - {b_3}\theta }},\qquad {\tilde p_c}(\theta ) = {B_1}{e^{ - {b_1}\theta }} + {B_2}{e^{ - {b_2}\theta }}{.}\end{equation}

As before, the initial oxygen concentration, c_init(θ,ϕ), is the steady-state oxygen concentration corresponding to the initial photoreceptor profile (that is, the steady-state solution to Equations 5 and 8 with Display Formula\({p_r}(\theta ,\phi ) = {p_{{r_{init}}}}(\theta ,\phi )\)[photoreceptors m^–2] and Display Formula\({p_c}(\theta ,\phi ) = {p_{{c_{init}}}}(\theta ,\phi )\) [photoreceptors m^–2]).  

Mutation-induced Cone Degeneration

Similarly, we may account for mutation-induced cone degeneration (corresponding to the cone–rod form of RP), obtaining Equation 8 together with  

\begin{equation}\tag{13}{{\partial {p_r}} \over {\partial t}} = \mu {p_r}\left( {1 - {{p_r} \over {{{\tilde p}_r}(\theta )}}} \right){\hat \lambda _1}(c) - \delta {p_r}{\hat \lambda _2}(c){\rm{,}}\end{equation}

 

\begin{equation}\tag{14}{{\partial {p_c}} \over {\partial t}} = - ({\phi _c} + \delta {\hat \lambda _2}(c)){p_c}{\rm{,}}\end{equation}

where ϕ_c (s^–1) is the rate of mutation-induced cone degeneration. The value chosen for ϕ_c is twice that used for ϕ_r in the present study (see Table 2), and twice that used for ϕ_r and ϕ_c in Roberts et al.,⁵⁹ in order to reduce the computation time, while remaining within the range of biologically realistic values for this parameter, the spatial pattern of degeneration remaining the same as for the lower value. We close the system by imposing the boundary and initial conditions given by Equations 5 and 11 respectively.  

Treatment With Antioxidants or Trophic Factors

For each of the cases above, we may apply treatment in the form of antioxidants or trophic factors. These treatments are assumed to increase the hyperoxic threshold, such that c_crit becomes a function of time, t. We set Display Formula\({c_{crit}}(t) = {c_{cri{t_1}}} + {c_{cri{t_2}}}H(t - {t_{crit}})\), where Display Formula\({c_{cri{t_1}}}\) (mol m^–3), Display Formula\({c_{cri{t_2}}}\) (mol m^–3) and t_crit (s) are positive constants and H is a Heaviside step function, such that Display Formula\({c_{crit}}(t) = {c_{cri{t_1}}}\) for t < t_crit and Display Formula\({c_{crit}}(t) = {c_{cri{t_1}}} + {c_{cri{t_2}}}\) for t ≥ t_crit. We use a Heaviside function rather than a hyperbolic (tanh) function here because smoothing is not required in this case for numerical stability. 

Results

Degenerate Annulus

Our model predicts that the photoreceptor distribution may evolve in one of four ways following the removal of an annulus (ring) of photoreceptors, in the absence of mutation-induced rod or cone degeneration (see Supplementary Fig. S2). The degenerate annulus may (1) remain stationary (Supplementary Fig. S2a), (2) spread centrally (Supplementary Fig. S2b), (3) spread peripherally (Supplementary Fig. S2c), or (4) expand both centrally and peripherally (Supplementary Fig. S2d). The evolution of the annulus is determined by the initial positions of the inner and outer boundaries of the annulus. Boundaries remain stationary where they border regions with a high local photoreceptor density and propagate outward where they border regions with a low local photoreceptor density. We note that these results are consistent with earlier results from a simpler 1D model (see Roberts et al.⁵⁹). 

Degenerate Disc

Simulations predict that, in the absence of mutation-induced rod or cone degeneration, the evolution of the photoreceptor distribution following the removal of a disc of photoreceptors is determined mainly by the disc's eccentricity, θ_c (where θ_c = 0 [rad] at the foveal center and θ_c = Θ [rad] at the ora serrata), and, to a lesser extent, by its radius, governed by ψ (larger values of ψ corresponding to a larger radius, see Fig. 4a). (See Roberts et al.⁵⁹ for an asymptotic analysis showing that the qualitative behavior of the model will be maintained within realistic parameter ranges due to the separation of scales between various parameter groupings.) As with the case of a degenerate annulus, a degenerate disc expands where its boundary borders regions in which the local photoreceptor density is low, and remains stationary otherwise. Photoreceptors recover for 0.3 × Θ ≤ θ_c ≤ 0.4 × Θ (rad) (Fig. 4d), while photoreceptor loss within the disc, which is only partial initially, is consolidated (that is photoreceptor loss becomes complete in this region) for θ_c = 0.7 × Θ (rad) (Fig. 4e) for all disc radii examined. The disc expands into the central retina for θ_c = 0 (rad) (Fig. 4b) and around the retinal periphery for θ_c > 0.7 × Θ (rad) (Fig. 4f) for all disc radii examined. Degeneration is consolidated within the disc for ψ = 0.025 × Θ (rad) and ψ = 0.0125 × Θ (rad) when θ_c = 0.1 × Θ (rad); however, when ψ = 0.05 × Θ (rad), degeneration spreads into the central retina (Fig. 4c). Lastly, the disc either completely recovers, consolidates completely or recovers in some areas while consolidating in others for θ_c = 0.2 × Θ (rad), θc = 0.5 × Θ (rad) and θc = 0.6 × Θ (rad), depending upon its radius. 

Figure 4

 

Sensitivity analysis demonstrating the degenerate retinal area resulting from an initial disc of photoreceptor loss. (a) Graph to show the effects of degenerate disc eccentricity, θ_c, and radius, governed by ψ, upon the resultant degenerate area at t = 3 × 10⁸ (s) (≈ 9.5 [yr]). The system has reached steady-state by this time for θ_c ≤ 0.7 × Θ (rad) and is close to steady-state for θ_c > 0.7 × Θ (rad). The inset shows a magnification of the bottom-left corner of the graph, demonstrating the degenerate area resulting from discs for which θ_c/Θ ≤ 0.2. An area is considered degenerate when photoreceptors have dropped beneath 10% of their average density, where p_average = 1.11 × 10⁵ (photoreceptors mm^–2). Degeneration is most severe for θ_c < 0.2 × Θ (rad) and θ_c > 0.7 × Θ (rad), while degeneration recovers for 0.3 × Θ ≤ θ_c ≤ 0.4 × Θ (rad). Degeneration is also more severe for larger values of ψ. (b–f) Panels show snapshots of the simulations at different stages in the disease progression as the initial radius and location of the degenerate patch vary. The plots show the projection of the solution from the retinal surface onto the (x,y) -plane (cm) (see Fig. 3b). Note that plots are magnified in (b, c) to clarify the dynamics. (b) θ_c = 0 (rad), (c) θ_c = 0.1 × Θ (rad), (d) θ_c = 0.3 × Θ (rad), (e) θ_c = 0.7 × Θ (rad), and (f) θ_c = 0.8 × Θ (rad). (b–e) ψ = 0.05 × Θ (rad) and (f) ψ = 0.0125 × Θ (rad). See Table 2 for the remaining parameter values.

Mutation-induced Rod Degeneration

Our model simulations suggest that mutation-induced rod degeneration will cause hyperoxic photoreceptor degeneration to initiate first at the ora serrata and later in the parafoveal/perifoveal region (see Fig. 5). Waves of hyperoxic degeneration propagate centrally from the ora serrata and peripherally from the macula, meeting in the midperiphery, such that only the foveal cones remain at steady-state. 

Figure 5

 

Simulation results illustrating the pattern of retinal degeneration associated with mutation-induced rod degeneration. Panels show the evolving rod (a) and cone (b) distributions over time. Note that we plot the natural logarithm of the (dimensionless) scaled cone density, \({\tilde p_c} = {p_c}{\rm{/}}{p_{average}}\), where p_average = 1.11 × 10⁵ (photoreceptors mm^–2), in (b), with the lower bound of the color range fixed to the minimum value of \(\log ({\tilde p_c}(\theta ,\phi ,t))\) at t = 0 (yr) to improve visual clarity. Hyperoxia-induced photoreceptor loss initiates first at the ora serrata and later in the parafoveal/perifoveal region, spreading into the midperiphery. Only the central cone island remains at steady-state. See Table 2 for parameter values.

Simulations further predict that mutation-induced rod degeneration can cause expansion of degenerate regions such as annuli and discs that would otherwise have remained fixed in size. Comparison of the results presented in Figures 4e and 6 illustrates this behavior: a degenerate disc which remains fixed in size in the absence of mutation-induced rod degeneration (see Fig. 4e) expands in size when mutation-induced rod degeneration is active (see Fig. 6). 

Figure 6

 

Simulation results showing how mutation-induced rod degeneration can stimulate expansion of a pre-existing degenerate disc. Disc eccentricity θ_c = 0.7 × Θ (rad) and radius parameter ψ = 0.05 × Θ (rad) are the same as in Figure 4e, where the disc does not expand. Panels show the evolving rod (a) and cone (b) distributions over time. (See Fig. 5 legend for details of the plotting procedure.) See Table 2 for the remaining parameter values.

Mutation-induced Cone Degeneration

Simulations suggest that mutation-induced cone degeneration will stimulate hyperoxia-driven degeneration of photoreceptors in the central retina (see Fig. 7). Degeneration initiates in the fovea and spreads peripherally, coming to a halt as it enters the midperiphery. 

Figure 7

 

Simulation results illustrating the pattern of retinal degeneration associated with mutation-induced cone degeneration. Panels show the evolving rod (a) and cone (b) distributions over time. (See Fig. 5 legend for details of the plotting procedure.) Hyperoxia-induced photoreceptor loss initiates in the fovea, spreading outward and stalling as it enters the midperiphery. See Table 2 for parameter values.

Treatment With Antioxidants or Trophic Factors

Our models predict that treatment with antioxidants and/or trophic factors may either delay, halt, or partially reverse hyperoxia-induced photoreceptor degeneration depending upon the strength (Display Formula\({c_{cri{t_2}}}\)) and timing (t_crit) of treatment, the extent of the degenerate region(s) when treatment is applied and the presence or absence of mutation-induced rod or cone degeneration (see Fig. 8). In the absence of mutation-induced rod or cone degeneration, treatment may halt degeneration (Figs. 8a, 8c) or enable partial recovery (Fig. 8b). By contrast, in the presence of mutation-induced cone degeneration treatment may at best halt (Fig. 8e) degeneration and in the presence of mutation-induced rod degeneration delay it (Fig. 8d). 

Figure 8

 

Graphs to show the effect of treatment with antioxidants and/or trophic factors upon the degenerate retinal area. An area is considered degenerate when photoreceptors have dropped beneath 10% of their average density, where p_average = 1.11 × 10⁵ (photoreceptors mm^–2). (a) Degenerate annulus with (θ₁,θ₂) = (0.08, 0.9) ×Θ (rad), (b) degenerate disc with (θ_c,ψ) = (0.1, 0.05) × Θ (rad), (c) degenerate disc with (θ_c,ψ) = (0.8, 0.0125) × Θ (rad), (d) mutation-induced rod degeneration, and (e) mutation-induced cone degeneration. Treatment is predicted to delay the spread of degeneration in (d), halt its spread in (a, c, e), and partially reverse degeneration in (b). We note that the step-like appearance of the curves in some panels is the result of multiple regions falling below the p_average threshold simultaneously. This effect is more pronounced for radially-symmetric solutions (a, d, e) and is exaggerated in (e) due to the magnification of the y-axis. Treatment is applied at t_crit = 0.63 (yr) in (a), t_crit = 1.27 (yr) in (b), t_crit = 3.17 (yr) in (c), t_crit = 9.51 (yr) in (d), and t_crit = 15.84 (yr) in (e). See Table 2 for the remaining parameter values.

Discussion

Mathematical modeling enables us to isolate biological mechanisms in a manner that would be difficult or impossible experimentally. In this paper we have used spatially resolved 2D models to predict those patterns of retinal degeneration that could be produced if oxygen toxicity were the sole cause for the spread of photoreceptor degeneration in RP. Comparing our results against the taxonomy of visual field loss patterns developed by Grover et al.,¹³ we find that our models recapitulate several patterns. Mutation-induced rod degeneration is predicted to produce a pattern in which degeneration initiates at the periphery and in the parafoveal/perifoveal regions and spreads into the midperiphery, similar to Pattern 1 and especially Pattern 1B. Likewise, if a patch of photoreceptor loss protrudes far enough into the peripheral retina, degeneration sweeps around the periphery of the retina, similar to the mid–late stage of Pattern 3. Patterns that involve a preferential loss of midperipheral photoreceptors, such as Pattern 2 and the initial stage of Pattern 3, are theoretically inaccessible to our models and hence to the oxygen toxicity mechanism. This is because the high density of photoreceptors and hence high rate of oxygen consumption in this region causes the propagation of hyperoxia-induced photoreceptor degeneration to halt (see Roberts et al.⁵⁹ for a detailed analysis). For the same reason, the central cone island is preserved, in agreement with Grover et al.,¹³ though its eventual loss cannot be explained unless mutation-induced cone loss is active. Similarly, the hyperoxic photoreceptor degeneration triggered by mutation-induced cone loss stalls at it enters the midperiphery because this region is rod-dominated. Thus our theoretical models demonstrate, in a novel manner, the strengths and weaknesses of the oxygen toxicity hypothesis: while this mechanism is sufficient to generate some of the patterns of degeneration seen in vivo, other mechanisms are necessary to explain the remaining patterns. 

Our models illustrate how treatment with antioxidants and/or trophic factors could halt (or partially reverse) the spread of photoreceptor degeneration in those cases where a patch of photoreceptors has been lost and mutation-induced rod and cone degeneration are inactive. In these cases the treatments are effective because they raise the hyperoxic threshold and, hence, decrease the critical photoreceptor density at which degeneration ceases to propagate. As a result, degenerate patches may remain stable, without expanding, in a larger proportion of the retina. While partial recovery may be feasible where photoreceptor IS are still intact, regeneration is likely to be limited, as many photoreceptors will have degenerated past the point of recovery. In the case of mutation-induced rod degeneration, treatment can at best delay degeneration, as the steady attrition of rods will eventually bring the local photoreceptor density below the critical level required for the propagation of hyperoxic degeneration. Treatment may, however, halt hyperoxic degeneration associated with mutation-induced cone degeneration. This is because cone density is low outside of the fovea, the remainder of the retina being rod-dominated. These results suggest that, where possible, antioxidant/trophic factor treatment should be targeted to regions of the retina in which the photoreceptor density is low (e.g., using ocular gene therapy^4,43) as these regions are most susceptible to hyperoxic degeneration. 

We plan to extend this work in a number of directions in future studies. We will incorporate heterogeneity in photoreceptor density in the azimuthal dimension and consider patterns of degeneration in other mammals such as rats, mice, rabbits, cats, pigs, and monkeys. This will allow us to explore further the strengths and weaknesses of the oxygen toxicity hypothesis and may help us to explain the distinct spatio-temporal patterns of photoreceptor degeneration associated with different species. We will also consider a 2D model, spanning the depth of the retina in 1D (as in Roberts et al.⁷⁴) and the region between the foveal center and the ora serrata in the other (neglecting the azimuthal dimension), to account for the oxygen distribution across the retinal depth and the effects of photoreceptor OS shrinkage upon the propagation of retinal degeneration. Lastly, we will use the present modeling framework to consider other disease mechanisms in RP, in particular the trophic factor, toxic substance, microglia, and mTOR hypotheses, both in isolation and in combination. 

More experimental and clinical data are needed to validate and inform future modeling studies. The most useful data would be those derived from more detailed longitudinal studies in animals and humans, recording the spatio-temporal development of degenerate patches (using optical coherence tomography) and photoreceptor topography (using adaptive optics scanning laser ophthalmoscopy) across the retina. This would enable us to conduct a more detailed quantitative validation of our models. Measurements of the oxygen distribution across the retinal depth at a range of eccentricities and azimuthal angles at different stages of the disease would also be valuable in validating and extending our models, as would detailed quantitative data describing the effects of antioxidant and trophic factor treatment over time. 

Our theoretical study has clarified the explanatory scope of the oxygen toxicity hypothesis as a disease mechanism underlying RP. Hyperoxia may well play a role in the spread of retinal degeneration, but other mechanisms are required to complete the picture. Our models further demonstrate the potential of antioxidant and trophic factor therapies in reducing hyperoxic disease progression, while highlighting the complementary insights that can be attained by combining a mathematical modeling approach with experimental and clinical studies. 

Acknowledgments

Supported by grants from Engineering and Physical Sciences Research Council (EPSRC) in Swindon, United Kingdom for funding through a studentship at the Systems Biology programme of the University of Oxford's Doctoral Training Centre (PAR; EP/G500029/1 and EP/G03706X/1). 

Disclosure: P.A. Roberts, None; E.A. Gaffney, None; J.P. Whiteley, None; P.J. Luthert, None; A.J.E. Foss, None; H.M. Byrne, None