Optimal Control with RdCVFL for Degenerating Photoreceptors

Both the rod and cone photoreceptors, along with the retinal pigment epithelium have been experimentally and mathematically shown to work interdependently to maintain vision. Further, the theoredoxin-like rod-derived cone viability factor (RdCVF) and its long form (RdCVFL) have proven to increase photoreceptor survival in experimental results. Aerobic glycolysis is the primary source of energy production for photoreceptors and RdCVF accelerates the intake of glucose into the cones. RdCVFL helps mitigate the negative effects of reactive oxidative species and has shown promise in slowing the death of cones in mouse studies. However, this potential treatment and its effects have never been studied in mathematical models. In this work, we examine an optimal control with the treatment of RdCVFL. We mathematically illustrate the potential this treatment might have for treating degenerative retinal diseases such as retinitis pigmentosa, as well as compare this to the results of an updated control model with RdCVF.


Background of Retinitis Pigmentosa and RdCVFL
Vision impairment is a serious problem that affects at least 2.2 billion people worldwide.Cataracts and glaucoma are the two most common causes of blindness in the world because developing countries rarely have access to the medical treatments needed to stop their progression and vision damage (Steinmetz et al. 2021).In contrast, the biggest contributors to vision loss in developed countries with no cure are due to photoreceptor degeneration, such as in age-related macular degeneration (AMD), untreated retinal detachment, and retinitis pigmentosa (RP) (Steinmetz et al. 2021).While a lot is known about these maladies, they still remain uncurable.The leading cause of inherited retinal degeneration that can cause complete loss of vision is RP, which affects about 1 in 4000 people (Besharse and Bok 2011;Shintani et al. 2009).RP is a heterogenous group of diseases characterized by rod photoreceptor death due to one of many different genetic mutations in the rods followed by the death of otherwise healthy cones.The evolving belief among experimentalists is that a combination of starvation and hyperoxia lead to this secondary cone death (Kanan et al. 2022;Léveillard and Sahel 2017;Petit et al. 2018).Recent experimental work has shown that the rod-derived cone viability factor (RdCVF), which is exclusively produced by the rods, increases glucose uptake by the cones; its long form (RdCVFL), produced by both rods and cones, has been shown to reduce photoreceptor degeneration by protecting against oxidative stress (Aït-Ali et al. 2015;Byrne et al. 2015;Elachouri et al. 2015;Léveillard et al. 2004).Thus, RdCVF and RdCVFL together provide a concrete link between starvation and hyperoxia once the rods have degenerated (Byrne et al. 2015;Elachouri et al. 2015;Léveillard and Sahel 2017).Mathematical work has confirmed the respective roles of RdCVF and RdCVFL as well (Aparicio et al. 2022;Camacho and Wirkus 2013;Camacho et al. 2014Camacho et al. , 2016Camacho et al. , ?, 2019Camacho et al. , 2021;;Dobreva et al. 2022;Wifvat et al. 2021).
Both RdCVF and RdCVFL are encoded in the cDNA Nucleoredoxin-like (Nxnl1) gene (Chalmel et al. 2007;Léveillard et al. 2004).RdCVFL is a member of the thioredoxin family and includes a C-terminal extension that confers enzymatic activities (Brennan et al. 2010;Funato and Miki 2007) whereas RdCVF is truncated and has no enzymatic activity (Chalmel et al. 2007).However RdCVF is necessary to stabilize GLUT1 in the plasma membrane of photoreceptors, by binding to basigin-1 (Bsg1), in order to enhance glucose uptake (Aït-Ali et al. 2015).Thioredoxin is an important regulator of cellular redox homeostasis, which catalyzes the reduction of disulfide bonds and in humans it has been implicated in a wide variety of intracellular and extracellular redox regulations (Matsuo and Kimura-Yoshida 2013).
The rods produce the thioredoxin enzymatic protein RdCVFL as well as the trophic factor RdCVF through a process called alternative splicing.Cones on the other hand, only produce RdCVFL (Aït-Ali et al. 2015;Byrne et al. 2015).Therefore, both cones and rods express the RdCVFL protein and it is also helpful to both types of photoreceptor in protecting them from photo-oxidative damage (such as the oxidation of polyunsaturated fatty acids).In mice studies, it has been shown that RdCVFL reduces the effects of photoreceptor degeneration (Byrne et al. 2015).Just like any other protein from the thioredoxin family, the thiol oxidoreductase activity of RdCVFL depends on how high or reduced its oxidized status is.The thioredoxin reductases step in to reduce RdCVFL whenever it becomes oxidized (Léveillard et al. 2017).Usually produced from glucose in the pentose phosphate pathway, nicotinamide adenine dinucleotide phosphate (NADPH) is the cofactor of these reductases and thus most likely RdCVF provides a reducing power to RdCVFL (Dobreva et al. 2023;Léveillard et al. 2017).With a retinal disease such as RP, if the rods die then RdCVF is lost, which contributes to the secondary wave of cone death.It is speculated that part of the reason for the secondary cone death could be a result of oxidative damage due to the loss of RdCVFL from the rods (Elachouri et al. 2015) and thus it is necessary to study both RdCVF and RdCVFL to enhance our understanding of possible therapies that may prevent secondary cone deaths.However the first step towards investigating the joint effect of RdCVFL and RdCVF is to investigate the effects of each separately.
According to Mei et al. (2016), RdCVF is secreted by rods and protects cones by stimulating aerobic glycolysis when it interacts with glucose transporter 1 (GLUT1) and a complex containing basigin-1.The role of Nxnl1 is studied in cones in this paper.Because RdCVFL is encoded by the Nxnl1 gene, we consider Nxnl1 to be a proxy for RdCVFL and thus the absence of RdCVFL is represented by the Nxnl1-/-mouse (a "knockout" mouse created to not have the Nxnl1 gene (Elachouri et al. 2015)).The studies with mice confirm that damage produced by oxidative stress in cones can be reduced through the administration of RdCVFL to mice that have the genetic abnormality of a deletion of the Nxnl1 gene (Elachouri et al. 2015).Further, RdCVFL is also shown to have cell-autonomous protection since cell viability is reduced when the RdCVFL expression is silenced in a cone-enriched culture (Elachouri et al. 2015).Thus both RdCVF and RdCVFL protect the cones in different ways and are each important and can both be incorporated into a treatment or therapy for RP.
RdCVFL protects rod and cone photoreceptors against photo-oxidative stress and damage as shown in mice studies (Byrne et al. 2015;Elachouri et al. 2015).The rd10 mice are models of RP (e.g., see Phillips et al. 2010) and it is shown that RdCVFL reduces the oxidation of polyunsaturated fatty acids caused by the degeneration of photoreceptors (Byrne et al. 2015).Since RdCVFL includes a thioredoxin fold and RdCVF does not (however, it contains the thioredoxic catalytic site sequence CXXC), RdCVF has no thiol oxidoreductase activity (Léveillard et al. 2017).It is shown in Mei et al. (2016) that cones protect themselves against oxidative stress during aging as they produce RdCVFL.It has further been shown in mouse studies (Elachouri et al. 2015) and previous mathematical work (Wifvat et al. 2021) that RdCVFL can slow the death of cones by mitigating the negative effects of reactive oxidative species (ROS) (Clérin et al. 2011).Knowing how RdCVF and RdCVFL both interact with the cones sheds light on RP and suggests the potential for a double-therapeutic approach involving administration of both RdCVF and RdCVFL.Since potential theraputic effects of RdCVF have been studied previously, in this work, we will study the individual treatment with RdCVFL.
It has been experimentally shown previously by Mei et al. (2016), that RdCVFL causes a rescue effect for the cones.It has also been shown that RdCVFL is not expressed in the retinal pigment epithelium (RPE) of the rd1 mouse (Mei et al. 2016;Léveillard and Sahel 2017, another mouse model of RP) but it is expressed in the photoreceptors of the eye.In lab studies with mice and a control group, RdCVFL was shown to reduce the damage done to cones by as much as 63% (Elachouri et al. 2015).The RdCVFL was subretinally injected after rods had degenerated and then the cone function was recorded.It was found that RdCVFL was able to protect cones from the secondary wave of death after rods had died by as much as 47% (Elachouri et al. 2015).One hypothesis of how this works is that after the rods die off and thus RdCVF is no longer being expressed, there is then less RdCVFL and a higher amount of oxidative stress which causes the cones to die (Léveillard and Sahel 2010).Thus, this protein shows a lot of promise for protecting cones even after rods have died, such as in later stages of RP.
A previous mathematical model looked at the administration of RdCVF and used optimal control theory to obtain the same rescue effect that was observed experimentally (Camacho et al. 2014).To the knowledge of the authors, this was the first use of optimal control utilizing treatment to try to prevent photoreceptor degeneration; a subsequent paper incorporating a different treatment via mesencephalic astrocytederived neurotrophic factor (MANF) again incorporated optimal control to reproduce the observed experimental data (Camacho et al. 2020).Optimal control has been used in a wide range of applications in mechanical engineering, economics, and a range of other disciplines (see, e.g., Athans and Falb 2013;Bryson 1975;Kirk 2004;Pontryagin 1987).In recent decades, optimal control in biological applications has become increasingly utilized (Lenhart and Workman 2007).Some applications to biological problems have included drug administration in chemotherapy, treatment strategies in epidemics, and population harvesting models to name a few (see, e.g., Lenhart and Workman (2007)).
This will be the first mathematical model of photoreceptors that includes optimal control treatment with RdCVFL to counteract photoreceptor degeneration.This information can potentially be used for developing treatments for trials.With optimal control, we maximize the amount of photoreceptors while minimizing the toxicity due to the potential treatments, to find the optimal amount of potential treatment for maintaining vision.Specifically, it is important to understand how the presence of RdCVFL impacts the cones and rods so that if a patient has a disease that affects their vision, potential treatments including RdCVFL can be used to help the patient maintain their night or peripheral vision as well as their central vision.

