Inverse Problem Reveals Conditions for Characteristic Retinal Degeneration Patterns in Retinitis Pigmentosa Under the Trophic Factor Hypothesis

Retinitis pigmentosa (RP) is the most common inherited retinal dystrophy with a prevalence of about 1 in 4,000, affecting approximately 1.5 million people worldwide. Patients with RP experience progressive visual field loss as the retina degenerates, destroying light-sensitive photoreceptor cells (rods and cones), with rods affected earlier and more severely than cones. Spatio-temporal patterns of retinal degeneration in human RP have been well characterised; however, the mechanism(s) giving rise to these patterns have not been conclusively determined. One such mechanism, which has received a wealth of experimental support, is described by the trophic factor hypothesis. This hypothesis suggests that rods produce a trophic factor necessary for cone survival; the loss of rods depletes this factor, leading to cone degeneration. In this article, we formulate a partial differential equation mathematical model of RP in one spatial dimension, spanning the region between the retinal centre (fovea) and the retinal edge (ora serrata). Using this model we derive and solve an inverse problem, revealing for the first time experimentally testable conditions under which the trophic factor mechanism will qualitatively recapitulate the spatio-temporal patterns of retinal regeneration observed in human RP.

1. Introduction

The group of inherited retinal diseases known as retinitis pigmentosa (RP) causes the progressive loss of visual function (Hamel, 2006; Hartong et al., 2006). The patterns of visual field loss associated with the human version of this condition have been well characterised (Grover et al., 1998); however, the mechanisms underpinning these patterns have yet to be conclusively determined (Newton and Megaw, 2020). In this article, we use mathematical models to predict the conditions under which a trophic factor mechanism could explain these patterns.

The retina is a tissue layer lining the back of the eye containing light-sensitive cells known as photoreceptors, which come in two varieties: rods and cones (Figure 1A). Rods confer monochromatic vision under low-light (scotopic) conditions, while cones confer colour vision under well-lit (photopic) conditions (Oyster, 1999). In RP, rod function and health are typically affected earlier and more severely than those of cones, with cone loss following rod loss. Rods are lost since either they or the neighbouring retinal pigment epithelium express a mutant version of one or both alleles (depending on inheritance mode) of a gene associated with RP (over 80 genes have been identified to date, see Gene Vision and Ge et al., 2015; Haer-Wigman et al., 2017; Birtel et al., 2018; Coussa et al., 2019). It is hypothesised that cones are lost following rods since they depend upon rods either directly or indirectly for their survival (Hamel, 2006; Hartong et al., 2006; Daiger et al., 2007).

FIGURE 1

Figure 1. Diagrams of the human eye and retinal photoreceptor distribution (reproduced, with permission, from Roberts et al., 2017). (A) Diagram of the (right) human eye, viewed in the transverse plane, illustrating the mathematical model geometry. The model is posed on a domain spanning the region between the foveal centre, at θ = 0, and the ora serrata, at θ = Θ, along the temporal horizontal meridian, where θ measures the eccentricity. Figure originally reproduced, with modifications, from http://www.nei.nih.gov/health/coloboma/coloboma.asp, courtesy: National Eye Institute, National Institutes of Health (NEI/NIH). (B) Measured and fitted photoreceptor profiles, along the temporal horizontal meridian, in the human retina. Cone profile: ${\tilde{p}}_c(\theta )=B_1{e}^{-b_1\theta }+B_2{e}^{-b_2\theta }$, and rod profile: ${\tilde{p}}_r(\theta )=B_3\theta e^{-b_3\theta }$ (see Table 2 for dimensionless parameter values). The photoreceptor profile is the sum of the rod and cone profiles (${\tilde{p}}_r(\theta )+{\tilde{p}}_c(\theta )$). Experimental data provided by Curcio and published in Curcio et al. (1990).

A number of mechanisms have been hypothesised to explain secondary cone loss, including trophic factor (TF) depletion (Léveillard et al., 2004; Aït-Ali et al., 2015; Mei et al., 2016), oxygen toxicity (Travis et al., 1991; Valter et al., 1998; Stone et al., 1999), metabolic dysregulation (Punzo et al., 2009, 2012), toxic substances (Ripps, 2002), and microglia (Gupta et al., 2003). While not typically related to spatio-temporal patterns of retinal degeneration in the literature, it is reasonable to infer that these mechanisms play an important role in determining spatio-temporal patterns of retinal degeneration.

Grover et al. (1998) have classified the spatio-temporal patterns of visual field loss in RP patients into three patterns and six sub-patterns (see Figure 2). Pattern 1A consists in a restriction of the peripheral visual field, while Pattern 1B also includes a para-/peri-foveal ring scotoma (blind spot); Pattern 2 (A, B and C) involves an initial loss of the superior visual field, winding nasally or temporally into the inferior visual field; lastly, Pattern 3 starts with loss of the mid-peripheral visual field, before spreading into the superior or inferior visual field and winding around the far-periphery. In all cases central vision is the best preserved, though it too is eventually lost (Hamel, 2006; Hartong et al., 2006). Patterns of visual field loss and photoreceptor degeneration (cell loss) are directly related (Escher et al., 2012), loss of the superior visual field corresponding to degeneration of photoreceptors in the inferior retina and vice versa, and loss of the temporal visual field corresponding to degeneration of photoreceptors in the nasal retina and vice versa.

FIGURE 2

Figure 2. Characteristic patterns of visual field loss in human RP (reproduced, with permission, from Roberts et al., 2018). Visual field loss patterns can be classified into three cases and six subcases (classification system described in Grover et al., 1998). Large gray arrows indicate transitions between stages of visual field loss and small red arrows indicate the direction of scotoma (blind spot) propagation. See text for details.

In this article, we explore the conditions under which the TF mechanism, in isolation, can replicate the patterns of cone degeneration observed in vivo. Isolating a mechanism in this way enables us to identify the effects for which it is sufficient to account, avoiding confusion with other mechanistic causes. Understanding the mechanisms of secondary cone degeneration is important since it is the cones that provide high-acuity color vision, and hence their loss, rather than the preceding rod loss, which is the most debilitating. Therefore, by elucidating these mechanisms, we can develop targeted therapies to prevent or delay cone loss, preserving visual function. The TF mechanism has been studied in detail. Rod photoreceptors have been shown to produce a TF called rod-derived cone viability factor (RdCVF), which is necessary for cone survival (Mohand-Saïd et al., 1997, 1998, 2000; Fintz et al., 2003; Léveillard et al., 2004; Yang et al., 2009). RdCVF increases cone glucose uptake, and hence aerobic glycolysis, by binding to the cone transmembrane protein Basigin-1, which consequently binds to the glucose transporter GLUT1 (Aït-Ali et al., 2015). Cones do not produce RdCVF, thus, when rods are lost, RdCVF concentration drops and cone degeneration follows (though it has been suggested that it may ultimately be oxygen toxicity which kills cones; Léveillard and Sahel, 2017).

Thus far, two groups have developed mathematical models operating under the TF hypothesis. Camacho et al. have developed a series of (non-spatial) dynamical systems ordinary differential equation models to describe the role of RdCVF in health and RP (Colón Vélez et al., 2003; Camacho et al., 2010, 2014, 2016a,b,c, 2019, 2020, 2021; Camacho and Wirkus, 2013; Wifvat et al., 2021). In Roberts (2022), we developed the first partial differential equation (PDE) models of the TF mechanism in RP, predicting the spatial spread of retinal degeneration. It was found that, assuming all cones are equally susceptible to RdCVF deprivation and that rods degenerate exponentially with a fixed decay rate, the mechanism is unable to replicate in vivo patterns of retinal degeneration. Previous modeling studies have also considered the oxygen toxicity (Roberts et al., 2017, 2018 and related Roberts et al., 2016b) and toxic substance (Burns et al., 2002) mechanisms, predicting the spatio-temporal patterns of retinal degeneration they would generate. For a review of these and other mathematical models of the retina in health, development and disease see (Roberts et al., 2016a).

In this study, we extend our work in Roberts (2022) by formulating and solving an inverse problem to determine the spatially heterogeneous cone susceptibility to RdCVF deprivation and rod exponential decay rate profiles that are required to qualitatively recapitulate observed patterns of spatio-temporal degeneration in human RP.

2. Materials and Methods

2.1. Model Formulation

We begin by formulating a reaction-diffusion PDE mathematical model (a simplified version of the model presented in Roberts, 2022). Reaction-diffusion PDE models describe the way in which the spatial distribution of cells and chemicals change over time as a result of processes such as movement (diffusion), production, consumption, death, and decay. We pose the model on a spherical geometry to replicate that of the human retina. This geometry is most naturally represented using a spherical polar coordinate system, (r,θ,ϕ), centred in the middle of the vitreous body, where r ≥ 0 (m) is the distance from the origin, 0 ≤ θ ≤ π (rad) is the polar angle and 0 ≤ ϕ < 2π (rad) is the azimuthal angle. To create a more mathematically tractable model, we simplify the geometry by assuming symmetry about the z-axis (directed outward from the origin through the foveal centre), eliminating variation in the azimuthal direction, and effectively depth-average through the retina, assuming that it lies at a single fixed distance, R > 0 (m), from the origin at all eccentricities, θ, leveraging the fact that the retinal width is two orders of magnitude smaller than the eye's radius (Oyster, 1999). Thus, we have reduced the coordinate system to (R,θ), where R is a positive constant parameter and 0 ≤ θ ≤ Θ is an independent variable, which we bound to range between the fovea (at θ = 0 rad) and the ora serrata (at θ = Θ = 1.33 rad; see Figure 1A). We further simplify the model by non-dimensionalising; scaling the dependent and independent variables so that they and the resultant model parameters are dimensionless and hence unitless. This reduces the number of parameters (including eliminating R) and allows us to identify the dominant terms of the governing equations in the ensuing asymptotic analysis. For this reason, there are no units to be stated in Figures 3–10. For the full dimensional model and non-dimensionalisation see Roberts (2022).

FIGURE 3

