Optimal control of a two-equation model of radiotherapy

This paper deals with the optimal control of a mathematical model for the evolution of a low-grade glioma (LGG). We will consider a model of the Fischer-Kolmogorov kind for two compartments of tumor cells, using ideas from Galochkina, Bratus and Perez-Garcia [ 10 ] and Perez-Garcia [ 17 ]. The controls are of the form $(t_1, \dots, t_n; d_1, \dots, d_n)$, where $t_i$ is the $i$-th administration time and $d_i$ is the $i$-th applied radiotherapy dose. In the optimal control problem, we try to find controls that maximize, in an admissible class, the first time at which the tumor mass reaches a critical value $M_{*}$. We present an existence result and, also, some numerical experiments (in the previous paper [ 7 ], we have considered and solved a very similar control problem where tumoral cells of only one kind appear).

1. Introduction. The most frequent type of primary brain tumor is the glioma. These tumors remain a challenge for medicine, since they cannot be erradicated with usual therapies and patients diagnosed with gliomas typically die because of the complications related to evolution. Radiation therapy (RT) is one of the main techniques used to treat malignant gliomas.
In the last years, there has been a lot of activity on the mathematical modelling and analysis of the evolution of gliomas; see [6], [13], [17], [18], [19] and [20]). In this paper, we will consider a model of the Fischer-Kolmogorov kind introduced in [10] for the progression of two compartments of tumor cells in a glioma. The first compartment accounts for the relative density u = u(x, t) of cancer cells proliferating with a typical time t p . The second compartment accounts for the relative density of damaged cells v = v(x, t). After the action of radiation, they are not able to repair the DNA damage. Therefore, they die after some time that depends on t p and the amount of mitosis cycles that they are able to complete, denoted k in this paper (for more details, see [21]).
We assume that, in the absence of therapy, the governing model is the following: where Ω ⊂ R N is a bounded open connected set (N = 1, 2 or 3) and, as usual, is the N -dimensional Laplace operator.
In (1), the positive constant D is the diffusion coefficient associated to cellular motility, measured in mm 2 /day; ρ is the proliferation rate, measured in 1/day units (that is, ρ = 1/t p ); it is assumed that the equations have been normalized in such a way that the maximal cell density is equal to 1; the tumor evolves starting from an initial state cell density (c 0 , 0) with 0 ≤ c 0 ≤ 1 and c 0 ≡ 0.
In radiotherapy, the full treatment involves a total radiation dose D tot split in n instantaneous doses d j = d j (x) (measured in Gy units), administered at times t j . Typically, they are equispaced except for the weekends, where the treatment is stopped. The amount of radiation varies depending on the type and stage of the cancer process. Thus, a typical dose for a solid epithelial tumor ranges from 60 to 80 Gy, while lymphomas are treated with 20 to 40 Gy. Preventive doses are typically around 45 to 60 Gy in 1.8-2 Gy fractions (for breast, head, and neck cancers).
We will assume that the t j satisfy We will also assume that the survival fraction, i.e. the fraction of cells not lethally damaged by a dose d j is given by where α t (Gy −1 ) and β t (Gy −2 ) are respectively the linear and quadratic coefficients for tumor cell damage and α t /β t ≈ 10. The main limitation in any radiation protocol is damage to normal tissues. This restriction can be quantitatively expressed as follows: where α h , β h are the normal tissue parameters, with α h /β h ≈ 2. In clinical practice, it is also typical to assume that with d * ≈ 3.4 Gy.
Each stage of the treatment is modeled by an individual system. Accordingly, in each interval (t j , t j+1 ), we will impose the PDEs in (1), Neumann conditions on the boundary and initial data of the form if j = 1, . . . , n. Here, it must be understood that t n+1 = +∞ and u(·, t + j ), v(·, t + j ) (resp. u(·, t − j ), v(·, t − j )) are the cell densities just after (resp. just before) the j-th dose has been applied.
The main goal of this paper is to provide an answer to the following question: Which is the best radiotherapy strategy for (1), subject to the constraints (4) and (5)? In other words, how do we have to choose the t j and the d j satisfying (2), (4) and (5) in order to get an "optimal" evolution of (u, v) subject to the previous PDE system and initial conditions? In general, the notation will be abridged. Thus, we will denote by L 2 the space L 2 (Ω); H 1 will stand for the Hilbert space of the R m -valued (classes of) functions on Ω that belong to L 2 and possess partial derivatives in L 2 ; we will use L 2 (t i , t i+1 ; H 1 ) instead of L 2 (t i , t i+1 ; H 1 (Ω)), etc. The symbol C will denote, as usual, a generic positive constant.
In the previous paper [7], the first and third authors have analyzed (and solved) a similar problem where only one kind of tumor cells is considered.
The plan of the present paper is the following: Section 2 is devoted to present a detailed formulation of the optimal control problem. Then, we will see in Section 3 that the formulation is correct, i.e. the set U ad of admissible controls is non-empty and the objective functional J is welldefined. In Section 4, we prove that there exist optimal controls. To this purpose, we check that U ad is a bounded closed convex set of an appropriate Hilbert space U and J : U ad → R is weakly sequentially lower semi-continuous. Finally, in Section 5, the results of several numerical experiments are presented.
2. The optimal control problem. In order to present the full optimization problem, let us first indicate how we can measure the benefits of a prescribed radiotherapy strategy.
Thus, let us fix the number of radiotherapy doses equal to n, the specific irradiation times t j with 0 < t 1 ≤ . . . ≤ t n ≤T < +∞, the doses d j ∈ L ∞ (Ω) and the initial cell density c 0 = c 0 (x). We will always assume that c 0 is an nontrivial initial cell density, i.e. 0 ≤ c 0 ≤ 1 a.e. and c 0 ≡ 0, c 0 ≡ 1. (7) Let the u j and v j be defined as follows. First, (u 0 , v 0 ) is the solution to the pre-therapy system in Ω. Then, for j = 1, 2, . . . , n − 1, the j-th cell density (u j , v j ) is the solution to the system where S(d j ) is given by (3). This corresponds to the j-th therapy period. Finally, (u n , v n ) is the solution to the post-therapy system As mentioned above, the doses d j have to satisfy the constraints (4) and (5). Now, we can define the global in time tumor cell concentrations (u, v), with with the convention that t 0 = 0 and t n+1 = +∞. Then, the optimization problem is to find specific irradiation times t j and doses d j such that the overall survival time T * (t 1 , . . . , t n ; d 1 , . . . , d n ) attains a maximum. By definition, T * (t 1 , . . . , t n ; d 1 , . . . , d n ) is the smallest time for which the total tumor mass reaches the critical value M * from the left. More precisely, we want to maximize the pay-off function subject to the constraint (t 1 , . . . , t n ; d 1 , . . . , d n ) ∈ U ad , where U ad is the set of admissible controls Obviously, this is equivalent to minimize the functional J : where we have introduced the notation w = (t 1 , . . . , t n ; d 1 , . . . , d n ) and we have denoted by (u, v) the corresponding solution to (8)-(11), i.e. the associated state.
3. Well-posedness of the state equations (8)- (10). Let us introduce the space Y := C 0 ([0, +∞); L 2 ) and let us consider the mapping S : is the solution to (8)- (10). We will prove that S is well-defined. Let c 0 ∈ L ∞ (Ω) satisfy (7) and let w ∈ U ad be given, with The following holds: There exists exactly one solution to (8), with Finally, for any small δ > 0, one also has with estimates in the corresponding spaces.
Proof. The existence and uniqueness of a (weak) solution satisfying (16) and (17) and the estimates 0 ≤ u 0 ≤ 1 are classical. We have v 0 = 0 and the weak maximum principle shows that 0 ≤ u 0 ≤ 1 in Ω × (0, t 1 ), as a consequence of the assumptions on c 0 . The fact that u 0 satisfies (18) is implied by usual parabolic regularity results and the strong maximum principle applied to u 0 and 1 − u 0 . Lemma 3.2. Assume that 1 ≤ j ≤ n − 1 and 0 < t j < t j+1 . There exists exactly one solution to (9) with Finally, for any small δ > 0, one also has with estimates in the corresponding spaces.
From the weak maximum principle applied to u j , v j and 1 − (u j + v j ), we immediately get u j , v j ≥ 0, u j + v j ≤ 1. On the other hand, the strong maximum principle applied to u j and 1 − ( . Hence, the proof is done. for all T < +∞. Furthermore, u n , v n ≥ 0, u n + v n ≤ 1 and there exists M = M (ρ, D, k, T ) > 0 such that Also, for any δ > 0 and any T > t n + δ, one has with estimates in the corresponding spaces. Finally, if we introduce one has W (t) → |Ω| as t → +∞.
Proof. The proof of (21)- (22) is similar to the proof of (19)- (20). Let us check that W satisfies the announced property. Notice that 0 < u n +v n ≤ b n−1 < 1 at t = t n and, consequently, 0 < W (t n ) < |Ω|. Integrating in space the equations in (10) and taking into account (22), we see that in [t n , +∞). Let us consider the logistic ODE The solutions to (24) have the form z K (t) = Ke ρt 1 + K |Ω| e ρt for K ∈ R and the unique solution to (24) satisfying z(t n ) = W (t n ) corresponds to K = K n , where In view of (23), we have that for all t ≥ t n . Let us assume that Then, we can argue as follows. First, note that where we have set by definition Indeed, this is proved taking into account that From (26) and (27), we see that there exist T 1 , T 2 , . . . with T 1 > t n + 1, T +1 > T + 1 and dU dt (T ) ≥ ρα n (1)ε ∀ ≥ 1.
On the other hand, the second order time derivative of U is uniformly bounded in [t n + 1/2, +∞). This can be seen as follows. First, From the equation satisfied by u n , we get On the other hand, multiplying by −∆u n , we also have Combining these estimates, we see that and, in particular, we deduce that ∇u n (· , t) L 2 is uniformly bounded in [t n + 1/2, +∞). A similar estimate can be obtained for ∇v n (· , t) L 2 . Consequently, there exists κ ∈ (0, 1/2) (independent of ) such that From (30), we deduce that for all m ≥ 1, which is an absurd. Therefore, (26) cannot hold. In view of (25), we get that W (t) goes to |Ω| as t → +∞.
This ends the proof. 4. The existence of a solution to the control problem. This section is devoted to prove the existence of a solution to the optimal control problem considered in Section 2: where J is given by (15). The following result holds: Theorem 4.1. Let us assume that 0 < M * < |Ω|. Then there exists at least one solution to the optimal control problem (31).
Let us set S(T ) : Indeed, let us suppose that f m ∈ S(T ) for all m and f m → f weakly in R n × L 2 (Ω) n and let us see that f ∈ S(T ). Let us put f m = (t 1,m , . . . , t n,m ; e 1,m , . . . , e n,m ) and f = (t 1 , . . . , t n ; e 1 , . . . , e n ). By assumption, we have t j,m → t j and e j,m → e j weakly in L 2 (Ω) for all j. Moreover, we can assume that e j,ms → e j weakly- * in L ∞ (Ω).
Several cases are possible: as k → +∞. This is a straightforward consequence of the regularity results depicted in Section 3. Thus, and, consequently, f ∈ S(T ). , and, again, f ∈ S(T ). 3. Finally, if T = t j and all t j,m > T for infinitely many m, one has at least for a subsequence whence we also have f ∈ S(T ).
This argument proves that S(T ) is sequentially weakly closed. Furthermore,J is sequentially l.s.c. for the weak convergence in R n × L 2 (Ω) n . Indeed, let us introduce the indicator χ T , with Then χ T is sequentially weakly l.s.c. Furthermore, in view of the assumptions on M * , for each f ∈Ũ ad there exist times T > 0 (large enough) such that S(T ) f . Now, let us set Then, ϕ * is real-valued. Since it is the least upper bound of a family of sequentially weakly l.s.c. functions, ϕ * is also sequentially weakly l.s.c. But we obviously havẽ J(f ) = ϕ * (f ) for all f ; therefore, the same holds forJ. In view of the well known results of the calculus of variations, there exists at least one global minimizer ofJ inŨ ad . We immediately deduce that an optimal control exists for (31).
This ends the proof.

5.
Numerical experiments: Problem with spherical symmetry. In this section, we present the results of some numerical experiments concerning the solution of a simplified version of (31), where we have assumed that N = 3, Ω and c 0 are spherically symmetric, the administration times t j are prescribed and the applied doses d j are constant in space. The optimization problems are solved computationally with the MATLAB routine fmincon, using two different methods: an interior-point algorithm (IP) and a sequential quadratic programming algorithm (SQP). They are respectively described in [2,3,23] and [8,14]. The computed survival times are similar (slightly better for IP; see Table 1), but the "optimal" dose distributions differ, as shown below.
The total tumor mass is computed with standard integration methods.
For the experiments, we have taken Ω = B(0, R) ⊂ R 3 , with R = 2. The following initial cell concentration was fixed: We have taken a standard dose distribution d st = 1.8 Gy and a maximal admissible dose d * = 3.2 Gy. The total number of doses is n = 30 in the first experiment, n = 40 in the second one and n = 60 in the third one.
In our experiments, the motility and the proliferation parameters are D = 0.01 and ρ = 0.05. The critical tumor size parameter is M * = 0.16; the mitosis cycle parameter is k = 0.9; for the parameters related to therapy, taken from [17], the following values have been taken: α h = 0.1,β h = 0.05 and if n = 30, 15 if n = 40, 24 if n = 60.
Note that, if all d j = d st , we have For the spatial and time approximations of (8)-(10), we use a standard finite difference techniques (of course, after the problem is rewritten in spherical variables).
Taking d j = d st for all j, the survival time is T * = 196 days. By contrast, the computed optimal doses (see Figure 1) give T * = 214 days of survival. In Figure 2, we present the evolution in time of the sizes of the compartments u and v. We see that, during the therapy period, the densities decrease considerably, specially for the compartment u; however, in the post-therapy period, the density for u starts to grow again (slowly at the beginning) while, for v, the density is very close to zero and this tumor compartment is practically erradicated. Finally, 3D views of u and v as functions of |x| and t are displayed in Figure 3. It is an interesting fact that the SQP algorithm provides computed optimal doses that give T * = 213 days of survival (only one day less than before) but with a completely different distribution of doses, as can be seen in Figure 4. Also, note that a maximal constant dose strategy (which is admissible in this case for d max = 2 Gy) leads to a slightly worse T * ; see a comparison in Table 1.

5.2.
Second experiment: 40 doses in 8 weeks. The process is similar, but we assume now that 40 irradiations are applied during 8 weeks. We have kept the same constant values used in Section 5.1, except for E * ; now, E * = 15. Here, we havẽ T = 56. Now, if all d j = d st , the survival time is T * = 238 days; on the other hand, the computed optimal doses (shown in Figure 5) lead to T * = 254 days. Figure 6 shows that the cell densities and tumor sizes behave like in Section 5.1 although, obviously, other numerical values are obtained.   Again, with SQP, the computed optimal doses give T * = 251 days (only three days worse). The distribution of doses is essentially uniform; it is displayed in Figure 7. The constant dose therapy can be applied now with d max = 1.0155 Gy and again provides a slightly worse survival time; see Table 1. 5.3. Third experiment: 60 doses in 15 weeks. Here, we assume that 60 irradiations are applied during 15 weeks in the following way: first, an irradiation period of 6 weeks with 30 doses is performed (exactly as in Section 5.1); then, a rest period of 3 weeks is imposed; finally, a second irradiation period of 6 weeks and 30 doses follows.
As before, all constants except E * have been conserved; in this case, we take E * = 24. Now,T = 105.
If all d j = d st , the survival time is T * = 321 days; but the computed optimal doses distribution, displayed in Figure 8, gives T * = 358 days. Again, Figure 9 is similar to those in Sections 5.1 and 5.2. With a SQP algorithm, the computed optimal doses give T * = 353 days (five days less than with the IP algorithm); the corresponding distribution is shown   in Figure 10. With all d j equal to the maximal admissible constant dose d max = 2 Gy, T * is slightly smaller, as shown in Table 1. We present in Table 1 a comparison of the survival times computed by both algorithms and also the values corresponding to the standard (constant) choice d st and the maximal admissible constant dose d max .
6. Some additional comments. The results in this paper may be viewed as a first step in the solution of control problems oriented to radiotherapy strategy optimization for low-grade gliomas. Several generalizations and extensions of the problems considered here will be analyzed in the next future:  • Numerical experiments without symmetry assumptions.
• The case of unprescribed t j and/or spatially variable d j . In fact, it is much more interesting to assume that the d j are supported by a small set ω ⊂ Ω (irradiation must be performed locally). Any theoretical or numerical result established under this assumption will be very useful.