ODE Model with RdCVFL Control and Derivations
The results presented in this manuscript build from the system of equations modeling the ecosystem of the rod cells, cone cells, and retinal pigmented epithelium (RPE):   et al. 2021).
We now introduce a control function to the model ( 1)-(3) to study the effects that RdCVFL treatment has on its equilibrium solutions.The system of equations modified to include the control u ∈ [0, 1] has effectiveness of treatment parameters σ 1 and σ 2 where the value for σ 1 is 192, which is δ r R n0 rounded to the nearest whole number, and σ 2 is 234, δ c C 0 rounded to the nearest whole number where R n 0 and C 0 are the initial cumulative portion of full length of rod and cone outer segments, respectively.
This system can be defined as the following for the rest of the paper, for convenience, with x = (R n , C, T ).The treatment that administers RdCVFL is incorporated into the system in the second term of ( 4) and ( 5) as u, the control function of time t.This treatment is in the form of a subretinal injection of an adeno-associated virus (AAV) vector.Upon the injection of the AAV, the pumping mechanisms of the RPE cause the gradual and complete absorption of the fluid (Martin et al. 2002).
Standard methods in optimal control do not apply due to the nonlinear manner in which the control appears in the system.We thus linearize the system in order to show existence of a control for the linear system and then compare with the nonlinear results.Here the variables and parameters R n , C, B 1 and B 2 are rearranged, only in the terms that include the control u, to set up the equations to be linearized in the next step: The linearization (first order Taylor's approximation about This system can be defined as the following for the rest of the paper: for convenience, with x = (R n , C, T ).
Further, we will compare the control terms of the nonlinear ( 4)-( 6) and linearized ( 11)-( 13) models: Most of the standard general theorems for the existence of optimal controls such as Filippov (1962), Neustadt (1963) and Fleming and Rishel (1975) require the set of velocities, in this case F([u min , u max ]), to be convex.
In our system (4)-( 6), the controls enter nonlinearly, which is the most accurate biological model; however, we cannot ensure existence of optimal controls.Fortunately, due to the magnitude of the parameters σ i B i for i = 1, 2 and the size of the control u, the linearization F(u) of F(u) about 0 is extremely close to the original system; see Fig. 1.Note that the linearization F(u) about u min +u max 2 would be even closer.However for easier readability, we choose to use the linearization about 0, which is already very close.
Because the rods and cones are an order of magnitude apart in healthy conditions with rods to cones 30:1 in mice and 20:1 in humans, we do not use the Euclidean distance but instead report the errors for rods and cones separately.These distances are less than 10% error and less than 1% error respectively (Table 2).Because the control terms are of the form A B+Cu , the derivatives don't change sign between u min and u max and so the maximum distance is at the endpoints.This shows that the control found for the nonlinear system must be very similar to the optimal control that exists.And, in the results that follow, it is shown that the optimal control for the linearized system is indeed very similar.

Objective Functional
Next the objective functional is defined for both the nonlinear ( 4)-( 6) and linearized ( 11)-( 13) control models.This optimal control problem is formulated to take into account the administration of RdCVFL to the cones because of the fact that the secondary death of cones is a result of being exposed to oxidative stress and the inclusion of RdCVFL protects from that stress.While using the cumulative portion of full length of cone outer segments, to maximize cone longevity the objective functional also takes care to minimize the dosage of RdCVFL to the cones.The control is denoted by the function u : [0, t f ] → [0, 1] and is a representation of the percentage of RdCVFL which is administered to assist in cone survival over the time period [0, t f ].This optimal control problem will be compared to experimental results from Elachouri et al. (2015) and thus the dosage u will be considered as being administered over the course of ten days.The objective functional is defined as where the importance of minimizing u is represented as the weight and the set of controls is defined as and this set is closed.The set U is an interval, and hence it is convex.Note that u min = 0, u max = 0.9 to allow for flexibility in efficiency.

Existence of Unique Global Solutions for the Nonlinear System
Theorem 1 For every fixed measurable control u : [0, t f ] → [u min , u max ] and every initial condition (R n0 , C 0 , T 0 ) in ¯ = {(R, C, T ) : R ≥ 0, C ≥ 0, T ≥ 0}, the system defined in (4)-( 6) has a locally unique and bounded solution defined on [0, t f ].
Proof By hypothesis, the controls are bounded and measurable functions of time t.
For every fixed (R n , C, T ) and fixed control u(•), the right hand side f is a measurable function of time, as a composition of a measurable function and a continuous function in the "good" order (Royden and Fitzpatrick 1988).For each fixed t, the vector field f (•, u(t)) is a rational function of (R n , C, T ) without any singularities in the closed first octant, and hence is locally Lipschitz continuous on this set.Using the boundedness of the controls, one can show analogously to Theorem 6 (see "Appendix A") that solutions for this u stay bounded and are thus defined for all positive times t less than or equal to t f .The calculations are the same, except in (39), the term δ r R n is replaced with u max σ 1 , and δ c C is replaced with u max σ 2 .By the Carathéodory theorem (Coddington and Levinson 1955 Theorem 1.1), there exists a unique solution on the interval [0, t f ] for this control u.

Existence of Unique Global Solutions for the Linearized System
Theorem 2 For every fixed measurable control u : [0, 13) has a locally unique and bounded solution defined on [0, t f ].
Proof By hypothesis, the controls are bounded and measurable functions of time t.For every fixed (R n , C, T ) and fixed control u(•), the right hand side f l is a measurable function of time.For each fixed t, the vector field f l (•, u(t)) is a rational function of (R n , C, T ) without any singularities in the closed first octant, and hence is locally Lipschitz continuous on this set.Using the boundedness of the controls, one can show analogously to Theorem 6 that solutions for this u stay bounded and are thus defined for all positive times t less than or equal to t f .The calculations are the same, except in (39) the term . By the Carathéodory theorem (Coddington and Levinson 1955 Theorem 1.1), there exists a unique solution on the interval [0, t f ] for this control u.

Existence of Optimal Control for the Linearized System
Because most standard general theorems for the existence of optimal controls require the set of velocities to be convex, we only show the existence of optimal controls for the system (11)-( 13) where the controls enter linearly.
We will verify that the optimal control problem ( 11)-( 13) together with the objective functional (17) satisfies the conditions of Theorem 2.2 from Lenhart and Workman (2007) for the basic problem (1.2) (p.7 of the same reference).To keep this manuscript self-contained, we paraphrase the statement of that theorem using our terminology and notation, such as noting that the dynamics and cost functional do not explicitly depend on time.
Theorem 3 (Lenhart and Workman, paraphrased and customized) Let the set of controls be Lebesgue measurable on 0 ≤ t ≤ t 1 with values in R. Suppose that f (x, u) is convex in u, and there exists constants C 4 and C 1 , C 2 , C 3 > 0, and β > 1 such that for all t ∈ [0, t 1 ], and all x, x 1 , u ∈ R Then there exists an optimal control u * minimizing J (u) with J (u * ) finite.
To verify that this theorem applies to our system, first note that with Corollary 1 and Remark 1 in the appendix, we may restrict our considerations to a compact, forward invariant simplex S in the first octant.On this set, the dynamics f (x, u) is affine in the control u and rational in the state x = (R, C, T ) and without singularities in the first quadrant, and thus (i) and (ii) are satisfied.Since f is continuously differentiable it is Lipschitz continuous on the compact set S with Lipschitz constant independent of the control, and thus (iii) is satisfied.The Lagrangian in the functional (17) clearly is of the form (iv) with β = 2 and with finite C 4 since the term (α 1 R n + α 2 C) is linear in the state variables and thus bounded on the domain S of interest.

Optimality Conditions
There are first-order necessary optimality conditions which must be met for the optimal control problem that will be outlined in this section.The optimal control problem, subject to the nonlinear state equations ( 4)-( 6) with objective functional I (u) must meet the following conditions in order for the optimal control u * (t) to be characterized.
The previous equations will be rewritten as the following maximization problem for convenience: where Thus an optimal control u * must be as follows: Then continuing to Pontryagin's Maximum Principle, the Hamiltonian is first defined: where next, λ is defined as the vector of adjoint variables, which are represented as Then we let H ≡ H (x, λ, u) and λ i ≡ λ i (t) with i = 1, . . ., 3 for convenience.

Nonlinear System
Even though an optimal control is not guaranteed to exist, we proceed with the standard technique for finding one and then compare with the linear version.The optimal control that we find given the standard procedure will be referred to as the putative optimal control.Using the nonlinear system (4)-( 6) and the objective function, the Hamiltonian is defined as follows: Theorem 4 Pontryagin's Maximum Principle: If u* and x * = (R * n , C * , T * ) are optimal for problem (22) then there exist absolutely continuous and piecewise differentiable adjoint functions 123 The costate λ satisfies the transversality conditions being For all t ∈ [t 0 , t f ], the optimal control u * (t) maximizes the Hamiltonian, i.e.

Description of the Discretized Variables
In this section an approximate solution to the system of equations ( 4)-( 6), x, is described.The following nodes describe how the time interval [0, t f ] is discretized into N subintervals that are equally spaced with h = t i+1 − t i , for i = 0, . . ., N − 1, and Next, the discretized vector x at time t k is defined as It follows similarly that for k = 0, . . ., N , u k = u(t k ) and for all time steps t k with k = 0, 1, 2, . . .N , all discretized values are contained in the vector u = Next, an approximation to the solution of the system of equations ( 4)-( 6) for the state variable x is calculated given initial conditions x 0 ∈ R 3 with initial u over the time interval [0, t f ], by implementing the fourth-order Runge-Kutta method forward in time with k = 0, . . ., N .The forward Runge-Kutta method is as follows: It is important to note that there are several different ways to approximate the value of u k in calculating k 2 and k 3 .In this algorithm, we replace u k with the average, 1 2 (u k + u k+1 ), because this has been shown to be sufficient (Lenhart and Workman 2007).
Next, the Runge-Kutta method is implemented backward in time to numerically solve the adjoint system and compute the solution λ.The definitions of h and discretization of the time variable t ∈ [0, t f ] remains the same as stated previously.Define the adjoint vector λ k ∈ R 3 , for which the values of the adjoint variables at each discrete time t k are the components, as λ k = (λ 1 (t k ), λ 2 (t k ), λ 3 (t k )) T .The initial iterate is set to λ N = 0 in order to enforce the transversality condition.Then, the backward solve Runge-Kutta method is as follows for k = N , N − 1, . . ., 1:

Nonlinear System
Next, the characterization of the optimal control is derived using Pontryagin's Maximum Principle.By this principle, the Hamiltonian must satisfy (33) at optimality.Rearranging gives where we observe that since all the parameters in (34) are positive, there will be only one intersection of u with the right-hand-side, and this intersection will occur in the first quadrant.However, solving for u, and only taking the unique real root of (33), the analytic solution is too cumbersome to enter into MATLAB, and so it is solved for using a root finding method.Then, u becomes the average of the new value u and the old value for u, for each time step k.Then, a new vector of control variables, u = [u 0 , u 1 , u 2 , . . ., u N ], is formed to test convergence once all the values over the discretized time interval [0, t f ] have been calculated.

Linearized System
Now for the linearized system the characterization of the optimal control is derived using Pontryagin's Maximum Principle.For this system, the Hamiltonian must satisfy It follows that the estimated values of u from the previous iteration of FBSM is defined as u old ∈ R N +1 .See, for example, Lenhart and Workman (2007).
Once the estimated values of u are calculated at each iteration, a stopping criteria must be met for FBSM using relative errors for the control, state variables, and adjoint variables.This implementation is now demonstrated using a given tolerance δ and the values for the control u as an example that covers the other variables listed in general.The relative error δ is now defined in terms of the control u and the 1 vector norm: which can be rewritten as the following in order to take into account the fact that the control variable may have the value 0: Next, β 1 is defined following the work done by Lenhart and Workman (2007): Now, for each current iteration of the FBSM method, the relative errors for each of the state variables can be rewritten.The state variables defined as X = ( X 1 , X 2 , X 3 ) T ∈ R 3×(N +1) .These relative errors which are obtained from the current iteration are given in the matrix: For each previous iteration of FBSM, the estimated state variables are also stored, in a matrix that is constructed in the same fashion as X .This matrix is defined as +1) with the difference being that the ith row is denoted by X old i instead of X i .This ith row is defined as follows for both of these matrices: Next, the adjoint variables, defined as λ i , i = 1, . . ., 3, must be taken into consideration and the process must be repeated for them.The estimated values of adjoint variables are similarly calculated using the current iteration of FBSM, which are then stored in the matrix = ( 1 , 2 , 3 ) T ∈ R 4×(N +1) .Once again, old ∈ R 3×(N +1) is a matrix created using the previously estimated values of FBSM from the last iteration, and it follows that  4)-( 6).Recall that for this case, we do not have proof of existence of optimal control.However, for the sake of comparison and biological relevance, we follow the same methods to find u * , which we call the putative optimal control in this instance.The putative optimal control u * is depicted by the bold line.Observe that the control iterates converge to u * , which is achieved in 12 iterations.b Plot of the control estimates for the linearized model with control (11)-( 13).In this case, we do have a mathematical proof of existence for the optimal control.The optimal control u * is depicted by the bold line.Observe that the control iterates converge to u * , which is achieved in 12 iterations Then, the following is obtained: Next, for the stopping criteria to be met in order to obtain convergence, the following must hold: and if this inequality holds, then convergence has been obtained and the process is halted.If the stopping criteria has not been met, then the process is repeated with an additional run through the FBSM method.