Figure 3. Initial conditions, ratio of rods to cones and RdCVF reaction rates. (A) initial conditions used in all simulations, consisting of healthy rod and cone profiles and the corresponding RdCVF profiles under Scalings 1 and 2 (the legend applies to (A) only). (B) variation in the healthy rod:cone ratio, ${\tilde{p}}_r(\theta )/{\tilde{p}}_c(\theta )$, with eccentricity. (C) RdCVF production, consumption and decay rates under Scalings 1 and 2 (Equation (1), the legend applies to (C) only). To obtain f_init(θ) in (A,C), Equations (1) and (4) were solved at steady-state using the finite difference method, with 4001 mesh points, where p_r(θ) = p_rinit(θ) and p_c(θ) = p_cinit(θ). Under Scaling 1, α = 7.01 × 10⁴ and β = 1.79 × 10⁶, while under Scaling 2, α = 7.01 × 10² and β = 1.79 × 10⁴. Remaining parameter values as in Table 2.

We proceed directly to the dimensionless model, which consists of a system of PDEs in terms of the dependent variables: TF concentration, f(θ, t), rod photoreceptor density, p_r(θ, t), and cone photoreceptor density, p_c(θ, t); as functions of the independent variables: polar angle, scaled to lie in the range 0 ≤ θ ≤ 1, and time, t > 0 (see Table 1).

TABLE 1

Table 1. Variables employed in the non-dimensional mathematical model [Equations (1–5)].

The TF equation is as follows

$\begin{array}{ll}\frac{\partial f}{\partial t}={\displaystyle \underset{\text{diffusion}}{{\displaystyle \underbrace{\frac{D_f}{sin(\Theta \theta )}\frac{\partial }{\partial \theta }\left(sin(\Theta \theta )\frac{\partial f}{\partial \theta }\right)}}}}+{\displaystyle \underset{\text{production}}{{\displaystyle \underbrace{\alpha p_r}}}}-{\displaystyle \underset{\text{consumption}}{{\displaystyle \underbrace{\beta f{p}_c}}}}-{\displaystyle \underset{\text{decay}}{{\displaystyle \underbrace{\eta f}}}}\text{,}& (1)\end{array}$

where ∂f/∂t is the rate of change in TF concentration over time and the parameters, D_f, the TF diffusivity, α, the rate of TF production by rods, β, the rate of TF consumption by cones, and η, the rate of TF decay, are positive constants. Trophic factor is free to diffuse across the retina through the interphotoreceptor matrix (Aït-Ali et al., 2015). We assume, in the absence of experimental evidence to the contrary, that all rods produce TF at an equal and constant rate, independent of the local TF concentration, such that the rate of TF production is directly proportional to the local rod density. Similarly, in the absence of further experimental evidence, we assume that all cones consume TF at an equal and constant rate for a given local TF concentration. Applying the physiological version of the Law of Mass Action, which states that the rate of a reaction is directly proportional to the product of the concentrations/densities of the reactants (Murray, 2002, in this case TF and cones), we assume that TF is consumed by cones at a rate directly proportional to the product of the local TF concentration and the local cone density. Lastly, we assume that TF decays exponentially, decreasing at a rate directly proportional to its local concentration, as has been shown to occur for a range of other proteins in living human cells (Eden et al., 2011).

The rod equation takes the following form

$\begin{array}{ll}\frac{\partial p_r}{\partial t}=-{\displaystyle \underset{\begin{array}{c}\text{cell degeneration}\\ \text{(mutation-induced)}\end{array}}{{\displaystyle \underbrace{{\varphi }_r(\theta )p_r}}}}\text{,}& (2)\end{array}$

where ∂p_r/∂t is the rate of change in rod density over time and we allow the variable ϕ_r(θ), the rate of mutation-induced rod degeneration, to vary spatially (functional forms defined in the Results section), or take a constant positive value, ϕ_r. Rods degenerate due to their expression of a mutant gene (Hamel, 2006; Hartong et al., 2006) and are assumed to do so exponentially, at a rate directly proportional to their local density, consistent with measurements of photoreceptor degeneration kinetics in mouse, rat and canine models of RP (Clarke et al., 2000). Unlike with cones, RdCVF does not serve a protective function for rods (Aït-Ali et al., 2015); therefore, their rate of degeneration is independent of the TF concentration. We note that Equation (2) can be solved to yield $p_r(\theta ,t)=p_{r_{\text{init}}}(\theta )e^{-{\varphi }_r(\theta )t}$ (where p_rinit(θ), the initial value of p_r(θ, t), is defined below), provided there is no delay in onset or interruption of degeneration.

The cone equation is as follows

$\begin{array}{ll}\frac{\partial p_c}{\partial t}=-{\displaystyle \underset{\begin{array}{c}\text{cell degeneration}\\ \text{(TF starvation)}\end{array}}{{\displaystyle \underbrace{p_c{\lambda }_2(f)}}}},& (3)\end{array}$

where ∂p_c/∂t is the rate of change in cone density over time. We define the Heaviside step function, H(·), such that

