**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.

^{1}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.

^{3}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.

^{4}In this paper, we use mathematical models to investigate one such hypothesis, namely the oxygen toxicity hypothesis.

^{5–7}

^{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.

^{13}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}

^{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,

^{31}microglia,

^{32}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}

^{37}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}

^{48–54}mTOR,

^{54}and toxic substance

^{55–57}hypotheses (see Roberts et al.

^{58}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.

^{59}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.,

^{13}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.

^{59}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

^{60}for further details).

^{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.

*θ*= 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]

^{61}). The ora serrata is assumed to lie at

*θ*= Θ ≈ 1.33 (rad), based on data supplied by Curcio et al.

^{62}(noting that, unlike the retina depicted in Fig. 3a, the retinas measured by Curcio et al.

^{62}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.

*c*(

*θ*,

*ϕ*,

*t*) (mol m

^{–3}), and the photoreceptor density,

*p*(

*θ*,

*ϕ*,t) (photoreceptors m

^{–2}), over time,

*t*(s), as follows

*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).

*D*(m

^{2}

*s*

^{–1}) is the (constant

^{63}) 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),

^{64}assuming that rods and cones have the same oxygen demand,

^{65}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*(mol m

_{ch}^{–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.

^{59}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

^{60}for details). OS biomass regrowth is assumed to occur logistically with intrinsic growth rate

*μ*(s

^{–1}) and carrying capacity

^{–2}), where

*B*

_{1}(photoreceptors m

^{–2}),

*B*

_{2}(photoreceptors m

^{–2}),

*B*

_{3}(photoreceptors m

^{–2}rad

^{–1}),

*b*

_{1}(rad

^{–1}),

*b*

_{2}(rad

^{–1}), and

*b*

_{3}(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.

^{62}(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**

*c*(mol m

_{crit}^{–3}) represents the ‘

*hyperoxic threshold*', above which photoreceptors degenerate and below which they remain healthy or regenerate biomass. The parameter

*S*(rad m

_{c}^{3}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,

^{59}in order to make the governing equations smooth, improving the convergence of the numerical scheme. This remains a biologically realistic choice for these terms.

*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:

*S*(dimensionless) and

^{–1}) determine the sharpness of the edge of an annulus or disc respectively, while

*θ*

_{1}(rad) and

*θ*

_{2}(rad) are the eccentricities of the inner and outer boundaries of the annulus. The disc is centered at

*ψ*(rad). The sin

^{2}(

*θ*) 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)

^{59}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.

*p*(

_{r}*θ*,

*ϕ*,

*t*) (photoreceptors m

^{–2}) and cone density by

*p*(

_{c}*θ*,

*ϕ*,

*t*) (photoreceptors m

^{–2}), to obtain the following system:

*ϕ*(s

_{r}^{–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:

*F*(

*θ*,

*ϕ*) is defined in Equation 7 and the healthy rod and cone distributions,

^{–2}) and

^{–2}), are given by

*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

^{–2}] and

^{–2}]).

*ϕ*(s

_{c}^{–1}) is the rate of mutation-induced cone degeneration. The value chosen for

*ϕ*is twice that used for

_{c}*ϕ*in the present study (see Table 2), and twice that used for

_{r}*ϕ*and

_{r}*ϕ*in Roberts et al.,

_{c}^{59}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.

*c*becomes a function of time,

_{crit}*t*. We set

^{–3}),

^{–3}) and

*t*(s) are positive constants and

_{crit}*H*is a Heaviside step function, such that

*t*<

*t*and

_{crit}*t*≥

*t*. We use a Heaviside function rather than a hyperbolic (tanh) function here because smoothing is not required in this case for numerical stability.

_{crit}^{59}).

*θ*(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

_{c}*ψ*(larger values of

*ψ*corresponding to a larger radius, see Fig. 4a). (See Roberts et al.

^{59}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 × Θ ≤

*θ*≤ 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

_{c}*ψ*= 0.025 × Θ (rad) and

*ψ*= 0.0125 × Θ (rad) when

*θ*= 0.1 × Θ (rad); however, when

_{c}*ψ*= 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

*θ*= 0.2 × Θ (rad),

_{c}*θc*= 0.5 × Θ (rad) and

*θc*= 0.6 × Θ (rad), depending upon its radius.

*t*) 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).

_{crit}^{13}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.

^{59}for a detailed analysis). For the same reason, the central cone island is preserved, in agreement with Grover et al.,

^{13}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.

^{4,43}) as these regions are most susceptible to hyperoxic degeneration.

^{74}) 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.

**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

*. 2009; 80: 384–401.*

*Optometry**. 2006; 1: 40.*

*Orphanet J Rare Dis**. 2006; 368: 1795–1809.*

*Lancet**. 2011; 2011: 753547.*

*J Ophthalmol**. 1999; 18: 689–735.*

*Prog Retin Eye Res**. 1991; 6: 61–70.*