Numerical Results of Control Applied to Rods and Cones with RdCVFL
Our nonlinear model ( 4)-( 6) is based on the experimental results from Elachouri et al. (2015) for the Nxnl1-/-mouse.While we don't have a mathematical proof

Estimated
for existence of the optimal control, this experiment shows a decrease in rods and cones over the course of 10 days, with long term behavior approaching E 2 , where all the rods and cones have died off while the nutrients remain.Table 3 gives the parameter values that fit the data and give E 5 stable, with B 1 , B 2 , δ r and δ c from (1)-( 3) having altered values to reflect the case of mice that do not have RdCVFL and have experienced light damage, where we consider Nxnl1 to be a proxy for RdCVFL (Table 4).The α-values are relative weights compared with each other and, in our case, chosen to match the experimental data of the degeneration of R n and C. The values for α were chosen to be α 1 = 1 and α 2 = 1.2 to show the weights of control for rods and cones, respectively.Recall that the values for σ 1 is 192, which is δ r R n0 rounded to the nearest whole number, and σ 2 is 234, δ c C 0 rounded to the nearest whole number.The lack of RdCVFL is illustrated by the δ r and δ c values being set equal to zero, and the light damage is reflected in the lowered values of B 1 and B 2 .The numerical results show that with an value of 7.0e5 in (21), the numerical results agree with the experimental data; see Fig. 3.This shows that the progression of retinal diseases such as RP may be slowed using treatments that involve RdCVFL.
With the linearized control model ( 11)-( 13), and using the same values for all the parameters and initial conditions, there is a higher saving rate but adjusting to be 1.8E+6, the savings rate is similar to the results of the nonlinear version.for the Nxnl1-/-mouse), without control compared to the results of adding control to the rods and cones.Note that no experimental results are available for treatment with RdCVFL for these initial conditions (so there are no data points at t = 10 as there were in Fig. 3)