$\begin{array}{l}H(x):=\{\begin{array}{cc}0& \text{if }x<0\text{,}\\ 1& \text{if }x\ge 0\text{,}\end{array}\end{array}$

the function λ₂(f) is given by

$\begin{array}{l}{\lambda }_2(f)=1-H(f-f_{\text{crit}}(\theta ))\text{,}\end{array}$

where we allow the variable f_crit(θ), the TF threshold concentration, to vary spatially (functional forms defined in the Results section), or take a constant positive value, f_crit. Cone density is assumed to remain constant provided the local TF concentration, f(θ, t), remains in the healthy range at or above the critical threshold, f_crit, while cones are assumed to decay exponentially (due to TF starvation) at a rate directly proportional to their local density if f(θ, t) drops below this threshold, again consistent with Clarke et al. (2000)'s measurements of photoreceptor degeneration kinetics.

Having defined the governing [Equations (1–3)], we close the system by imposing boundary and initial conditions. We apply zero-flux boundary conditions at both ends of the domain,

$\begin{array}{ll}\frac{\partial f}{\partial \theta }(0,t)=0=\frac{\partial f}{\partial \theta }(1,t)\text{,}& (4)\end{array}$

where ∂f/∂θ is the TF concentration gradient in the polar direction, such that there is no net flow of TF into or out of the domain. This is justified by symmetry at θ = 0, while we assume that TF cannot escape from the retina where it terminates at the ora serrata (θ = 1). The healthy rod and cone distributions are given by the following functions

$\begin{array}{l}{\tilde{p}}_r(\theta )=B_3\theta e^{-b_3\theta }\text{,}\\ {\tilde{p}}_c(\theta )=B_1{e}^{-b_1\theta }+B_2{e}^{-b_2\theta }\text{,}\end{array}$

where the values of the positive constants B₁, B₂, B₃, b₁, b₂, and b₃ are found by fitting to the mean of Curcio et al. (1990)'s measurements of healthy human rod and cone distributions along the temporal horizontal meridian using the Trust-Region Reflective algorithm in Matlab's curve fitting toolbox (see Figure 1B). Lastly, we impose the initial conditions

$\begin{array}{l}f(\theta ,0)=f_{\text{init}}(\theta ),p_r(\theta ,0)=p_{r_{\text{init}}}(\theta )={\tilde{p}}_r(\theta )\text{,}\\ p_c(\theta ,0)=p_{c_{\text{init}}}(\theta )={\tilde{p}}_c(\theta )\text{,}& (5)\end{array}$

where f_init(θ) is the steady-state solution to Equations (1) and (4) with p_r = p_rinit(θ) and p_c = p_cinit(θ) (see Figure 3A). Thus, the retina starts in the healthy state in all simulations. See Table 2 for the dimensionless parameter values [see Roberts (2022) for dimensional values and justification of parameter values]. The model presented here simplifies that in Roberts (2022) in the following ways: it does not include treatment, cone outer segment regeneration, or initial patches of rod or cone loss, while mutation-induced rod loss is active for all simulations in this study. The present model also adds two new features to the previous model: allowing the rate of mutation-induced rod degeneration, ϕ_r(θ), and the TF threshold concentration, f_crit(θ), to vary spatially, where before they were constant (or piecewise constant in the high f_crit subcase).

TABLE 2

Table 2. Parameters employed in the non-dimensional mathematical model [Equations (1–5)]. Values are given to three significant figures (radians are dimensionless).

2.2. Numerical Solutions

Numerical (computational) solutions to Equations (1–5) were obtained using the method of lines (as in Roberts, 2022), discretising in space and then integrating in time. The time integration was performed using the Matlab routine ode15s, a variable-step, variable-order solver, designed to solve problems involving multiple timescales such as this (Matlab version R2020a was used here and throughout the paper). We used a relative error tolerance of 10⁻⁶ and an absolute error tolerance of 10⁻¹⁰, with the remaining settings at their default values. The number of spatial mesh points employed varies between simulations, taking values of 26, 51, 101, 401, or 4,001. The upper bound of 4,001 mesh points was chosen such that the distance between mesh points corresponds to the average width of a photoreceptor. In each case the maximum computationally feasible mesh density was employed, all mesh densities being sufficient to achieve accurate results. The initial TF profile, f(θ, 0) = f_init(θ), was calculated by discretising Equations (1) and (4) at steady-state, using a finite difference scheme, and solving the consequent system of nonlinear algebraic equations using the Matlab routine fsolve (which employs a Trust–Region–Dogleg algorithm) with p_r = p_rinit(θ) and p_c = p_cinit(θ).

2.3. Inverse Problem

Our previous modeling study of the TF hypothesis predicted patterns of cone degeneration which failed to match any known patterns in human RP (Roberts, 2022). In that study, we made the simplifying assumption that model parameters are spatially uniform, such that they do not vary with retinal eccentricity. While this is a reasonable assumption in most cases, we have reason to believe that two of the parameters—the rate of mutation-induced rod loss, ϕ_r, and the TF threshold concentration, f_crit—may vary spatially (see below), which could help account for in vivo patterns of retinal degeneration.

Rates of rod degeneration in human RP have not been studied in great detail. Thus far, histopathological examination of human RP retinas has revealed that rod degeneration is most severe in the mid-peripheral retina, with relative sparing of rods in the macula and far-periphery until later in the disease (Milam et al., 1998). It may be that this pattern varies depending upon the mutation involved and between individuals (cf. Huang et al., 2012 for which different spatial patterns of rod function loss occur in patients, all of whom have a mutation in the RPGR gene). The rate of decay of rod photoreceptors has also been shown to vary with retinal eccentricity in mouse and pig models of RP (Carter-Dawson et al., 1978; Li et al., 1998). Further, under healthy conditions, the RdCVF concentration at the centre of the retina (near θ = 0) is much lower (f(θ, t) ~ O(10⁻⁵)) than in the remainder of the retina (where f(θ, t) ~ O(0.1) to O(1), see Figure 3A). Therefore, it is reasonable to assume that central retinal cones are able to cope with lower RdCVF concentrations than those toward the periphery, and hence that f_crit is also heterogeneous. To determine whether these heterogeneities could account for cone degeneration patterns in human RP, we formulate and solve something known as an inverse problem.

In an inverse problem we seek to determine the model input required to attain a known/desired output. In this case, the known output is the target cone degeneration profile, t_degen(θ), while the input is either the rate of mutation-induced rod loss profile, ϕ_r(θ), or the TF threshold concentration profile, f_crit(θ), with corresponding inverses denoted as ϕ_r(θ) = ϕ_rinv(θ) and f_crit(θ) = f_critinv(θ), respectively. When searching for ϕ_rinv(θ), we hold the TF threshold concentration constant at $f_{\text{crit}}(\theta )=f_{\text{crit}}=3\times 10^{-5}$, while, when searching for f_critinv(θ), we hold the rate of mutation-induced rod loss constant at ${\varphi }_r(\theta )={\varphi }_r=7.33\times 10^{-2}$. The constant value of f_crit is chosen to lie just below the minimum steady-state value of f(θ), such that, in the absence of rod loss, cones remain healthy, while the constant value of ϕ_r is chosen to be one hundred times higher than the value that can be inferred from measurements in the healthy human retina (Curcio et al., 1993), placing the timescale of the resultant cone loss on the order of decades, in agreement with in vivo RP progression rates (Hamel, 2006; Hartong et al., 2006).

We consider a range of target cone degeneration profiles, summarized in Figure 5 and Table 3, which qualitatively replicate visual field loss Patterns 1A, 1B, and 3 seen in vivo (and hence the corresponding in vivo cone degeneration patterns; taking the degeneration of the far-peripheral retina to occur in a radially symmetric manner in Pattern 3—see Figure 2 and Grover et al., 1998). We do not consider patterns of type 2 (to be explored in a future study) as these cannot be replicated by a 1D model (since the radial symmetry, assumed by the 1D model, is broken by variation in the azimuthal/circumferential direction). For each degeneration pattern, we consider a set of sub-patterns to examine how this affects the shape of the inverses, allowing us to confirm that a modest change in the degeneration pattern results in a modest change in the inverses, while exploring both linear/piecewise linear profiles and more biologically realistic nonlinear (quadratic/cubic/exponential) patterns. We also consider a uniform target cone degeneration profile for comparison.

TABLE 3

Table 3. Target cone degeneration profiles, t_degen(θ).

For each pattern, we consider the effect of two (biologically realistic) scalings for the rate of TF production by rods, α, and the rate of TF consumption by cones, β, upon the inverse profiles: (i) Scaling 1—for which α = 7.01 × 10⁴ and β = 1.79 × 10⁶ as in Roberts (2022); and (ii) Scaling 2—for which α = 7.01 × 10² and β = 1.79 × 10⁴. Under Scaling 1, production and consumption of TF dominate over decay (with rate constant η), such that decay has a negligible effect upon the TF profile and model behavior. Under Scaling 2, TF production and consumption occur at a similar rate to decay, such that they balance each other, resulting in a different TF profile and model behavior (see Figures 3A,C). As discussed in Roberts (2022), none of α, β, or η have been measured. The decay rate, η, was chosen to match the measured decay rate of proteins in living human cells (Eden et al., 2011). Under Scaling 1, the consumption rate, β, is chosen such that it dominates over the decay rate (being a factor ϵ⁻¹ = O(10²) larger), while the production rate, α, is chosen to balance consumption (see the Analytical Inverse Section). This is a sensible scaling as it is likely that cones consume RdCVF at a much faster rate than that at which it decays. It is possible, however, that cones consume RdCVF at a similar rate to its decay rate, which is the scenario we consider in Scaling 2; reducing α and β by a factor of 100 (~ ϵ⁻¹) to bring consumption and production into balance with decay (see the Analytical Inverse Section).

We solve the inverse problem both analytically and numerically (computationally), as described in the Analytical Inverse and Numerical Inverse sections below. Analytical approximations are computationally inexpensive and provide deeper insight into model behavior, while numerical solutions, though computationally intensive, are more accurate.

2.3.1. Analytical Inverse

Less mathematically inclined readers may wish to skip over the following derivation and proceed to the resulting Equations (6–11) and surrounding explanatory text. To derive analytical (algebraic) approximations for the inverses, ϕ_rinv(θ) and f_critinv(θ), we perform an asymptotic analysis, seeking the leading order behavior of Equations (1–5). In other words, we are simplifying our equations, making it possible to solve them algebraically (by hand), by only including those terms (corresponding to specific biological processes, e.g., TF production) which dominate the behavior of the solution, where the method known as “asymptotic analysis” enables us to rationally identify these dominant terms. Proceeding as in Roberts (2022) (where we considered a steady-state problem), we choose ϵ = O(10⁻²) and scale the parameters η = ϵ⁻¹η′ and $b_1={ϵ}^{-1}{b}_1^{\prime }$, introducing the new scaling ${\varphi }_r(\theta )=ϵ{\varphi }_r^{\prime }(\theta )$, as we study the time-dependent problem here (where dashes ′ denote scaled variables and parameters). We consider two possible (biologically realistic) scalings on α and β: (i) Scaling 1—for which α = ϵ⁻²α′ and β = ϵ⁻³β′ as in Roberts (2022) (corresponding to α = 7.01 × 10⁴ and β = 1.79 × 10⁶); and (ii) Scaling 2—for which α = ϵ⁻¹α′ and β = ϵ⁻²β′ (corresponding to α = 7.01 × 10² and β = 1.79 × 10⁴). All remaining parameters are assumed to be O(1). We also scale the dependent variable $p_c(\theta ,t)=ϵp_c^{\prime }(\theta ,t)$, and assume f(θ, t) = O(1) and p_r(θ, t) = O(1).

Applying the above scalings and dropping the dashes (working with the scaled versions of the variables and parameters, but omitting the dashes ′ for notational convenience), Equation (2) becomes

$\begin{array}{l}\frac{\partial p_r}{\partial t}=-ϵ{\varphi }_r(\theta )p_r\text{.}\end{array}$

Thus, on this (fast) timescale, the rod density is constant. Since we are interested in the timescale upon which rods degenerate, we scale time as t′ = ϵt such that the decay term enters the dominant balance. Thus, on this slow timescale, after dropping the dashes, we have that

$\begin{array}{l}\frac{\partial p_r}{\partial t}=-{\varphi }_r(\theta )p_r\text{,}\end{array}$

such that, at leading order, $p_{r_0}(\theta ,t)=p_{r_{\text{ini}{\text{t}}_0}}(\theta )e^{-{\varphi }_r(\theta )t}=B_3\theta e^{-b_3\theta }{e}^{{\varphi }_r(\theta )t}$.

We are interested here in the regime in which cones have not yet degenerated, thus, we assume the leading order cone density remains constant at $p_{c_0}(\theta )=p_{c_{\text{ini}{\text{t}}_0}}(\theta )=B_2{e}^{-b_2\theta }$.

Applying Scaling 1 and the slow timescale to Equation (1), we obtain

$\begin{array}{l}ϵ\frac{\partial f}{\partial t}=D_f\frac{{\partial }^{2}f}{\partial {\theta }^2}+D_f\Theta cot(\Theta \theta )\frac{\partial f}{\partial \theta }+{ϵ}^{-2}\alpha p_r-{ϵ}^{-2}\beta p_{c}f-{ϵ}^{-1}\eta f\text{.}\end{array}$

Since the TF dynamics occur on a faster timescale than mutation-induced rod loss, we make a quasi-steady-state approximation (QSSA), assuming that the TF concentration instantaneously takes its steady-state profile, for any given rod density profile, as the rods degenerate (ϵ∂_tf ~ 0). Thus, at leading order, we obtain

$\begin{array}{l}{f}_{0_{\text{QSSA}}}(\theta )=\frac{\alpha p_{r_0}(\theta ,t)}{\beta p_{c_0}(\theta )}\text{.}\end{array}$

Rearranging this expression and assuming that cone degeneration initiates when f_0QSSA(θ) = f_crit(θ), we obtain the cone degeneration time profile,

$\begin{array}{ll}{t}_{\text{degen}}(\theta )=\frac{1}{{\varphi }_r(\theta )}\left(log\left(\frac{\alpha B_3}{\beta B_2{f}_{\text{crit}}(\theta )}\theta \right)-(b_3-b_2)\theta \right)\text{,}& (6)\end{array}$

the inverse mutation-induced rod degeneration rate profile,

$\begin{array}{ll}{\varphi }_{r_{\text{inv}}}(\theta )=\frac{1}{t_{\text{degen}}(\theta )}\left(log\left(\frac{\alpha B_3}{\beta B_2{f}_{\text{crit}}}\theta \right)-(b_3-b_2)\theta \right)\text{,}& (7)\end{array}$

and the inverse TF threshold concentration profile,

$\begin{array}{ll}{f}_{\text{cri}{\text{t}}_{\text{inv}}}(\theta )=\frac{\alpha B_3}{\beta B_2}\theta e^{-((b_3-b_2)\theta +{\varphi }_r{t}_{\text{degen}}(\theta ))}\text{.}& (8)\end{array}$

Alternatively, if we apply Scaling 2 and the slow timescale to Equation (1) we obtain

$\begin{array}{l}ϵ\frac{\partial f}{\partial t}=D_f\frac{{\partial }^{2}f}{\partial {\theta }^2}+D_f\Theta cot(\Theta \theta )\frac{\partial f}{\partial \theta }+{ϵ}^{-1}\alpha p_r-{ϵ}^{-1}\beta p_{c}f-{ϵ}^{-1}\eta f\text{,}\end{array}$

with the TF decay term, ηf, now entering the dominant balance. Applying the QSSA and proceeding as above we find

$\begin{array}{l}{f}_{0_{\text{QSSA}}}(\theta )=\frac{\alpha p_{r_0}(\theta ,t)}{\beta p_{c_0}(\theta )+\eta }\text{,}\end{array}$

with cone degeneration time profile,

$\begin{array}{ll}{t}_{\text{degen}}(\theta )=\frac{1}{{\varphi }_r(\theta )}\left(log\left(\frac{\alpha B_3}{(\beta B_2+\eta e^{b_2\theta })f_{\text{crit}}(\theta )}\theta \right)-(b_3-b_2)\theta \right)\text{,}& (9)\end{array}$

inverse mutation-induced rod degeneration rate profile,

$\begin{array}{ll}{\varphi }_{r_{\text{inv}}}(\theta )=\frac{1}{t_{\text{degen}}(\theta )}\left(log\left(\frac{\alpha B_3}{(\beta B_2+\eta e^{b_2\theta })f_{\text{crit}}}\theta \right)-(b_3-b_2)\theta \right)\text{,}& (10)\end{array}$

and inverse TF threshold concentration profile,

$\begin{array}{ll}{f}_{\text{cri}{\text{t}}_{\text{inv}}}(\theta )=\frac{\alpha B_3}{(\beta B_2+\eta e^{b_2\theta })}\theta e^{-((b_3-b_2)\theta +{\varphi }_r{t}_{\text{degen}}(\theta ))}\text{.}& (11)\end{array}$

These equations reveal how the inverses, ϕ_rinv(θ) and f_critinv(θ), are influenced by our choices for fixed values of f_crit and ϕ_r, respectively. As can be seen from Equations (7) and (10), ϕ_rinv(θ) is inversely and monotonically related to f_crit, such that as f_crit increases, ϕ_rinv(θ) decreases. Similarly, f_critinv(θ) and ϕ_r are inversely and monotonically related in Equations (8) and (11), such that as ϕ_r increases, f_critinv(θ) decreases. Lastly, as would be expected intuitively, t_degen(θ), ϕ_rinv(θ) and f_critinv(θ) all increase monotonically with increasing TF production, α, and decrease monotonically with increasing TF consumption, β, and TF decay η [Equations (6–8) and (9–11)].

2.3.2. Numerical Inverse

The numerical inverse is calculated by repeatedly solving the forward problem [Equations (1–5)] for different values of the input (ϕ_r(θ) or f_crit(θ)), with the aim of converging upon the inverse (ϕ_rinv(θ) or f_critinv(θ)). To find ϕ_rinv(θ), we use the Matlab routine fminsearch (which uses a simplex search method), while to obtain f_critinv(θ) the Matlab routine patternsearch (which uses an adaptive mesh technique) was found to be more effective. In both cases, the objective function (the quantity we are seeking to minimise) was taken as the sum of squares of the difference between the target cone degeneration profile, t_degen(θ), and the contour described by $p_c(\theta ,t)/{\tilde{p}}_c(\theta )=0.99$ (along which cone degeneration is deemed to have initiated). Equations (1–5) were solved at each iteration as described in the Numerical Solutions section. Numerical inverses were calculated only at those locations (eccentricities) where the analytical inverse failed to generate a t_degen(θ) profile matching the target profile, the analytical inverse being assumed to hold at all other eccentricities.

3. Results

We begin by calculating the cone degeneration profiles, t_degen(θ), in the case where both the rate of mutation induced rod degeneration, ϕ_r, and the TF threshold concentration, f_crit, are spatially uniform (or piecewise constant). We set the standard value for ${\varphi }_r=7.33\times 10^{-2}$ and consider the subcases (i) $f_{\text{crit}}=3\times 10^{-5}$ for 0 ≤ θ ≤ 1 (Figure 4A), and (ii) f_crit = 0.3 for θ > 0.13 while $f_{\text{crit}}=3\times 10^{-5}$ for θ ≤ 0.13 (Figure 4B), as were explored in Roberts (2022). These subcases correspond to the situation in which the TF threshold concentration lies beneath the minimum healthy TF value at all retinal locations (i), and the situation in which foveal cones are afforded special protection compared to the rest of the retina, such that they can withstand lower TF concentrations (ii). For notational simplicity, we shall refer to subcase (ii) simply as f_crit = 0.3 in what follows. As with Figures 6–9, we consider both Scaling 1 and Scaling 2 (see Inverse Problem) on the rate of TF production by rods, α, and the rate of TF consumption by cones, β, calculating both analytical and numerical solutions.

FIGURE 4

Figure 4. Cone degeneration profiles. Graphs show the time, t_degen(θ), at which cones degenerate due to RdCVF deprivation, with constant rate of mutation-induced rod degeneration, ${\varphi }_r=7.33\times 10^{-2}$, and constant TF threshold concentrations: $f_{\text{crit}}=3\times 10^{-5}$ (A) and f_crit = 0.3 (B). The solid black and dashed green curves correspond to Scaling 1 (α = 7.01 × 10⁴ and β = 1.79 × 10⁶), while the solid blue and dashed red curves correspond to Scaling 2 (α = 7.01 × 10² and β = 1.79 × 10⁴). The black and blue solid curves are analytical approximations, obtained by plotting Equations (6) and (9), respectively, while the green and red dashed curves are $p_c(\theta ,t)/{\tilde{p}}_c(\theta )=0.99$ contours, obtained by solving Equations (1–5) using the method of lines with 401 mesh points. (A) Simulation spans ~17.7 years in dimensional variables; (B) simulation spans ~2.8 years in dimensional variables. Insets show magnified portions of each graph. Cone degeneration initiates at the fovea (θ = 0) in (A) and at θ = 0.13 in (B), spreading peripherally (rightwards) in both cases. Degeneration occurs earlier in (B) than in (A) and for Scaling 2 than for Scaling 1 (except near the fovea in (A)). Remaining parameter values as in Table 2.

Cone degeneration initiates at the fovea (θ = 0) in Figure 4A and at θ = 0.13 in Figure 4B, spreading peripherally (rightwards) in both cases, while degeneration also initiates at the ora serrata (θ = 1) under Scaling 2 in both Figures 4A,B, spreading centrally. Degeneration occurs earlier in Figure 4B than in Figure 4A and earlier for Scaling 2 than for Scaling 1 (except near the fovea in Figure 4A). Numerical and analytical solutions agree well, only diverging close to the fovea in Figure 4A, where the analytical solution breaks down. None of these patterns of degeneration match those seen in vivo (see Figure 2).

In Figures 6–9, we calculate the ϕ_r(θ) = ϕ_rinv(θ) and f_crit(θ) = f_critinv(θ) profiles required to qualitatively replicate the cone degeneration profiles, t_degen(θ), observed in vivo (Figure 5), by solving the associated inverse problems (see Inverse Problem). As noted in the Inverse Problem section, when searching for ϕ_rinv(θ), we hold the TF threshold concentration constant at $f_{\text{crit}}(\theta )=f_{\text{crit}}=3\times 10^{-5}$, while, when searching for f_critinv(θ), we hold the rate of mutation-induced rod loss constant at ${\varphi }_r(\theta )={\varphi }_r=7.33\times 10^{-2}$. Analytical inverses are plotted across the domain (0 ≤ θ ≤ 1), while numerical inverses are calculated and plotted only at those locations (eccentricities) where the analytical inverse fails to generate a t_degen(θ) profile matching the target profile (as determined by visual inspection, the t_degen(θ) and target profiles being visually indistinguishable outside of these regions).

FIGURE 5

Figure 5. Target cone degeneration profiles. Panels (left) show cone degeneration profiles, t_degen(θ), qualitatively replicating typical spatio-temporal patterns of visual field loss in RP: (A) Uniform, (B) Pattern 1A, (C) Pattern 1B and (D) Pattern 3. Visual field loss patterns directly correspond to cone degeneration patterns in these radially symmetric cases. We seek to replicate these patterns by finding appropriate ϕ_rinv(θ) and f_critinv(θ) profiles in Figures 6–9. Diagrams on the right show the corresponding 2D patterns of visual field loss — white regions: preserved vision, black regions: scotomas (blind spots), and red arrows: direction of scotoma propagation. Parameters: t₀ = 100 (~ 11.0 years), t₁ = 150 (~ 16.6 years), t₂ = 200 ~ 22.1 years, θ_d1 = 0.1 (~ 7.6 degrees), θ_d2 = 0.55 (~ 41.9 degrees), θ_d3 = 0.3 (~ 22.9 degrees), θ_d4 = 0.4 (~ 30.5 degrees) and θ_d5 = 0.6 (~ 45.7 degrees). Cone degeneration profile formulas and parameters are given in Table 3. Remaining parameter values as in Table 2.

In Figure 6, we calculate inverses for a Uniform degeneration profile. While this pattern is not typically observed in humans, we consider this case as a point of comparison with the non-uniform patterns explored in Figures 7–9. Both inverses, ϕ_rinv(θ) and f_critinv(θ), are monotone increasing for Scaling 1, and increase initially for Scaling 2 before reaching a maximum and decreasing toward the ora serrata (at θ = 1). Consequently, Scaling 1 and 2 inverses, while close near the fovea (θ = 0), diverge toward the ora serrata, this effect being more prominent for f_critinv(θ). The inverse profiles have a similar shape to the t_degen(θ) profiles in Figure 4 (see Discussion). Numerical solutions reveal lower values of the inverses near the fovea, where the analytical approximations break down.

FIGURE 6

Figure 6. Inverse mutation-induced rod degeneration rate and TF threshold concentration—Uniform target cone degeneration profile. (A) inverse mutation-induced rod degeneration rate, ϕ_rinv(θ) ($f_{\text{crit}}=3\times 10^{-5}$); (B) inverse TF threshold concentration, f_critinv(θ) (${\varphi }_r=7.33\times 10^{-2}$). The solid black and dashed green curves correspond to Scaling 1 (α = 7.01 × 10⁴ and β = 1.79 × 10⁶), while the solid blue and dashed red curves correspond to Scaling 2 (α = 7.01 × 10² and β = 1.79 × 10⁴). The black and blue solid curves are analytical approximations to the inverses, obtained by plotting Equations (7) and (10), respectively (A), and Equations (8) and (11) respectively (B). The green and red dashed curves are numerical inverses, obtained by using the Matlab routines fminsearch (A) and patternsearch (B) to calculate the ϕ_r and f_crit profiles for which the contour described by $p_c(\theta ,t)/{\tilde{p}}_c(\theta )=0.99$ matches the target cone degeneration profile, t_degen(θ). Equations (1–5) were solved at each iteration using the method of lines, with 101 mesh points. Insets show magnified portions of each graph. Numerical inverses are calculated and plotted only at those locations (eccentricities) where the analytical inverse fails to generate a t_degen(θ) profile matching the target profile. Inverses are monotone increasing for Scaling 1, and increase initially for Scaling 2 before reaching a maximum and decreasing toward the ora serrata (θ = 1). Numerical solutions reveal lower values of the inverses near the fovea (θ = 0) than the analytical approximations suggest. Cone degeneration profile formulas and parameters are given in Table 3. Remaining parameter values as in Table 2.

FIGURE 7

Figure 7. Inverse mutation-induced rod degeneration rate and TF threshold concentration—Pattern 1A target cone degeneration profiles. (A,C,E) inverse mutation-induced rod degeneration rate, ϕ_rinv(θ) ($f_{\text{crit}}=3\times 10^{-5}$); (B,D,F) inverse TF threshold concentration, f_critinv(θ) (${\varphi }_r=7.33\times 10^{-2}$). (A,B) linear target cone degeneration profile, t_degen(θ); (C,D) concave up quadratic t_degen(θ) profile; (E,F) concave down quadratic t_degen(θ) profile. The solid black and dashed green curves correspond to Scaling 1 (α = 7.01 × 10⁴ and β = 1.79 × 10⁶), while the solid blue and dashed red curves correspond to Scaling 2 (α = 7.01 × 10² and β = 1.79 × 10⁴). The black and blue solid curves are analytical approximations to the inverses, obtained by plotting Equations (7) and (10), respectively, (A,C,E) and Equations (8) and (11), respectively, (B,D,F). The green and red dashed curves are numerical inverses, obtained by using the Matlab routines fminsearch (A,C,E), and patternsearch (B,D,F) to calculate the ϕ_r and f_crit profiles for which the contour described by $p_c(\theta ,t)/{\tilde{p}}_c(\theta )=0.99$ matches the target cone degeneration profile, t_degen(θ). Equations (1–5) were solved at each iteration using the method of lines, with 26, 51, or 101 mesh points. Insets show magnified portions of each graph. Numerical inverses are calculated and plotted only at those locations (eccentricities) where the analytical inverse fails to generate a t_degen(θ) profile matching the target profile. Inverses are monotone increasing functions for both scalings in (A, B, E, F), and for Scaling 1 in (C,D) while the inverses increase initially for Scaling 2 before reaching a maximum and decreasing toward the ora serrata (θ = 1) in (C,D). Numerical solutions reveal lower values of the inverses near the fovea (θ = 0) than the analytical approximations suggest. Cone degeneration profile formulas and parameters are given in Table 3. Remaining parameter values as in Table 2.

Inverses for linear (Figures 7A,B), concave up (quadratic) (Figures 7C,D) and concave down (quadratic) (Figures 7E,F) Pattern 1A degeneration profiles are shown in Figure 7. Inverses are monotone increasing functions for both Scalings 1 and 2 in Figures 7A–F and for Scaling 1 in Figures 7C,D, while the inverses increase initially for Scaling 2 before reaching a maximum and decreasing toward the ora serrata in Figures 7C,D. Numerical solutions reveal lower values of the inverses near the fovea, where the analytical approximations break down.

Figure 8 shows inverses for linear (Figures 8A,B), quadratic (Figures 8C,D) and exponential (Figures 8E,F) Pattern 1B degeneration profiles. Inverses resemble vertically flipped versions of the t_degen(θ) profiles in Figure 5C (see Discussion). Numerical solutions reveal lower values of the inverses near the fovea, where the analytical approximations break down, and higher values in some regions away from the fovea in Figures 8A–D. The discontinuities in the linear and quadratic cases are biologically unrealistic, though consistent with the idealised qualitative target cone degeneration patterns in Figure 5C.

FIGURE 8

Figure 8. Inverse mutation-induced rod degeneration rate and TF threshold concentration—Pattern 1B target cone degeneration profiles. (A,C,E) inverse mutation-induced rod degeneration rate, ϕ_rinv(θ) ($f_{\text{crit}}=3\times 10^{-5}$); (B,D,F) inverse TF threshold concentration, f_critinv(θ) (${\varphi }_r=7.33\times 10^{-2}$). (A,B) linear target cone degeneration profile, t_degen(θ); (C,D) quadratic t_degen(θ) profile; (E,F) exponential t_degen(θ) profile. The solid black and dashed green curves correspond to Scaling 1 (α = 7.01 × 10⁴ and β = 1.79 × 10⁶), while the solid blue and dashed red curves correspond to Scaling 2 (α = 7.01 × 10² and β = 1.79 × 10⁴). The black and blue solid curves are analytical approximations to the inverses, obtained by plotting Equations (7) and (10), respectively, (A,C,E), and Equations (8) and (11), respectively, (B,D,F). The green and red dashed curves are numerical inverses, obtained by using the Matlab routines fminsearch (A,C,E), and patternsearch (B,D,F) to calculate the ϕ_r and f_crit profiles for which the contour described by $p_c(\theta ,t)/{\tilde{p}}_c(\theta )=0.99$ matches the target cone degeneration profile, t_degen(θ). Equations (1–5) were solved at each iteration using the method of lines, with 51 or 101 mesh points. Insets show magnified portions of each graph. Numerical inverses are calculated and plotted only at those locations (eccentricities) where the analytical inverse fails to generate a t_degen(θ) profile matching the target profile. Inverses resemble vertically flipped versions of the t_degen(θ) profiles. Numerical solutions reveal lower values of the inverses near the fovea (θ = 0) than the analytical approximations suggest and higher values in some regions away from the fovea in (A–D). Cone degeneration profile formulas and parameters are given in Table 3. Remaining parameter values as in Table 2.

In Figure 9, we calculate inverses for linear 1 (Figures 9A,B), linear 2 (Figures 9C,D), quadratic (Figures 9E,F), and cubic (Figures 9G,H) Pattern 3 degeneration profiles. Inverses resemble vertically flipped versions of the t_degen(θ) profiles in Figure 5D (see Discussion). Numerical solutions reveal lower values of the inverses near the fovea, where the analytical approximations break down, and higher values in some regions away from the fovea in Figures 9C–F,H. Similarly to Figure 8, the discontinuities in the linear 2 and quadratic cases are biologically unrealistic, though consistent with the idealised qualitative target cone degeneration patterns in Figure 5D.

FIGURE 9

Figure 9. Inverse mutation-induced rod degeneration rate and TF threshold concentration—Pattern 3 target cone degeneration profiles. (A,C,E,G) inverse mutation-induced rod degeneration rate, ϕ_rinv(θ) ($f_{\text{crit}}=3\times 10^{-5}$); (B,D,F,H) inverse TF threshold concentration, f_critinv(θ) (${\varphi }_r=7.33\times 10^{-2}$). (A,B) linear 1 target cone degeneration profile, t_degen(θ); (C,D) linear 2 t_degen(θ) profile; (E,F) quadratic t_degen(θ) profile; (G,H) cubic t_degen(θ) profile. The solid black and dashed green curves correspond to Scaling 1 (α = 7.01 × 10⁴ and β = 1.79 × 10⁶), while the solid blue and dashed red curves correspond to Scaling 2 (α = 7.01 × 10² and β = 1.79 × 10⁴). The black and blue solid curves are analytical approximations to the inverses, obtained by plotting Equations (7) and (10) respectively (A,C,E,G), and Equations (8) and (11), respectively, (B,D,F,H). The green and red dashed curves are numerical inverses, obtained by using the Matlab routines fminsearch (A,C,E,G), and patternsearch (B,D,F,H) to calculate the ϕ_r and f_crit profiles for which the contour described by $p_c(\theta ,t)/{\tilde{p}}_c(\theta )=0.99$ matches the target cone degeneration profile, t_degen(θ). Equations (1–5) were solved at each iteration using the method of lines, with 26, 51 or 101 mesh points. Insets show magnified portions of each graph. Numerical inverses are calculated and plotted only at those locations (eccentricities) where the analytical inverse fails to generate a t_degen(θ) profile matching the target profile. Inverses resemble vertically flipped versions of the t_degen(θ) profiles. Numerical solutions reveal lower values of the inverses near the fovea (θ = 0) than the analytical approximations suggest and higher values in some regions away from the fovea in (C–F,H). Cone degeneration profile formulas and parameters are given in Table 3. Remaining parameter values as in Table 2.

Lastly, in Figure 10, we show simulation results of proportional cone loss for analytical and numerical ϕ_rinv(θ) and f_critinv(θ), for Uniform (Scaling 1, Figures 10A–D), concave up Pattern 1A (Scaling 1, Figures 10E–H), linear Pattern 1B (Scaling 2, Figures 10I–L) and quadratic Pattern 3 (Scaling 2, Figures 10M–P) target degeneration profiles. Cone degeneration profiles generally show good agreement with the target t_degen(θ) profiles. There is some divergence from t_degen(θ) for the analytical inverses near the fovea and at discontinuous or nonsmooth portions of t_degen(θ); this is mostly corrected by the numerical inverses. This correction is not perfect near the centre of the fovea, where cones still degenerate earlier than the target profiles. This occurs because it is necessary to replace the Heaviside step function in λ₂(f) [see Equation (3)] with a hyperbolic tanh function to satisfy the smoothness requirements for the numerical solver, with the result that the initiation of cone degeneration is sensitive to the low TF concentrations (f(θ, t) < 10⁻⁴) in that region.

FIGURE 10

Figure 10. Simulations of proportional cone loss for a range of inverse mutation-induced rod degeneration rates and TF threshold concentrations. Plots show the proportion of cones remaining compared to local healthy values, $p_c(\theta ,t)/{\tilde{p}}_c(\theta )$, across space and over time. (A,E,I,M) analytical inverse mutation-induced rod degeneration rate, ϕ_rinv(θ) ($f_{\text{crit}}=3\times 10^{-5}$); (B,F,J,N) numerical ϕ_rinv(θ) ($f_{\text{crit}}=3\times 10^{-5}$); (C,G,K,O) analytical inverse TF threshold concentration, f_critinv(θ) (${\varphi }_r=7.33\times 10^{-2}$); (D,H,L,P) numerical f_critinv(θ) (${\varphi }_r=7.33\times 10^{-2}$). (A–D) Uniform target cone degeneration profile, t_degen(θ), with Scaling 1 (α = 7.01 × 10⁴ and β = 1.79 × 10⁶); (E–H) Pattern 1A quadratic concave up t_degen(θ) profile with Scaling 1; (I–L) Pattern 1B linear t_degen(θ) profile with Scaling 2 (α = 7.01 × 10² and β = 1.79 × 10⁴); (M–P) Pattern 3 quadratic t_degen(θ) profile with Scaling 2. Equations (1–5) were solved using the method of lines, with 26, 51 or 101 mesh points. Analytical and numerical ϕ_rinv(θ) and f_critinv(θ) are as plotted in Figures 5–8. Solid red curves denote the contours along which $p_c(\theta ,t)/{\tilde{p}}_c(\theta )=0.99$, while dashed green curves show the target t_degen(θ) profiles. Cone degeneration profiles generally show good agreement with the target t_degen(θ) profiles. There is some divergence from t_degen(θ) for the analytical inverses near the fovea (θ = 0) and at discontinuous or nonsmooth portions of t_degen(θ); this is mostly corrected by the numerical inverses. Cone degeneration profile formulas and parameters are given in Table 3. Remaining parameter values as in Table 2.

4. Discussion

The spatio-temporal patterns of retinal degeneration observed in human retinitis pigmentosa (RP) are well characterised; however, the mechanistic explanation for these patterns has yet to be conclusively determined. In this paper, we have formulated a one-dimensional (1D) reaction-diffusion partial differential equation (PDE) model (modified from Roberts, 2022) to predict RP progression under the trophic factor (TF) hypothesis. Using this model, we solved inverse problems to determine the rate of mutation-induced rod loss profiles, ϕ_r(θ) = ϕ_rinv(θ), and TF threshold concentration profiles, f_crit(θ) = f_critinv(θ), that would be required to generate spatio-temporal patterns of cone degeneration qualitatively resembling those observed in vivo, were the TF mechanism solely responsible for RP progression. In reality, multiple mechanisms (including oxidative damage and metabolic dysregulation, Travis et al., 1991; Valter et al., 1998; Stone et al., 1999; Punzo et al., 2009, 2012) likely operate in tandem to drive the initiation and propagation of retinal degeneration in RP. By using mathematics to isolate the TF mechanism, in a way that would be impossible to achieve experimentally, we are able to determine the conditions under which the TF mechanism alone would recapitulate known phenotypes. Having identified these conditions, this paves the way for future biomedical and experimental studies to test our predictions.

Other mechanisms may give rise to spatio-temporal patterns of retinal degeneration different from those predicted for the TF mechanism and may do so using fewer assumptions. For example, our previous work on oxygen toxicity in RP demonstrated that this mechanism can replicate visual field loss Pattern 1 (especially 1B) and the late far-peripheral degeneration stage of Pattern 3, without imposing heterogeneities on the rod decay rate or photoreceptor susceptibility to oxygen toxicity (Roberts et al., 2017, 2018). Further, we hypothesise that the toxic substance hypothesis (in which dying rods release a chemical which kills neighbouring photoreceptors) is best able to explain the early mid-peripheral loss of photoreceptors in Patterns 2 and 3, given the high density of rods in this region. In future work, we will explore the toxic substance and other hypotheses, ultimately combining them together in a more comprehensive modeling framework, aimed at explaining and predicting all patterns of retinal degeneration in RP.

Spatially uniform ϕ_r(θ) and f_crit(θ) profiles fail to replicate any of the in vivo patterns of degeneration (Figure 4), showing that heterogenous profiles are required, all else being equal. Throughout this article, we have considered two scalings on the rate of TF production by rods, α, and the rate of TF consumption by cones, β (denoted as Scalings 1 and 2, see the Inverse Problem section for details). Under Scaling 1, the rod:cone ratio (Figure 3B) dominates the model behavior [see Equation (6)], leading to a monotone, central to peripheral pattern of degeneration, while under Scaling 2, the trophic factor decay term, ηf, enters the dominant balance [see Equation (9)], such that degeneration initiates at both the fovea and (later) at the ora serrata, the degenerative fronts meeting in the mid-/far-periphery (Figure 4).

As discussed in the Inverse Problem section, the rate of mutation-induced rod loss, ϕ_r(θ), is known to be spatially heterogeneous in humans with RP (Milam et al., 1998). The ϕ_r(θ) profile predicted for Pattern 3 is consistent with the preferential loss of rods in the mid-peripheral retina noted by Milam et al. (1998) for human RP. A more extensive biomedical investigation is required to characterise quantitatively the diversity of ϕ_r(θ) profiles across individuals and for different mutations. This would make it possible to determine if the ϕ_r(θ) profiles predicted by our model for cone degeneration Patterns 1A and 1B are realised in human RP patients with those cone degeneration patterns. To the best of our knowledge, we are the first to suggest that the intrinsic susceptibility of cones to RdCVF deprivation, characterised in our models by the TF threshold concentration, f_crit(θ), may vary across the retina. Assuming it does vary, what might account for this phenomenon? There is a precedent for special protection being provided to localised parts of the retina. For example, experiments in mice have found that production of basic fibroblast growth factor (bFGF) and glial fibrillary acidic protein (GFAP) is permanently upregulated along the retinal edges, at the ora serrata and optic disc, to protect against elevated stress in these regions (Mervin and Stone, 2002; Stone et al., 2005). Similarly, in the human retina, rods (though not cones) contain bFGF, with a concentration gradient increasing toward the periphery (Li et al., 1997, potentially explaining the relative sparing of rods often observed at the far-periphery). By analogy, we speculate that, in the human retina, cone protective factors may be upregulated at the fovea to compensate for the low RdCVF concentrations in that region, lowering the local value of f_crit(θ). This hypothesis awaits experimental confirmation.

We solved the inverse functions, ϕ_rinv(θ) and f_critinv(θ), both analytically (algebraically) and numerically (computationally). Analytical solutions are approximations; however, they have the advantage of being easier to compute (increasing their utility for biomedical researchers) and provide a more intuitive understanding of model behavior, while numerical solutions are more accurate, though computationally expensive. We calculated the inverses for a range of target cone degeneration profiles, consisting of a Uniform profile and profiles which qualitatively replicate those found in vivo: Pattern 1A, Pattern 1B and Pattern 3 (Pattern 2 being inaccessible to a 1D model; see Figure 5 and Table 3).

The shapes of the inverse functions are determined partly by the rod and cone distributions, ${\tilde{p}}_r(\theta )$ and ${\tilde{p}}_c(\theta )$, and partly by the target cone degeneration profile, t_degen(θ) [see Equations (7,8,10,11)]. As such, in the Uniform case (Figure 6), the Scaling 1 inverse profiles take a similar shape to the rod:cone ratio (Figure 3B), the inverses being lower toward the fovea to compensate for the smaller rod:cone ratio and hence lower supply of TF to each cone. The Scaling 2 inverse profiles follow a similar trend but decrease toward the ora serrata after peaking in the mid-/far-periphery due to the greater influence of the trophic factor decay term under this scaling. Interestingly, the shapes of these inverse profiles bear a striking resemblance to the cone degeneration profiles for spatially uniform ϕ_r(θ) and f_crit(θ) (Figure 4). This is because lower values of the inverses are required to delay degeneration, in those regions where cones would otherwise degenerate earlier, to achieve a uniform degeneration profile. The inverse functions resemble vertically flipped versions of the target cone degeneration profiles for Patterns 1A, 1B and 3 (Figures 7, 8), this being more apparent for Patterns 1B and 3 due to their more distinctive shapes. This makes sense since lower inverse values are required for later degeneration times. Scaling 2 inverses typically lie below Scaling 1 inverses, compensating for the fact that degeneration generally occurs earlier under Scaling 2 than under Scaling 1 for any given ϕ_r(θ) and f_crit(θ).

Analytical inverses give rise to cone degeneration profiles that accurately match the target cone degeneration profiles, except near the fovea (centred at θ = 0, where the validity of the analytical approximation breaks down) and where the target t_degen(θ) profile is nonsmooth or discontinuous (i.e. linear and quadratic Pattern 1B, and linear 1, linear 2 and quadratic Pattern 3; see Figure 10 for examples). Numerical inverses improve accuracy in these regions, consistently taking lower values near the fovea, delaying degeneration where it occurs prematurely under the analytical approximation.

We have assumed throughout this study that at least one of ϕ_r(θ) and f_crit(θ) is spatially uniform. It is possible, however, that both vary spatially. In this case there are no unique inverses; however, if the profile for one of these functions could be measured experimentally, then the inverse problem for the remaining function could be solved as in this paper.

This work could be extended both experimentally and theoretically. Experimental and biomedical studies could measure how the rate of mutation-induced rod loss and TF threshold concentration vary with location in the retina, noting the spatio-temporal pattern of cone degeneration and comparing with the inverse ϕ_rinv(θ) and f_critinv(θ) profiles predicted by our models. Curcio et al. (1993) have previously measured variation in the rate of rod loss in normal (non-RP) human retinas (where rods degenerated most rapidly in the central retina); a similar approach could be taken to quantify the rate of rod loss in human RP retinas. The parameter f_crit is less straightforward to measure. Léveillard et al. (2004) incubated cone-enriched primary cultures from chicken embryos with glutathione S-transferase-RdCVF (GST-RdCVF) fusion proteins, doubling the number of living cells per plate compared with GST alone. If experiments of this type could be repeated for a range of controlled RdCVF concentrations, then the value of f_crit could be identified. Determining the spatial variation of f_crit(θ) in a foveated human-like retina may not be possible presently; however, the recent development of retinal organoids provides promising steps in this direction (O'Hara-Wright and Gonzalez-Cordero, 2020; Fathi et al., 2021). If organoids could be developed with a specialised macular region, mirroring that found in vivo, then the minimum RdCVF concentration required to maintain cones in health could theoretically be tested at a variety of locations. Further, the distribution of RdCVF, predicted in our models, could theoretically be measured in post-mortem human eyes using fluorescent immunohistochemistry, as was done for the protein neuroglobin by Ostojić et al. (2008) and Rajendram and Rao (2007), and perhaps also fluorescent immunocytochemistry as was done for bFGF by Li et al. (1997). In particular, it would be interesting to see if RdCVF concentration varies with retinal eccentricity as starkly as our model predicts, with extremely low levels in the fovea.

In future work, we will extend our mathematical model to two spatial dimensions, accounting for variation in the azimuthal/circumferential dimension (allowing us to capture radially asymmetric aspects of visual field loss Patterns 2 and 3, and to account for azimuthal variation in the rod and cone distributions), and use quantitative target cone degeneration patterns derived from SD-OCT imaging of RP patients (e.g., as in Escher et al., 2012). We will also adapt the model to consider animal retinas for which the photoreceptor distribution has been well characterised (e.g., rats, mice and pigs, Chandler et al., 1999; Gaillard et al., 2009; Ortín-Martínez et al., 2014).

In conclusion, we have formulated and solved a mathematical inverse problem to determine the rate of mutation-induced rod loss and TF threshold concentration profiles required to explain the spatio-temporal patterns of retinal degeneration observed in human RP. Inverse profiles were calculated for a set of qualitatively distinct degeneration patterns, achieving a close match with the target cone degeneration profiles. Predicted inverse profiles await future experimental verification.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

PR: conceptualisation, methodology, software, validation, formal analysis, investigation, data curation, writings–original draft, writings–review and editing, visualisation, and project administration.

Funding

PR was funded by BBSRC (BB/R014817/1).

Conflict of Interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

PR thanks Tom Baden for allowing the time to pursue this research. PR also thanks the reviewers for their helpful and insightful comments.

References

Aït-Ali, N., Fridlich, R., Millet-Puel, G., Clérin, E., Delalande, F., Jaillard, C., et al. (2015). Rod-derived cone viability factor promotes cone survival by stimulating aerobic glycolysis. Cell 161, 817–832. doi: 10.1016/j.cell.2015.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Birtel, J., Gliem, M., Mangold, E., Müller, P. L., Holz, F. G., Neuhaus, C., et al. (2018). Next-generation sequencing identifies unexpected genotype-phenotype correlations in patients with retinitis pigmentosa. PLoS ONE 13, e0207958. doi: 10.1371/journal.pone.0207958

PubMed Abstract | CrossRef Full Text | Google Scholar

Burns, J., Clarke, G., and Lumsden, C. J. (2002). Photoreceptor death: Spatiotemporal patterns arising from one-hit death kinetics and a diffusible cell death factor. Bull. Math. Biol. 64, 1117–1145. doi: 10.1006/bulm.2002.0320

PubMed Abstract | CrossRef Full Text | Google Scholar

Camacho, E. T., Brager, D., Elachouri, G., Korneyeva, T., Millet-Puel, G., Sahel, J. A., et al. (2019). A mathematical analysis of aerobic glycolysis triggered by glucose uptake in cones. Sci. Rep. 9, 4162. doi: 10.1038/s41598-019-39901-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Camacho, E. T., Colón Vélez, M. A., Hernández, D. J., Bernier, U. R., van Laarhoven, J., and Wirkus, S. (2010). A mathematical model for photoreceptor interactions. J. Theor. Biol. 267, 638–646. doi: 10.1016/j.jtbi.2010.09.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Camacho, E. T., Dobreva, A., Larripa, K., Rǎdulescu, A., Schmidt, D., and Trejo, I. (2021). Mathematical Modeling of Retinal Degeneration: Aerobic Glycolysis in a Single Cone, volume 22 of Association for Women in Mathematics Series Cham: Springer International Publishing, 135–178.

Google Scholar

Camacho, E. T., Lenhart, S., Melara, L. A., Villalobos, M. C., and Wirkus, S. (2020). Optimal control with MANF treatment of photoreceptor degeneration. Math. Med. Biol. 37, 1–21. doi: 10.1093/imammb/dqz003

PubMed Abstract | CrossRef Full Text | Google Scholar

Camacho, E. T., Léveillard, T., Sahel, J. A., and Wirkus, S. (2016a). Mathematical model of the role of RdCVF in the coexistence of rods and cones in a healthy eye. Bull. Math. Biol. 78, 1394–1409. doi: 10.1007/s11538-016-0185-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Camacho, E. T., Melara, L. A., Villalobos, M. C., and Wirkus, S. (2014). Optimal control in the treatment of retinitis pigmentosa. Bull. Math. Biol. 76, 292–313. doi: 10.1007/s11538-013-9919-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Camacho, E. T., Punzo, C., and Wirkus, S. A. (2016b). Quantifying the metabolic contribution to photoreceptor death in retinitis pigmentosa via a mathematical model. J. Theor. Biol. 408, 75–87. doi: 10.1016/j.jtbi.2016.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Camacho, E. T., Radulescu, A., and Wirkus, S. (2016c). Bifurcation analysis of a photoreceptor interaction model for retinitis pigmentosa. Commun. Nonlin. Sci. Numer. Simulat. 38, 267–276. doi: 10.1016/j.cnsns.2016.02.030

CrossRef Full Text | Google Scholar

Camacho, E. T., and Wirkus, S. (2013). Tracing the progression of retinitis pigmentosa via photoreceptor interactions. J. Theor. Biol. 317, 105–118. doi: 10.1016/j.jtbi.2012.09.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Carter-Dawson, L. D., LaVail, M. M., and Sidman, R. L. (1978). Differential effect of the rd mutation on rods and cones in the mouse retina. Invest. Ophthalmol. Vis. Sci. 17, 489–498.

PubMed Abstract | Google Scholar

Chandler, M. J., Smith, P. J., Samuelson, D. A., and Mackay, E. O. (1999). Photoreceptor density of the domestic pig retina. Vet. Ophthalmol. 2, 179–184. doi: 10.1046/j.1463-5224.1999.00077.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Clarke, G., Collins, R. A., Leavitt, B. R., Andrews, D. F., Hayden, M. R., Lumsden, C. J., and McInnes, R. R. (2000). A one-hit model of cell death in inherited neuronal degenerations. Nature 406, 195–199. doi: 10.1038/35018098

PubMed Abstract | CrossRef Full Text | Google Scholar

Colón Vélez, M. A., Hernández, D. J., Bernier, U. R., van Laarhoven, J., and Camacho, E. T. (2003). Mathematical Models for Photoreceptor Interactions. Technical Report, Cornell University, Department of Biological Statistics and Computational Biology.

Google Scholar

Coussa, R. G., Basali, D., Maeda, A., DeBenedictis, M., and Traboulsi, E. I. (2019). Sector retinitis pigmentosa: report of ten cases and a review of the literature. Mol. Vis. 25, 869–889.

PubMed Abstract | Google Scholar

Curcio, C. A., Millican, C. L., Allen, K. A., and Kalina, R. E. (1993). Aging of the human photoreceptor mosaic: evidence for selective vulnerability of rods in central retina. Invest. Ophthalmol. Vis. Sci. 34, 3278–3296.

PubMed Abstract | Google Scholar

Curcio, C. A., Sloan, K. R., Kalina, R. E., and Hendrickson, A. E. (1990). Human photoreceptor topography. J. Comp. Neurol. 292, 497–523. doi: 10.1002/cne.902920402

PubMed Abstract | CrossRef Full Text | Google Scholar

Daiger, S. P., Bowne, S. J., and Sullivan, L. S. (2007). Perspective on genes and mutations causing retinitis pigmentosa. Arch. Ophthalmol. 125, 151–158. doi: 10.1001/archopht.125.2.151

PubMed Abstract | CrossRef Full Text | Google Scholar

Eden, E., Geva-Zatorsky, N., Issaeva, I., Cohen, A., Dekel, E., Danon, T., et al. (2011). Proteome half-life dynamics in living human cells. Science 331, 764–768. doi: 10.1126/science.1199784

PubMed Abstract | CrossRef Full Text | Google Scholar

Escher, P., Tran, H. V., Vaclavik, V., Borruat, F. X., Schorderet, D. F., and Munier, F. L. (2012). Double concentric autofluorescence ring in NR2E3-p.G56R-linked autosomal dominant retinitis pigmentosa. Invest. Ophthalmol. Vis. Sci. 53, 4754–4764. doi: 10.1167/iovs.11-8693

PubMed Abstract | CrossRef Full Text | Google Scholar

Fathi, M., Ross, C. T., and Hosseinzadeh, Z. (2021). Functional 3-dimensional retinal organoids: technological progress and existing challenges. Front. Neurosci. 15, 668857. doi: 10.3389/fnins.2021.668857

PubMed Abstract | CrossRef Full Text | Google Scholar

Fintz, A. C., Audo, I., Hicks, D., Mohand-Saïd, S., Léveillard, T., and Sahel, J. (2003). Partial characterization of retina-derived cone neuroprotection in two culture models of photoreceptor degeneration. Invest. Ophthalmol. Vis. Sci. 44, 818–825. doi: 10.1167/iovs.01-1144

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaillard, F., Kuny, S., and Sauvé, Y. (2009). Topographic arrangement of S-cone photoreceptors in the retina of the diurnal nile grass rat (arvicanthis niloticus). Invest. Ophthalmol. Vis. Sci. 50, 5426–5434. doi: 10.1167/iovs.09-3896

PubMed Abstract | CrossRef Full Text | Google Scholar

Ge, Z., Bowles, K., Goetz, K., Scholl, H. P. N., Wang, F., Wang, X., et al. (2015). NGS-based molecular diagnosis of 105 eyegene^® probands with retinitis pigmentosa. Sci. Rep. 5, 18287. doi: 10.1038/srep18287

PubMed Abstract | CrossRef Full Text | Google Scholar

Grover, S., Fishman, G. A., and Brown Jr, J. (1998). Patterns of visual field progression in patients with retinitis pigmentosa. Ophthalmology 105, 1069–1075. doi: 10.1016/S0161-6420(98)96009-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Gupta, N., Brown, K. E., and Milam, A. H. (2003). Activated microglia in human retinitis pigmentosa, late-onset retinal degeneration, and age-related macular degeneration. Exp. Eye Res. 76, 463–471. doi: 10.1016/S0014-4835(02)00332-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Haer-Wigman, L., van Zelst-Stams, W. A. G., Pfundt, R., van den Born, L. I., Klaver, C. C. W., Verheij, J. B. G. M., et al. (2017). Diagnostic exome sequencing in 266 dutch patients with visual impairment. Eur. J. Hum. 25, 591–599. doi: 10.1038/ejhg.2017.9

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamel, C. (2006). Retinitis pigmentosa. Orphanet. J. Rare Dis. 1, 40. doi: 10.1186/1750-1172-1-40

PubMed Abstract | CrossRef Full Text | Google Scholar

Hartong, D. T., Berson, E. L., and Dryja, T. P. (2006). Retinitis pigmentosa. Lancet 368, 1795–1809. doi: 10.1016/S0140-6736(06)69740-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, W. C., Wright, A. F., Roman, A. J., Cideciyan, A. V., Manson, F. D., Gewaily, D. Y., et al. (2012). RPGR-associated retinal degeneration in human X-linked RP and a murine model. Invest. Ophthalmol. Vis. Sci. 53, 5594–5608. doi: 10.1167/iovs.12-10070

PubMed Abstract | CrossRef Full Text | Google Scholar

Léveillard, T., Mohand-Saïd, S., Lorentz, O., Hicks, D., Fintz, A. C., Clérin, E., et al. (2004). Identification and characterization of rod-derived cone viability factor. Nat. Genet. 36, 755–759. doi: 10.1038/ng1386

PubMed Abstract | CrossRef Full Text | Google Scholar

Léveillard, T., and Sahel, J. A. (2017). Metabolic and redox signaling in the retina. Cell. Mol. Life Sci. 74, 3649–3665. doi: 10.1007/s00018-016-2318-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z. Y., Chang, J. H., and Milam, A. H. (1997). A gradient of basic fibroblast growth factor in rod photoreceptors in the normal human retina. Vis. Neurosci. 14, 671–679.

PubMed Abstract | Google Scholar

Li, Z. Y., Wong, F., Chang, J. H., Possin, D. E., Hao, Y., Petters, R. M., and Milam, A. H. (1998). Rhodopsin transgenic pigs as a model for human retinitis pigmentosa. Invest. Ophthalmol. Vis. Sci. 39, 808–819.

PubMed Abstract | Google Scholar

Mei, X., Chaffiol, A., Kole, C., Yang, Y., Millet-Puel, G., Clérin, E., et al. (2016). The thioredoxin encoded by the rod-derived cone viability factor gene protects cone photoreceptors against oxidative stress. Antiox. Redox Signal. 24, 909–923. doi: 10.1089/ars.2015.6509

PubMed Abstract | CrossRef Full Text | Google Scholar

Mervin, K., and Stone, J. (2002). Developmental death of photoreceptors in the C57BL/6J mouse: association with retinal function and self-protection. Exp. Eye Res. 75, 703–713. doi: 10.1006/exer.2002.2063

PubMed Abstract | CrossRef Full Text | Google Scholar

Milam, A. H., Zong, Y. L., and Fariss, R. N. (1998). Histopathology of the human retina in retinitis pigmentosa. Prog. Retin. Eye Res. 17, 175–205.

PubMed Abstract | Google Scholar

Mohand-Saïd, S., Deudon-Combe, A., Hicks, D., Simonutti, M., Forster, V., Fintz, A. C., et al. (1998). Normal retina releases a diffusible factor stimulating cone survival in the retinal degeneration mouse. Proc. Natl. Acad. Sci. 95, 8357–8362.

PubMed Abstract | Google Scholar

Mohand-Saïd, S., Hicks, D., Dreyfus, H., and Sahel, J. A. (2000). Selective transplantation of rods delays cone loss in a retinitis pigmentosa model. Arch. Ophthalmol. 118, 807–811. doi: 10.1001/archopht.118.6.807

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohand-Saïd, S., Hicks, D., Simonutti, M., Tran-Minh, D., Deudon-Combe, A., Dreyfus, H., et al. (1997). Photoreceptor transplants increase host cone survival in the retinal degeneration (rd) mouse. Ophthalmic Res. 29, 290–297.

PubMed Abstract | Google Scholar

Murray, J. D. (2002). Mathematical Biology I: An Introduction, Interdisciplinary Applied Mathematics. 3rd Edn. New York, NY; Berlin; Heidelberg; Barcelona; Hong Kong; London; Milan; Paris; Singapore; Tokyo: Springer.

Google Scholar

Newton, F., and Megaw, R. (2020). Mechanisms of photoreceptor death in retinitis pigmentosa. Genes 11, 1120. doi: 10.3390/genes11101120

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Hara-Wright, M., and Gonzalez-Cordero, A. (2020). Retinal organoids: a window into human retinal development. Development 147, dev189746. doi: 10.1242/dev.189746

PubMed Abstract | CrossRef Full Text | Google Scholar

Ortín-Martínez, A., Nadal-Nicolás, F. M., Jiménez-López, M., Alburquerque-Béjar, J. J., Nieto-López, L., García-Ayuso, D., et al. (2014). Number and distribution of mouse retinal cone photoreceptors: differences between an albino (Swiss) and a pigmented (C57/BL6) strain. PLoS ONE 9, e102392. doi: 10.1371/journal.pone.0102392

PubMed Abstract | CrossRef Full Text | Google Scholar

Ostojić, J., Grozdanić, S. D., Syed, N. A., Hargrove, M. S., Trent, J. T., Kuehn, M. H., et al. (2008). Patterns of distribution of oxygen-binding globins, neuroglobin and cytoglobin in human retina. Arch. Ophthalmol. 126, 1530–1536. doi: 10.1001/archopht.126.11.1530

PubMed Abstract | CrossRef Full Text | Google Scholar

Oyster, C. W. (1999). The Human Eye: Structure and Function. Sunderland, MA: Sinauer Associates Inc.

Google Scholar

Punzo, C., Kornacker, K., and Cepko, C. L. (2009). Stimulation of the insulin/mTOR pathway delays cone death in a mouse model of retinitis pigmentosa. Nat. Neurosci. 12, 44–52. doi: 10.1038/nn.2234

PubMed Abstract | CrossRef Full Text | Google Scholar

Punzo, C., Xiong, W., and Cepko, C. L. (2012). Loss of daylight vision in retinal degeneration: are oxidative stress and metabolic dysregulation to blame? J. Biol. Chem. 287, 1642–1648. doi: 10.1074/jbc.R111.304428

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajendram, R., and Rao, N. A. (2007). Neuroglobin in normal retina and retina from eyes with advanced glaucoma. Br. J. Ophthalmol. 91, 663–666. doi: 10.1136/bjo.2006.093930

PubMed Abstract | CrossRef Full Text | Google Scholar

Ripps, H. (2002). Cell death in retinitis pigmentosa: gap junctions and the ‘bystander’ effect. Exp. Eye Res. 74, 327–336. doi: 10.1006/exer.2002.1155

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts, P. A. (2022). Mathematical models of retinitis pigmentosa: The trophic factor hypothesis. J. Theor. Biol. 534, 110938. doi: 10.1016/j.jtbi.2021.110938

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts, P. A., Gaffney, E. A., Luthert, P. J., Foss, A. J. E., and Byrne, H. M. (2016a). Mathematical and computational models of the retina in health, development and disease. Prog. Retin. Eye. Res. 53, 48–69. doi: 10.1016/j.preteyeres.2016.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts, P. A., Gaffney, E. A., Luthert, P. J., Foss, A. J. E., and Byrne, H. M. (2016b). Retinal oxygen distribution and the role of neuroglobin. J. Math. Biol. 73, 1–38. doi: 10.1007/s00285-015-0931-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts, P. A., Gaffney, E. A., Luthert, P. J., Foss, A. J. E., and Byrne, H. M. (2017). Mathematical models of retinitis pigmentosa: the oxygen toxicity hypothesis. J. Theor. Biol. 425, 53–71. doi: 10.1016/j.jtbi.2017.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts, P. A., Gaffney, E. A., Whiteley, J. P., Luthert, P. J., Foss, A. J. E., and Byrne, H. M. (2018). Predictive mathematical models for the spread and treatment of hyperoxia-induced photoreceptor degeneration in retinitis pigmentosa. Invest. Ophthalmol. Vis. Sci. 59, 1238–1249. doi: 10.1167/iovs.17-23177

PubMed Abstract | CrossRef Full Text | Google Scholar

Stone, J., Maslim, J., Valter-Kocsi, K., Mervin, K., Bowers, F., Chu, Y., et al. (1999). Mechanisms of photoreceptor death and survival in mammalian retina. Prog. Retin. Eye Res. 18, 689–735. doi: 10.1016/S1350-9462(98)00032-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Stone, J., Mervin, K., Walsh, N., Valter, K., Provis, J. M., and Penfold, P. L. (2005). “Photoreceptor stability and degeneration in mammalian retina: lessons from the edge,” in Macular Degeneration, eds P. Penfold, and J. Provis Berlin: Springer, 149–165.

Google Scholar

Travis, G. H., Sutcliffe, J. G., and Bok, D. (1991). The retinal degeneration slow (rds) gene product is a photoreceptor disc membrane-associated glycoprotein. Neuron 6, 61–70.

PubMed Abstract | Google Scholar

Valter, K., Maslim, J., Bowers, F., and Stone, J. (1998). Photoreceptor dystrophy in the RCS rat: roles of oxygen, debris, and bFGF. Invest. Ophthalmol. Vis. Sci. 39, 2427–2442.

PubMed Abstract | Google Scholar

Wifvat, K., Camacho, E. T., Wirkus, S., and Léveillard, T. (2021). The role of RdCVFL in a mathematical model of photoreceptor interactions. J. Theor. Biol. 520, 110642. doi: 10.1016/j.jtbi.2021.110642

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Mohand-Said, S., Danan, A., Simonutti, M., Fontaine, V., Clérin, E., et al. (2009). Functional cone rescue by RdCVF protein in a dominant model of retinitis pigmentosa. Mol. Ther. 17, 787–795. doi: 10.1038/mt.2009.28

PubMed Abstract | CrossRef Full Text | Google Scholar