*Neuron**. 1998; 39: 2427–2442.*

*Invest Ophthalmol Vis Sci**. 1998; 95: 7103–7108.*

*Proc Natl Acad Sci U S A**. 2013; 54: 5888–5900.*

*Invest Ophthalmol Vis Sci**. 2012; 520: 874–888.*

*J Comp Neurol**. 2011; 59: 1107–1117.*

*Glia**. 2013; 351: 29–40.*

*Cell Tissue Res**, 1998; 105: 1069–1075.*

*Ophthalmology**. 2012; 53: 4754–4764.*

*Invest Ophthalmol Vis Sci**. 2009; 29: 1025–1031.*

*Retina**. 2012; 153: 718–727.*

*Am J Ophthalmol**. 2008; 145: 687–694.*

*Am J Ophthalmol**. 2005; 243: 1018–1027.*

*Graefes Arch Clin Exp Ophthalmol**. 2003; 44: 3544–3550.*

*Invest Ophthalmol Vis Sci**. 2004; 45: 4119–4125.*

*Invest Ophthalmol Vis Sci**. 2006; 90: 472–479.*

*Br J Ophthalmol**. 2008; 116: 79–89.*

*Doc Ophthalmol**. 2011; 10: 1–10.*

*Retina**. 1997; 29: 290–297.*

*Ophthalmic Res**. 1998; 95: 8357–8362.*

*Proc Natl Acad Sci U S A**. 2000; 118: 807–811.*

*Arch Ophthalmol**. 2003; 44: 818–825.*

*Invest Ophthalmol Vis Sci**. 2004; 36: 755–759.*

*Nat Genet**. 2015; 161: 817–832.*

*Cell**. 2016; 24: 909–923.*

*Antioxid Redox Signal**. 2002; 74: 327–336.*

*Exp Eye Res**. 2003; 76: 463–471.*

*Exp Eye Res**. 2009; 12: 44–52.*

*Nat Neurosci**. 2012; 287: 1642–1648.*

*J Biol Chem**. 2004; 45: 2013–2019.*

*Invest Ophthalmol Vis Sci**. 2005; 80: 745–751.*

*Exp Eye Res**. 2002; 30: 620–650.*

*Toxicol Pathol**. 2006; 103: 11300–11305.*

*Proc Natl Acad Sci U S A**. 2007; 213: 809–815.*

*J Cell Physiol**. 2002; 43: 3292–3298.*

*Invest Ophthalmol Vis Sci**. 2004; 10: 555–565.*

*Mol Vis**. 2006; 6: 717–726.*

*Expert Opin Biol Ther**. 2006; 208: 516–526.*

*J Cell Physiol**. 2012; 31: 136–151.*

*Prog Retin Eye Res**. 2003; 23: 4164–4172.*

*J Neurosci**. 2007; 145: 1120–1129.*

*Neuroscience**. 2001; 159: 1113–1120.*

*Am J Pathol**. 2003: BU-1640-M.*

*Biomnetrics Unit Technical Reports**. 2010; 267: 638–646.*

*J Theor Biol**. 2013; 317: 105–118.*

*J Theor Biol**. 2014; 76: 292–313.*

*Bull Math Biol**. 2016; 38: 267–276.*

*Commun Nonlinear Sci Numer Simulat**. 2016; 78: 1394–1409.*

*Bull Math Biol**. 2016; 408: 75–87.*

*J Theor Biol**. 2000; 406: 195–199.*

*Nature**. 2001; 10: 2269–2275.*

*Hum Mol Genet**. 2202; 64: 1117–1145.*

*Bull Math Biol**. 2016; 53: 48–69.*

*Prog Retin Eye Res**. 2017; 425: 53–71.*

*J Theor Biol**. Oxford, UK: University of Oxford, 2015.*

*Mathematical Models of the Retina in Health and Disease*[DPhil thesis]*. Sunderland: Sinauer Associates Inc., 1999.*

*The Human Eye: Structure and Function**. 1990; 292: 497–523.*

*J Comp Neurol**. 1990; 277: 127–136.*

*Adv Exp Med Biol**. 2003; 121: 547–557.*

*Arch Ophthalmol**. 2002; 19: 395–407.*

*Vis Neurosci**. 2007; 293: H1696–H1704.*

*Am J Physiol Heart Circ Physiol**. 1997; 273: C852–C858.*

*Am J Physiol**. 2008; 15: 795–811.*

*Microcirculation**. 2001; 91: 2255–2265.*

*J Appl Physiol**. 1999; 277: H1831–H1840.*

*Am J Physiol Heart Circ Physiol**. 1994; 267: H2498–H2507.*

*Am J Physiol Heart Circ Physiol**. 1990; 4: 303–309.*

*Eye**. 1993; 34: 3278–3296.*

*Invest Ophthalmol Vis Sci**. 2016; 73: 1–38.*

*J Math Biol*