Sensitivity Analysis
Next in order to see more clearly how the addition of the control impacts the rest of the model and the other parameters in the model, a sensitivity analysis is performed.
The PRCC method is used in order to do a global sensitivity analysis and see which parameters have the strongest effects, and understand how variation in the parameter values affects the outcomes of our biological system.We first analyze the model without control (1)-(3), the nonlinear RdCVFL control model ( 4)-( 6) and the linearized RdCVFL control model ( 11)-( 13).Using the initial conditions presented in the control results in Fig. 3, R n = 6.4e6,C = 1.8e5,T = 8e4, we see from Fig. 5a-d the difference between the photoreceptor PRCC values without control (a)-(b), and with control (c)-(d).When there is not a control added to the model, the initial values of the state variables for rods and cones are important.However, when there is control added, these initial values are no longer significant.Notice that with or without con-trol, the values B 1 and B 2 are significant.These parameters are correlated to the ROS, which RdCVFL has been experimentally shown to counteract.Similarly, the amount of RdCVFL present also affects shedding rates, which are represented by μ n and μ c , which also show strong significance.Notice also that the quantities are also significant, because the related parameters , κ, μ n and a n are all significant as well.These quantities are important because they have been shown in previous mathematical work (Camacho et al. 2016, ?) to be key ratios that determine the stability of the equilibria.This shows that even in the case where RdCVFL is not present, key factors relating to energy uptake and the metabolic processes which are affected by RdCVFL are still significant.
For the initial values R n = 0.5e6, C = 1.8e5,T = 8e4, which are the initial conditions in the control model shown in Fig. 4, all of the PRCC plots are the same for both the nonlinear model without control (1)-( 3), the nonlinear control model ( 4)-( 6) and the linearized control model ( 11)-( 13) except for four cases which are shown in Fig. 6.Only for the cones, we can see that several parameters are now more significant.For cones in the nonlinear model without control (1)-( 3), the initial value of the rods is now significant, along with a, δ n , μ n , K m and B 1 .Then for the nonlinear control model, the new significant parameters are a, δ n , K m and B 1 .For the linearized control model, when = 7e5, which is the same that's used for the nonlinear control model, the newly significant parameters are a, δ n , μ n and K m .Finally, when we use the value that is needed to get similar results from the linear control model as the nonlinear control model, = 1.8e6, we see the same significant parameters as for the previous value, except now B 1 is also significant as well.Additionally, the PRCC plots were also tested at the initial values R n = 3.6e6, C = 1.8e5,T = 8e4 and the results were the same as for the first initial values R n = 6.4e6,C = 1.8e5,T = 8e4 for both the rods and the cones.

