Abstract
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., ; Coussa et al., ). 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., ).
Figure 1
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., ; 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., ), 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
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., ; 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., ). 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., ; Camacho et al., , , ,,, , , ; Camacho and Wirkus, ; 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., ) 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
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, pr(θ, t), and cone photoreceptor density, pc(θ, 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
| Variable | Description |
|---|---|
| θ | Eccentricity |
| t | Time |
| f(θ, t) | Trophic factor concentration |
| pr(θ, t) | Rod density |
| pc(θ, t) | Cone density |
Variables employed in the non-dimensional mathematical model [Equations (1–5)].
The TF equation is as follows
where ∂f/∂t is the rate of change in TF concentration over time and the parameters, Df, 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., ). 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., ).
The rod equation takes the following form
where ∂pr/∂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., ). Unlike with cones, RdCVF does not serve a protective function for rods (Aït-Ali et al., ); therefore, their rate of degeneration is independent of the TF concentration. We note that Equation (2) can be solved to yield (where prinit(θ), the initial value of pr(θ, t), is defined below), provided there is no delay in onset or interruption of degeneration.
The cone equation is as follows
where ∂pc/∂t is the rate of change in cone density over time. We define the Heaviside step function, H(·), such that
the function λ2(f) is given by
where we allow the variable fcrit(θ), the TF threshold concentration, to vary spatially (functional forms defined in the Results section), or take a constant positive value, fcrit. 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, fcrit, 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. ()'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,
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
where the values of the positive constants B1, B2, B3, b1, b2, and b3 are found by fitting to the mean of Curcio et al. ()'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
where finit(θ) is the steady-state solution to Equations (1) and (4) with pr = prinit(θ) and pc = pcinit(θ) (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, fcrit(θ), to vary spatially, where before they were constant (or piecewise constant in the high fcrit subcase).
Table 2
| Parameter | Description | Value |
|---|---|---|
| Θ | Eccentricity of the ora serrata | 1.33 rad |
| Df | Trophic factor diffusivity | 0.237 |
| α | Rate of trophic factor production by rods | 7.01 × 102 or 7.01 × 104 |
| β | Rate of trophic factor consumption by cones | 1.79 × 104 or 1.79 × 106 |
| η | Rate of trophic factor decay | 1.79 × 102 |
| ϕr | Rate of mutation-induced rod degeneration | 7.33 × 10−2 |
| fcrit | Trophic factor threshold concentration | 3 × 10−5 or 0.3 |
| B1 | Cone profile parameter | 1.56 |
| B2 | Cone profile parameter | 0.158 |
| B3 | Rod profile parameter | 10.6 |
| b1 | Cone profile parameter | 71.8 |
| b2 | Cone profile parameter | 2.67 |
| b3 | Rod profile parameter | 3.06 |
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−6 and an absolute error tolerance of 10−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) = finit(θ), 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 pr = prinit(θ) and pc = pcinit(θ).
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, fcrit—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., ; 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−5)) 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 fcrit 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, tdegen(θ), while the input is either the rate of mutation-induced rod loss profile, ϕr(θ), or the TF threshold concentration profile, fcrit(θ), with corresponding inverses denoted as ϕr(θ) = ϕrinv(θ) and fcrit(θ) = fcritinv(θ), respectively. When searching for ϕrinv(θ), we hold the TF threshold concentration constant at , while, when searching for fcritinv(θ), we hold the rate of mutation-induced rod loss constant at . The constant value of fcrit 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., ), 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
| Degeneration pattern | Sub-pattern | Cone degeneration time (tdegen(θ)) |
|---|---|---|
| Uniform | — | t1 |
| Pattern 1A | Linear | t2 − (t2 − t0)θ |
| Quadratic (concave up) | ||
| Quadratic (concave down) | ||
| Pattern 1B | Linear | |
| Quadratic | ||
| Pattern 3 | Linear 1 | |
| Parameter values* | ||
| C0 = t2 = 200 | ||
Target cone degeneration profiles, tdegen(θ).
We choose θd1 and θd2 such that θd2 = (θd1 + 1)/2, so that θd2 lies halfway between θ = θd1 and θ = 1.
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 × 104 and β = 1.79 × 106 as in Roberts (2022); and (ii) Scaling 2—for which α = 7.01 × 102 and β = 1.79 × 104. 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., ). Under Scaling 1, the consumption rate, β, is chosen such that it dominates over the decay rate (being a factor ϵ−1 = O(102) 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 (~ ϵ−1) 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 fcritinv(θ), 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−2) and scale the parameters η = ϵ−1η′ and , introducing the new scaling , 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 α = ϵ−2α′ and β = ϵ−3β′ as in Roberts (2022) (corresponding to α = 7.01 × 104 and β = 1.79 × 106); and (ii) Scaling 2—for which α = ϵ−1α′ and β = ϵ−2β′ (corresponding to α = 7.01 × 102 and β = 1.79 × 104). All remaining parameters are assumed to be O(1). We also scale the dependent variable , and assume f(θ, t) = O(1) and pr(θ, 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
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
such that, at leading order, .
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 .
Applying Scaling 1 and the slow timescale to Equation (1), we obtain
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
Rearranging this expression and assuming that cone degeneration initiates when f0QSSA(θ) = fcrit(θ), we obtain the cone degeneration time profile,
the inverse mutation-induced rod degeneration rate profile,
and the inverse TF threshold concentration profile,
Alternatively, if we apply Scaling 2 and the slow timescale to Equation (1) we obtain
with the TF decay term, ηf, now entering the dominant balance. Applying the QSSA and proceeding as above we find
with cone degeneration time profile,
inverse mutation-induced rod degeneration rate profile,
and inverse TF threshold concentration profile,
These equations reveal how the inverses, ϕrinv(θ) and fcritinv(θ), are influenced by our choices for fixed values of fcrit and ϕr, respectively. As can be seen from Equations (7) and (10), ϕrinv(θ) is inversely and monotonically related to fcrit, such that as fcrit increases, ϕrinv(θ) decreases. Similarly, fcritinv(θ) and ϕr are inversely and monotonically related in Equations (8) and (11), such that as ϕr increases, fcritinv(θ) decreases. Lastly, as would be expected intuitively, tdegen(θ), ϕrinv(θ) and fcritinv(θ) 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 fcrit(θ)), with the aim of converging upon the inverse (ϕrinv(θ) or fcritinv(θ)). To find ϕrinv(θ), we use the Matlab routine fminsearch (which uses a simplex search method), while to obtain fcritinv(θ) 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, tdegen(θ), and the contour described by (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 tdegen(θ) 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, tdegen(θ), in the case where both the rate of mutation induced rod degeneration, ϕr, and the TF threshold concentration, fcrit, are spatially uniform (or piecewise constant). We set the standard value for and consider the subcases (i) for 0 ≤ θ ≤ 1 (Figure 4A), and (ii) fcrit = 0.3 for θ > 0.13 while 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 fcrit = 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
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 fcrit(θ) = fcritinv(θ) profiles required to qualitatively replicate the cone degeneration profiles, tdegen(θ), 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 , while, when searching for fcritinv(θ), we hold the rate of mutation-induced rod loss constant at . 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 tdegen(θ) profile matching the target profile (as determined by visual inspection, the tdegen(θ) and target profiles being visually indistinguishable outside of these regions).
Figure 5
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 fcritinv(θ), 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 fcritinv(θ). The inverse profiles have a similar shape to the tdegen(θ) 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 7
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 tdegen(θ) 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
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 tdegen(θ) 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
Lastly, in Figure 10, we show simulation results of proportional cone loss for analytical and numerical ϕrinv(θ) and fcritinv(θ), 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 tdegen(θ) profiles. There is some divergence from tdegen(θ) for the analytical inverses near the fovea and at discontinuous or nonsmooth portions of tdegen(θ); 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 λ2(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−4) in that region.
Figure 10
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, fcrit(θ) = fcritinv(θ), 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 fcrit(θ) 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, fcrit(θ), 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 fcrit(θ). This hypothesis awaits experimental confirmation.
We solved the inverse functions, ϕrinv(θ) and fcritinv(θ), 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, and , and partly by the target cone degeneration profile, tdegen(θ) [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 fcrit(θ) (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 fcrit(θ).
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 tdegen(θ) 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 fcrit(θ) 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 fcritinv(θ) profiles predicted by our models. Curcio et al. () 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 fcrit 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 fcrit could be identified. Determining the spatial variation of fcrit(θ) 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., ). 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., ). 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., ; Gaillard et al., ; 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.
Funding
PR was funded by BBSRC (BB/R014817/1).
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.
Statements
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.
Acknowledgments
PR thanks Tom Baden for allowing the time to pursue this research. PR also thanks the reviewers for their helpful and insightful comments.
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.
References
1
Aït-AliN.FridlichR.Millet-PuelG.ClérinE.DelalandeF.JaillardC.et al. (2015). Rod-derived cone viability factor promotes cone survival by stimulating aerobic glycolysis. Cell161, 817–832. 10.1016/j.cell.2015.03.023
2
BirtelJ.GliemM.MangoldE.MüllerP. L.HolzF. G.NeuhausC.et al. (2018). Next-generation sequencing identifies unexpected genotype-phenotype correlations in patients with retinitis pigmentosa. PLoS ONE13, e0207958. 10.1371/journal.pone.0207958
3
BurnsJ.ClarkeG.LumsdenC. J. (2002). Photoreceptor death: Spatiotemporal patterns arising from one-hit death kinetics and a diffusible cell death factor. Bull. Math. Biol.64, 1117–1145. 10.1006/bulm.2002.0320
4
CamachoE. T.BragerD.ElachouriG.KorneyevaT.Millet-PuelG.SahelJ. A.et al. (2019). A mathematical analysis of aerobic glycolysis triggered by glucose uptake in cones. Sci. Rep.9, 4162. 10.1038/s41598-019-39901-z
5
CamachoE. T.Colón VélezM. A.HernándezD. J.BernierU. R.van LaarhovenJ.WirkusS. (2010). A mathematical model for photoreceptor interactions. J. Theor. Biol.267, 638–646. 10.1016/j.jtbi.2010.09.006
6
CamachoE. T.DobrevaA.LarripaK.RǎdulescuA.SchmidtD.TrejoI. (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.
7
CamachoE. T.LenhartS.MelaraL. A.VillalobosM. C.WirkusS. (2020). Optimal control with MANF treatment of photoreceptor degeneration. Math. Med. Biol.37, 1–21. 10.1093/imammb/dqz003
8
CamachoE. T.LéveillardT.SahelJ. A.WirkusS. (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. 10.1007/s11538-016-0185-x
9
CamachoE. T.MelaraL. A.VillalobosM. C.WirkusS. (2014). Optimal control in the treatment of retinitis pigmentosa. Bull. Math. Biol.76, 292–313. 10.1007/s11538-013-9919-1
10
CamachoE. T.PunzoC.WirkusS. A. (2016b). Quantifying the metabolic contribution to photoreceptor death in retinitis pigmentosa via a mathematical model. J. Theor. Biol.408, 75–87. 10.1016/j.jtbi.2016.08.001
11
CamachoE. T.RadulescuA.WirkusS. (2016c). Bifurcation analysis of a photoreceptor interaction model for retinitis pigmentosa. Commun. Nonlin. Sci. Numer. Simulat.38, 267–276. 10.1016/j.cnsns.2016.02.030
12
CamachoE. T.WirkusS. (2013). Tracing the progression of retinitis pigmentosa via photoreceptor interactions. J. Theor. Biol.317, 105–118. 10.1016/j.jtbi.2012.09.034
13
Carter-DawsonL. D.LaVailM. M.SidmanR. L. (1978). Differential effect of the rd mutation on rods and cones in the mouse retina. Invest. Ophthalmol. Vis. Sci.17, 489–498.
14
ChandlerM. J.SmithP. J.SamuelsonD. A.MackayE. O. (1999). Photoreceptor density of the domestic pig retina. Vet. Ophthalmol.2, 179–184. 10.1046/j.1463-5224.1999.00077.x
15
ClarkeG.CollinsR. A.LeavittB. R.AndrewsD. F.HaydenM. R.LumsdenC. J.McInnesR. R. (2000). A one-hit model of cell death in inherited neuronal degenerations. Nature406, 195–199. 10.1038/35018098
16
Colón VélezM. A.HernándezD. J.BernierU. R.van LaarhovenJ.CamachoE. T. (2003). Mathematical Models for Photoreceptor Interactions. Technical Report, Cornell University, Department of Biological Statistics and Computational Biology.
17
CoussaR. G.BasaliD.MaedaA.DeBenedictisM.TraboulsiE. I. (2019). Sector retinitis pigmentosa: report of ten cases and a review of the literature. Mol. Vis.25, 869–889.
18
CurcioC. A.MillicanC. L.AllenK. A.KalinaR. E. (1993). Aging of the human photoreceptor mosaic: evidence for selective vulnerability of rods in central retina. Invest. Ophthalmol. Vis. Sci.34, 3278–3296.
19
CurcioC. A.SloanK. R.KalinaR. E.HendricksonA. E. (1990). Human photoreceptor topography. J. Comp. Neurol.292, 497–523. 10.1002/cne.902920402
20
DaigerS. P.BowneS. J.SullivanL. S. (2007). Perspective on genes and mutations causing retinitis pigmentosa. Arch. Ophthalmol.125, 151–158. 10.1001/archopht.125.2.151
21
EdenE.Geva-ZatorskyN.IssaevaI.CohenA.DekelE.DanonT.et al. (2011). Proteome half-life dynamics in living human cells. Science331, 764–768. 10.1126/science.1199784
22
EscherP.TranH. V.VaclavikV.BorruatF. X.SchorderetD. F.MunierF. L. (2012). Double concentric autofluorescence ring in NR2E3-p.G56R-linked autosomal dominant retinitis pigmentosa. Invest. Ophthalmol. Vis. Sci.53, 4754–4764. 10.1167/iovs.11-8693
23
FathiM.RossC. T.HosseinzadehZ. (2021). Functional 3-dimensional retinal organoids: technological progress and existing challenges. Front. Neurosci.15, 668857. 10.3389/fnins.2021.668857
24
FintzA. C.AudoI.HicksD.Mohand-SaïdS.LéveillardT.SahelJ. (2003). Partial characterization of retina-derived cone neuroprotection in two culture models of photoreceptor degeneration. Invest. Ophthalmol. Vis. Sci.44, 818–825. 10.1167/iovs.01-1144
25
GaillardF.KunyS.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. 10.1167/iovs.09-3896
26
GeZ.BowlesK.GoetzK.SchollH. P. N.WangF.WangX.et al. (2015). NGS-based molecular diagnosis of 105 eyegene® probands with retinitis pigmentosa. Sci. Rep.5, 18287. 10.1038/srep18287
27
GroverS.FishmanG. A.BrownJ.Jr (1998). Patterns of visual field progression in patients with retinitis pigmentosa. Ophthalmology105, 1069–1075. 10.1016/S0161-6420(98)96009-2
28
GuptaN.BrownK. E.MilamA. H. (2003). Activated microglia in human retinitis pigmentosa, late-onset retinal degeneration, and age-related macular degeneration. Exp. Eye Res.76, 463–471. 10.1016/S0014-4835(02)00332-9
29
Haer-WigmanL.van Zelst-StamsW. A. G.PfundtR.van den BornL. I.KlaverC. C. W.VerheijJ. B. G. M.et al. (2017). Diagnostic exome sequencing in 266 dutch patients with visual impairment. Eur. J. Hum.25, 591–599. 10.1038/ejhg.2017.9
30
HamelC. (2006). Retinitis pigmentosa. Orphanet. J. Rare Dis.1, 40. 10.1186/1750-1172-1-40
31
HartongD. T.BersonE. L.DryjaT. P. (2006). Retinitis pigmentosa. Lancet368, 1795–1809. 10.1016/S0140-6736(06)69740-7
32
HuangW. C.WrightA. F.RomanA. J.CideciyanA. V.MansonF. D.GewailyD. Y.et al. (2012). RPGR-associated retinal degeneration in human X-linked RP and a murine model. Invest. Ophthalmol. Vis. Sci.53, 5594–5608. 10.1167/iovs.12-10070
33
LéveillardT.Mohand-SaïdS.LorentzO.HicksD.FintzA. C.ClérinE.et al. (2004). Identification and characterization of rod-derived cone viability factor. Nat. Genet.36, 755–759. 10.1038/ng1386
34
LéveillardT.SahelJ. A. (2017). Metabolic and redox signaling in the retina. Cell. Mol. Life Sci.74, 3649–3665. 10.1007/s00018-016-2318-7
35
LiZ. Y.ChangJ. H.MilamA. H. (1997). A gradient of basic fibroblast growth factor in rod photoreceptors in the normal human retina. Vis. Neurosci.14, 671–679.
36
LiZ. Y.WongF.ChangJ. H.PossinD. E.HaoY.PettersR. M.MilamA. H. (1998). Rhodopsin transgenic pigs as a model for human retinitis pigmentosa. Invest. Ophthalmol. Vis. Sci.39, 808–819.
37
MeiX.ChaffiolA.KoleC.YangY.Millet-PuelG.ClérinE.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. 10.1089/ars.2015.6509
38
MervinK.StoneJ. (2002). Developmental death of photoreceptors in the C57BL/6J mouse: association with retinal function and self-protection. Exp. Eye Res.75, 703–713. 10.1006/exer.2002.2063
39
MilamA. H.ZongY. L.FarissR. N. (1998). Histopathology of the human retina in retinitis pigmentosa. Prog. Retin. Eye Res.17, 175–205.
40
Mohand-SaïdS.Deudon-CombeA.HicksD.SimonuttiM.ForsterV.FintzA. 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.
41
Mohand-SaïdS.HicksD.DreyfusH.SahelJ. A. (2000). Selective transplantation of rods delays cone loss in a retinitis pigmentosa model. Arch. Ophthalmol.118, 807–811. 10.1001/archopht.118.6.807
42
Mohand-SaïdS.HicksD.SimonuttiM.Tran-MinhD.Deudon-CombeA.DreyfusH.et al. (1997). Photoreceptor transplants increase host cone survival in the retinal degeneration (rd) mouse. Ophthalmic Res.29, 290–297.
43
MurrayJ. 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.
44
NewtonF.MegawR. (2020). Mechanisms of photoreceptor death in retinitis pigmentosa. Genes11, 1120. 10.3390/genes11101120
45
O'Hara-WrightM.Gonzalez-CorderoA. (2020). Retinal organoids: a window into human retinal development. Development147, dev189746. 10.1242/dev.189746
46
Ortín-MartínezA.Nadal-NicolásF. M.Jiménez-LópezM.Alburquerque-BéjarJ. J.Nieto-LópezL.García-AyusoD.et al. (2014). Number and distribution of mouse retinal cone photoreceptors: differences between an albino (Swiss) and a pigmented (C57/BL6) strain. PLoS ONE9, e102392. 10.1371/journal.pone.0102392
47
OstojićJ.GrozdanićS. D.SyedN. A.HargroveM. S.TrentJ. T.KuehnM. H.et al. (2008). Patterns of distribution of oxygen-binding globins, neuroglobin and cytoglobin in human retina. Arch. Ophthalmol.126, 1530–1536. 10.1001/archopht.126.11.1530
48
OysterC. W. (1999). The Human Eye: Structure and Function. Sunderland, MA: Sinauer Associates Inc.
49
PunzoC.KornackerK.CepkoC. L. (2009). Stimulation of the insulin/mTOR pathway delays cone death in a mouse model of retinitis pigmentosa. Nat. Neurosci.12, 44–52. 10.1038/nn.2234
50
PunzoC.XiongW.CepkoC. L. (2012). Loss of daylight vision in retinal degeneration: are oxidative stress and metabolic dysregulation to blame?J. Biol. Chem.287, 1642–1648. 10.1074/jbc.R111.304428
51
RajendramR.RaoN. A. (2007). Neuroglobin in normal retina and retina from eyes with advanced glaucoma. Br. J. Ophthalmol.91, 663–666. 10.1136/bjo.2006.093930
52
RippsH. (2002). Cell death in retinitis pigmentosa: gap junctions and the ‘bystander’ effect. Exp. Eye Res.74, 327–336. 10.1006/exer.2002.1155
53
RobertsP. A. (2022). Mathematical models of retinitis pigmentosa: The trophic factor hypothesis. J. Theor. Biol.534, 110938. 10.1016/j.jtbi.2021.110938
54
RobertsP. A.GaffneyE. A.LuthertP. J.FossA. J. E.ByrneH. M. (2016a). Mathematical and computational models of the retina in health, development and disease. Prog. Retin. Eye. Res.53, 48–69. 10.1016/j.preteyeres.2016.04.001
55
RobertsP. A.GaffneyE. A.LuthertP. J.FossA. J. E.ByrneH. M. (2016b). Retinal oxygen distribution and the role of neuroglobin. J. Math. Biol.73, 1–38. 10.1007/s00285-015-0931-y
56
RobertsP. A.GaffneyE. A.LuthertP. J.FossA. J. E.ByrneH. M. (2017). Mathematical models of retinitis pigmentosa: the oxygen toxicity hypothesis. J. Theor. Biol.425, 53–71. 10.1016/j.jtbi.2017.05.006
57
RobertsP. A.GaffneyE. A.WhiteleyJ. P.LuthertP. J.FossA. J. E.ByrneH. 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. 10.1167/iovs.17-23177
58
StoneJ.MaslimJ.Valter-KocsiK.MervinK.BowersF.ChuY.et al. (1999). Mechanisms of photoreceptor death and survival in mammalian retina. Prog. Retin. Eye Res.18, 689–735. 10.1016/S1350-9462(98)00032-9
59
StoneJ.MervinK.WalshN.ValterK.ProvisJ. M.PenfoldP. L. (2005). Photoreceptor stability and degeneration in mammalian retina: lessons from the edge, in Macular Degeneration, eds PenfoldP.ProvisJ.Berlin: Springer, 149–165.
60
TravisG. H.SutcliffeJ. G.BokD. (1991). The retinal degeneration slow (rds) gene product is a photoreceptor disc membrane-associated glycoprotein. Neuron6, 61–70.
61
ValterK.MaslimJ.BowersF.StoneJ. (1998). Photoreceptor dystrophy in the RCS rat: roles of oxygen, debris, and bFGF. Invest. Ophthalmol. Vis. Sci.39, 2427–2442.
62
WifvatK.CamachoE. T.WirkusS.LéveillardT. (2021). The role of RdCVFL in a mathematical model of photoreceptor interactions. J. Theor. Biol.520, 110642. 10.1016/j.jtbi.2021.110642
63
YangY.Mohand-SaidS.DananA.SimonuttiM.FontaineV.ClérinE.et al. (2009). Functional cone rescue by RdCVF protein in a dominant model of retinitis pigmentosa. Mol. Ther.17, 787–795. 10.1038/mt.2009.28
Summary
Keywords
partial differential equations, asymptotic analysis, retina, photoreceptors, rod-derived cone viability factor
Citation
Roberts PA (2022) Inverse Problem Reveals Conditions for Characteristic Retinal Degeneration Patterns in Retinitis Pigmentosa Under the Trophic Factor Hypothesis. Front. Aging Neurosci. 14:765966. doi: 10.3389/fnagi.2022.765966
Received
27 August 2021
Accepted
21 March 2022
Published
02 May 2022
Volume
14 - 2022
Edited by
Adam M. Dubis, University College London, United Kingdom
Reviewed by
Filip Van Den Broeck, Ghent University, Belgium; Moussa A. Zouache, The University of Utah, United States; Robert Linsenmeier, Northwestern University, United States
Updates
Copyright
© 2022 Roberts.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Paul A. Roberts p.a.roberts@univ.oxon.org
This article was submitted to Neurocognitive Aging and Behavior, a section of the journal Frontiers in Aging Neuroscience
Disclaimer
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.