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

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.

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.

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 ϕr=7.33×10-2 and consider the subcases (i) fcrit=3×10-5 for 0 ≤ θ ≤ 1 (Figure 4A), and (ii) f_crit = 0.3 for θ > 0.13 while fcrit=3×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.

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 fcrit(θ)=fcrit=3×10-5, while, when searching for f_critinv(θ), we hold the rate of mutation-induced rod loss constant at ϕr(θ)=ϕr=7.33×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).

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.

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.

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.

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.

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, p~r(θ) and p~c(θ), 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.

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

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

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

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.

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.