Conclusion
This is the first mathematical optimal control model of a light damaged retina and the potential effects of RdCVFL treatment.The RdCVFL control terms are nonlinear for biological accuracy but this results in the standard general theorems for existence of optimal controls failing to apply.We show the close agreement of the original model with the linearization of the control terms in the model and prove existence of an optimal control for the linearized model in Sect. 2. Because of this close agreement, we proceed under the assumption that a numerical optimal control obtained algorithmically and numerically may accurately represent an optimal control of the nonlinear system.The parameters are fit to real world experimental data for Nxnl1-/-mice, and the results show a savings rate for rods and cones that agrees with experiemental data.Similar results are obtained with the linear system by choosing a different .We conclude that with the presence of RdCVFL added to  3).The parameters are in Table 3 except the four changed parameters in Table 4.The c, d PRCC plots taken at 10 days for rods and cones, respectively, for the -/-mouse under the same conditions except now with nonlinear RdCVFL control u added (4)-( 6).The difference is that the initial population of the state variable is no longer significant.Note that this is for the nonlinear model ( 4)-( 6) with initial conditions R n = 6.4e6,C = 1.8e5,T = 8e4.See Fig. 3 for the control plot.The linearized control model ( 11)-( 13) has the same qualitative results as shown in (c, d) the model of a diseased retina without any naturally occurring RdCVFL, the photoreceptor death is greatly reduced.In fact, it is potentially possible to retain the same amount of vision as a healthy eye, which is reflected in the control model reaching the same population levels as the experimental data of the nondiseased eye.This is an optimistic result for the further study of RdCVFL treatment in human trials.
Further, a global sensitivity and uncertainty analysis was performed to understand what the affect is of each input on the variability of the output of the model.This analysis of the control model was compared to the model without control to understand the differences which occur when control is added to the model.Without control, the PRCC plots show the significant parameters that affect the model output the most.These are the initial values R(0), C(0), the shedding rates and cell metabolism μ n , μ c , as well as other uptake and renewal factors such as a, and κ.The only difference when this analysis is performed on the optimal control model, is that the initial values R(0), C(0) are no longer significant.All models showed similar qualitative results for the listed initial conditions, meaning that it is especially important to look into  13) where c has the same value as the nonlinear model, = 7e5 and d has = 1.8e6.We observe that the qualitative results of the nonlinear and linearized models, b, d respectively, are similar and note that we considered 0.4 to be the value for statistical significance in the PRCC plots these significant parameters.These results show that the potential RdCVFL treatment may still be able to take effect and help save the remaining photoreceptor population, even with a small initial population of photoreceptors.This is important because many individuals who suffer from RP do not attain a diagnosis or treatment until the disease has already progressed and affected their vision somewhat.Thus it is hopeful that for all patients, there can still be some saving effect for what vision is left.
Based on both the optimal control and the sensitivity and uncertainty analysis, RdCVFL is shown to be an important factor for the survival of photoreceptor populations.In retinas that do not have RdCVFL being produced, a treatment involving RdCVFL may be able to help patients retain vision and slow the progression of the disease.Based on the sensitivity and uncertainty results, this type of potential treatment shows promise for helping at any stage of such retinal diseases.Overall, RdCVFL has been shown to be an important factor in helping the retinas maintain photoreceptor survival and may potentially help patients retain vision.which turns into the following after substituting the equation definitions: Then, expanding and rearranging the terms, we get: Note that Then, we have the following inequality: Then, substituting h, n, j and q as defined previously, we get the following: Which reduces to: Finally, substituting back W = R + C + ( h j )T , this turns into: Then, dW dt ≤ −nW + q, we get Therefore, W is bounded along solutions for positive times starting in ¯ .Thus, R n , C and T are also bounded on ¯ for positive times.Since is invariant, and solutions starting in ¯ stay bounded for positive times, the solutions exist for all positive times.
As a direct consequence of (40) we have the following corollary which we will use when arguing the existence of optimal controls.
Corollary 1 Using the notation and parameters as in Theorem 6, the compact simplex is positively forward invariant under the flow of the system (1)-(3).
Proof From equation ( 40) it immediately follows that if W > q n then dW dt < 0.

Fig. 1 (
Fig. 1 (Color figure online) a Left panel: the plots of F(u) and F(u).Observe that these are practically indistinguishable at this scale.Right panel: Zoomed in plots of F(u) and F(u)

Fig. 2
Fig. 2 Plots of the control iterates for the RdCVFL control models with = 7e5.a Plot of the control estimates for the nonlinear model with control (4)-(6).Recall that for this case, we do not have proof of existence of optimal control.However, for the sake of comparison and biological relevance, we follow the same methods to find u * , which we call the putative optimal control in this instance.The putative optimal control u * is depicted by the bold line.Observe that the control iterates converge to u * , which is achieved in 12 iterations.b Plot of the control estimates for the linearized model with control (11)-(13).In this case, we do have a mathematical proof of existence for the optimal control.The optimal control u * is depicted by the bold line.Observe that the control iterates converge to u * , which is achieved in 12 iterations

Fig. 3 (Fig. 4 (
Fig. 3 (Color figure online) The results of decay of both cones and rods, for the nonlinear RdCVFL control model (4)-(6).The initial conditions are R n = 6.4e6,C 1.8e5, T = 8e4, over the course of 10 days (reflected by the experimental results ofElachouri et al. (2015) for the Nxnl1-/-mouse, presented in the top panel), without control compared to the results of adding control to the rods and cones.The data points at time t = 10 for the rods and cones are from experimental results ofElachouri et al. (2015)

Fig. 5 (
Fig. 5 (Color figure online) a, b PRCC plots taken at 10 days for rods and cones, respectively, for the -/-mouse with LD for the model without control (1)-(3).The parameters are in Table3except the four changed parameters in Table4.The c, d PRCC plots taken at 10 days for rods and cones, respectively, for the -/-mouse under the same conditions except now with nonlinear RdCVFL control u added (4)-(6).The difference is that the initial population of the state variable is no longer significant.Note that this is for the nonlinear model (4)-(6) with initial conditions R n = 6.4e6,C = 1.8e5,T = 8e4.See Fig.3for the control plot.The linearized control model (11)-(13) has the same qualitative results as shown in(c, d)

Fig. 6
Fig. 6 PRCC plots taken at 10 days for cones, for the -/-mouse with LD and initial conditions R n = 0.5e6, C = 1.8e5,T = 8e4.The a is for the nonlinear model without control (1)-(3) and b is for the nonlinear control model (4)-(6).The c, d the PRCC plots for the linearized model with control (11)-(13) where c has the same value as the nonlinear model, = 7e5 and d has = 1.8e6.We observe that the qualitative results of the nonlinear and linearized models, b, d respectively, are similar and note that we considered 0.4 to be the value for statistical significance in the PRCC plots

Table 2
Percent error for R n and C for F(u) and F(u)