A distributed control problem for a fractional tumor growth model

In this paper, the authors study the distributed optimal control of a system of three evolutionary equations involving fractional powers of three selfadjoint, monotone, unbounded linear operators having compact resolvents. The system is a generalization of a Cahn-Hilliard type phase field system modeling tumor growth that goes back to Hawkins-Daarud et al. (Int. J. Numer. Math. Biomed. Eng. 28 (2012), 3-24). The aim of the control process, which could be realized by either administering a drug or monitoring the nutrition, is to keep the tumor cell fraction under control while avoiding possible harm for the patient. In contrast to previous studies, in which the occurring unbounded operators governing the diffusional regimes were all given by the Laplacian with zero Neumann boundary conditions, the operators may in our case be different; more generally, we consider systems with fractional powers of the type that were studied in the recent work Adv. Math. Sci. Appl. 28 (2019), 343-375 (see arXiv:1906.10874), by the present authors. In the analysis, the Fr\'echet differentiability of the associated control-to-state operator is shown, by also establishing the existence of solutions to the associated adjoint system, and deriving the first-order necessary conditions of optimality for a cost functional of tracking type.


Introduction
The recent paper [1] investigates the evolutionary system where the equations are understood to hold in Ω, a bounded, connected and smooth domain in R 3 , and in the time interval (0, T). In the above system, A 2ρ , B 2σ , and C 2τ , with r > 0, σ > 0, ρ > 0, denote fractional powers of the self-adjoint, monotone, and unbounded, linear operators A, B and C, respectively, which are supposed to be densely defined in H := L 2 (Ω) and to have compact resolvents. Moreover, α and β are positive real parameters. System (1)-(3) is a generalization of a diffuse interface model for tumor growth. Such models, which are usually established in the framework of the Cahn-Hilliard model originating from the theory of phase transitions, have drawn increasing attention in the past years among mathematicians and applied scientists. We cite here just [2][3][4][5][6][7][8] as a sample of pioneering papers in this direction. In this connection, ϕ stands for an order parameter that should attain its values in the interval [−1, 1], where the values −1 and +1 indicate the healthy cell and tumor cell cases, respectively. The variable S represents the nutrient extra-cellular water concentration, u stands for a source term that acts as a control to monitor the evolution of the tumor cell fraction ϕ, and the nonlinearity P occurring in Equations (1) and (3) is a nonnegative and smooth function modeling a proliferation rate. Finally, µ represents the chemical potential, which acts as the driving thermodynamic force of the evolution and is obtained as the variational derivative with respect to the order parameter ϕ of a suitable free energy functional. In this connection, the nonlinearity f denotes the derivative of a double-well potential F which plays the role of a specific local free energy and yields the main contribution to the total free energy. Important examples for F are the so-called classical regular potential and the logarithmic double-well potential, given by the formulas respectively. In definition (5), the constant c 1 is larger than 1, so that F log is nonconvex. Furthermore, the function P in (1) and (3) is nonnegative and smooth. Finally, the datum u appearing in (3) is given.
In the literature, the diffusional developments in the system have usually been modeled by the Laplacian, that is, the case A 2ρ = B 2σ = C 2τ = −∆, accompanied by zero Neumann boundary conditions, was assumed, where two main classes of models were considered. The first class of models regards the tumor and healthy cells as inertialess fluids; in such models, special fluid effects can be incorporated by postulating a Darcy or Stokes-Brinkman law, see, e.g., the works [7,[9][10][11][12][13][14][15][16][17][18], where we also refer to [19,20]. The other class of models, to which the model considered here belongs, neglects the velocity. Typical contributions in this direction were given in [21][22][23][24][25][26][27], to name just a few.
While the occurrence of more general diffusional regimes of fractional type has been studied for a long time in the mathematical literature, it was only recently (see, e.g., [28][29][30][31][32][33][34][35][36]) that fractional operators have been investigated in the framework of Cahn-Hilliard systems (for phase field systems of Caginalp type, see also [37]), and the only investigations of tumor growth models involving fractional diffusive regimes such as in the system (1)-(3) seem to be the recent papers [1,38] by the present authors. Concerning the motivation and biological meaningfulness of the systems (1)-(3), we point out that, in our approach, the three fractional operators, which may be considerably different one from the other, are used as major dynamics for the tumor growth and for the diffusion processes. These operators may be related and close to fractional Laplacians, but also to other elliptic operators, and can present different orders. Indeed, some components in tumor development, such as immune cells, do have anomalous diffusion dynamics observed by experiments [39], but other components, like chemical potential and nutrient concentration may be ruled by other fractional or non-fractional flows. However, taking all this into consideration, it is the case of emphasizing that fractional operators are more and more utilized in the world of biological applications. In this direction, let us mention some recent and related work: we quote [40] for growth of species and their movements; we quote [41] for the numerical study of a super-diffusive model for a benign brain tumor; we quote [39] for the application to cell movements, in particular for T cells migrating through chronically-infected brain tissue; we quote the monography [42] concerned with fractional diffusion processes in applied sciences; we quote [43] for a hyperbolic-parabolic model of chemotaxis related to tumor angiogenesis; we quote [44] for the study of the dynamics of growth bacteria; we quote [45] for the study of a space time fractional reaction-diffusion equation in Parkinson's disease; we quote [46] for the study of a model for collective cell movement due to chemical sensing; we quote [47] for applications to biological populations in competition; we quote [48] for the investigation of a fractional order tumor model; we quote [49] for the optimal control for a nonlinear mathematical model of tumor under immune suppression; we quote [50] for the study of a fractional reaction-diffusion equation with a nonlocal boundary condition modeling the invasion of a tumor and its growth. Now, coming back to our system (1)-(3) and referring to the paper [1], it turns out that, under rather general assumptions on the operators and the potentials, well-posedness and regularity results for the initial value problem for (1)- (3) were established in the case u = 0, under the assumption that α > 0 and β > 0. However, some remarks on more general cases including u ∈ L 2 (Q), where Q := Ω × (0, T), have been given in [1]. In particular, under suitable assumptions on the initial data, for every u ∈ L 2 (Q), there exists at least a solution (µ, ϕ, S) in a proper functional space to a weak version of (1)-(3) (namely, Equation (2) is replaced by a variational inequality involving the convex part F 1 of F rather than f , since F 1 is not supposed to be differentiable). Moreover, the solution is unique if the domains of the fractional operators A ρ and C τ satisfy suitable embeddings of Sobolev type. Finally, if B 2σ behaves like the Laplace operator with either Dirichlet or Neumann zero boundary conditions and f is single valued (like (4) and (5)), then the solution solves Equation (2) in a stronger sense, and it is even smoother under more restrictive assumptions on the initial data.
In this paper, we first establish similar results for system (1)-(3) by assuming α = 0 and β > 0 (in fact, we take β = 1 without loss of generality). In particular, we extend some results shown for this case in the recent paper [38]. Then, we discuss a distributed control problem for the modified system. Namely, given nonnegative constants κ i , i = 1, . . . , 5, and functions ϕ Q , S Q ∈ L 2 (Q) and ϕ Ω , S Ω ∈ L 2 (Ω), we consider the problem of minimizing the cost functional where ϕ and S are the components of the solution (µ, ϕ, S) corresponding to the control u, which is supposed to vary under restrictions of the type u min ≤ u ≤ u max . The choice of this tracking-type cost functional reflects the plan of a medical treatment via the application of drugs over some finite time interval (0, T) with the aim of monitoring the evolution of the tumor fraction ϕ under the restriction that no harm be inflicted on the patient. We remark at this place that it would be desirable to minimize the duration, i.e., the time T > 0, of the medical treatment as well, in order to prevent that the tumor cells develop a resistance against the drug. However, such an approach, which was possible (see, e.g., [21]) in the special case when A 2ρ = B 2σ = C 2τ = −∆, becomes very complicated in the situation considered here and was therefore not included.
The literature on optimal control problems for Cahn-Hilliard systems is still scarce. In this connection, we refer the reader to [31,51], where a number of references is given. Even less investigations have been made on optimal control problems for tumor growth models. About that, let us refer to the works [18,21,24,[52][53][54][55][56][57][58][59], for various models involving the Laplacian. Concerning the optimal control of Cahn-Hilliard systems with fractional operators, we just can cite [31,32], and, to the authors' best knowledge, the present paper is the first contribution on the optimal control of the tumor growth model with fractional operators.
The remainder of the paper is organized as follows. In the next section, we list our assumptions and notations and present our results on the state system. Section 3 is devoted to the study of the control-to-state mapping and of its Fréchet differentiability. In the last section, we deal with the control problem. Namely, the existence of an optimal control is proved and the first order necessary conditions involving a proper adjoint system are derived.

The State System
In this section, we first introduce the notations and the assumptions needed for the analysis of the state system. Then, we present our results. We closely follow [1]. First, of all, the set Ω ⊂ R 3 is assumed to be bounded, connected and smooth, with volume |Ω| and outward unit normal vector field ν on Γ := ∂Ω. Moreover, ∂ ν stands for the corresponding normal derivative. We set H := L 2 (Ω) (6) and denote by · and ( · , · ) the standard norm and inner product of H. As for the operators, we first postulate that unbounded monotone self-adjoint linear operators with compact resolvents.
Therefore, there are sequences {λ j }, {λ j }, {λ j } and {e j }, {e j }, {e j } of eigenvalues and of corresponding eigenvectors satisfying Ae j = λ j e j , Be j = λ j e j , and Ce j = λ j e j , with (e i , e j ) = (e i , e j ) = (e i , e j ) = δ ij for i, j = 1, 2, . . . , {e j }, {e j } and {e j } are complete systems in H.
As a consequence, we can define the powers of the above operators with arbitrary positive real exponents. As far as the first one is concerned, we have, for ρ > 0, and we endow V Similarly, we set with the graph norms From now on, we assume: ρ, σ and τ are fixed positive real numbers.
However, we need the further assumptions we list at once. It is understood that all of the embeddings below are assumed to be continuous.
The first eigenvalue λ 1 of A is strictly positive.
B and every monotone and Lipschitz continuous function ψ : R → R vanishing at the origin. (19) Due to the continuus embeddings (18), there exists a constant C * > 0 such that where, for p ∈ [1, +∞], the symbol · p denotes the norm in L p (Ω). The same symbol will also be used for the norm in L p (Q) provided that no confusion can arise.

Remark 1.
We have to make some comments on (17)- (19). The first of these assumptions is satisfied if A is, e.g., the Laplace operator −∆ with zero Dirichlet (or Robin) boundary conditions, while the case of zero Neumann boundary conditions is excluded unless one adds to the Laplace operator, e.g., some zero-order term ensuring coerciveness. However, it is clear that A could be a much more general operator. By still considering the Laplace operator with (zero) Dirichlet boundary conditions as A, we can also discuss the first two embeddings in (18). By noting that D(A) = H 2 (Ω) ∩ H 1 0 (Ω) and Ω is smooth, it results that D(A 2ρ ) ⊂ H 4ρ (Ω) and D(A ρ ) ⊂ H 2ρ (Ω). Hence, both embeddings hold true if ρ ≥ 3/8, since Ω is three-dimensional. Finally, we make a comment on (19). Assume, for instance, that B 2σ = −∆ with zero Neumann boundary conditions. Then, V 2σ B = {v ∈ H 2 (Ω) : ∂ ν v = 0 on Γ} and, for every v ∈ V 2σ B and ψ as in (19), we have that ψ(v) ∈ H 1 (Ω) (since v ∈ H 1 (Ω)) and The same argument works if we take the Dirichlet boundary conditions instead of the Neumann ones, since the functions ψ considered in (19) vanish at the origin. More generally, B 2σ can be the principal part of an elliptic operator in divergence form with smooth coefficients. In particular, even though some restrictions on A, B, and C have to be imposed in order to fulfill the properties (17)- (19), no relationship between them is needed, and the three operators can be completely independent from each other.

Remark 2. Assumption (17) allows us to consider an equivalent norm in
so that the function v → A ρ v defines a norm in V ρ A that is equivalent to the graph norm (13).
For the nonlinear functions entering our system, we postulate the following properties: for some constants C i > 0 and every s ∈ D(F).

Remark 3.
The hypotheses (23)-(26) on F ensure that the conditions required in [1], i.e., F = F 1 +F 2 | (a,b) , wherẽ F λ 1 (s) +F 2 (s) ≥ −C 0 for some constant C 0 and every s ∈ R, (30) are satisfied, as we show at once. We first split F by defining, for s ∈ (a, b), Notice that F 1 is nonnegative and convex and that F 1 (0) = 0. If (a, b) = R, we properly extend these functions F i to functionsF i defined in the whole of R. One can preserve the mentioned properties of F 1 , including its lower semicontinuity, by setting Moreover, one can ensure that the derivative of the extensionF 2 is Lipschitz continuous, by noting that F 2 = f 2 already is Lipschitz continuous in (a, b) since its derivative f 2 = −( f ) − is bounded by the assumption (25) on F . The last condition that we have to check is (30), whereF λ 1 is the Moreau-Yosida approximation ofF 1 at the level λ. We notice that this condition is not equivalent to an inequality of type F(s) ≥ −C 2 , which follows from (25) and looks rather natural in performing formal a priori estimates. On the other hand, one can prove that the inequality we need is implied by the full quadratic growth condition given in (25) (see (formula (3.1) [33] for some explanation). For this reason, we have postulated the latter.
Although some of the results to be presented will not require the whole set of hypotheses made so far, the statements will be greatly simplified if we do not each time recall the properties of the involved operators, spaces, and nonlinearities; we therefore make the following general assumption: All of the assumptions made above on the structure are in force from now on.
As mentioned in the Introduction, we only deal with the case α = 0 and β > 0 of system (1)-(3). Clearly, we can take β = 1 without loss of generality. Hence, the Cauchy problem forming the state system under investigation reads as follows: where ϕ 0 and S 0 are prescribed initial data that are supposed to satisfy In fact, we could solve a weak form of the above problem under milder assumption on the initial data; however, in order to guarantee a sufficient regularity level of the solution, we need the whole of (36). Given a final time T ∈ (0, +∞), the regularity we can ensure (besides some boundedness to be discussed later on) is the following: so that Equations (32)-(34) are satisfied a.e. in Q, where we recall that We notice at once that the first embedding in (18) yields that At this point, we are ready to present our results. To this end, it is convenient to introduce the following variational formulation of (32)- (34): This is based on obvious properties of the powers of the operators A, B and C, like the Green type Before stating our well-posedness theorem, we prove some auxiliary results. The first one is a separation property enjoyed by any solution under our assumptions on the data. Theorem 1. Assume (36) and u ∈ L 2 (0, T; H), and let (µ, ϕ, S) be a solution to problem (32)- (35) satisfying (37)- (40). Then, it holds for every This interval depends only on f , the initial datum ϕ 0 , and M.
Proof. Notice that µ is bounded thanks to (42). Thus, we fix a constant M and assume that µ ∞ < M. By the assumptions (26) on f and (36) on ϕ 0 , we can choose a M ∈ (a, Now, we notice that, for a.a. s ∈ (0, T), the value ϕ(s) belongs to V 2σ B by (38). Moreover, the function z → ψ(z) := (z − b M ) + is monotone and Lipschitz continuous on R and vanishes at the origin. Hence, we have that ψ(ϕ(s)) ∈ H by (19). Thus, we can multiply (33), written at the time s, by ψ(ϕ(s)) and integrate over (0, t) with respect to s. By noting that ψ(ϕ 0 ) = 0 (since in Ω), we obtain that Thanks to the inequality (19), the second term on the left-hand side is nonnegative. Moreover, the right-hand side is nonpositive since
Proof. We show the uniqueness at the end and first prove the estimate (47), noting that the assumption µ i ∞ < M is meaningful since the functions µ i are bounded due to (42). We apply Theorem 1 and Let L be the corresponding Lipschitz constant. After this preparation, we can start the proof. We set, for convenience, u := u 1 − u 2 , µ := µ 1 − µ 2 , ϕ := ϕ 1 − ϕ 2 , and S := S 1 − S 2 , writing (32)-(34) for both solutions and multiplying the differences by µ, ∂ t ϕ, and S, respectively, in the inner product of H. Then, we sum up and integrate with respect to time over (0, t). By noting that the terms involving (µ, ∂ t ϕ) cancel each other, and adding the same contributions (1/2) ϕ(t) 2 = t 0 (ϕ(s), ∂ t ϕ(s)) ds and t 0 S(s) 2 ds to both sides, we obtain the equation If we term I the first integral on the right-hand side, apply the Young inequality to the next three terms, use the Lipschitz continuity of f , and rearrange, we deduce that Now, we rewrite I as the sum of two terms. The first of these is nonpositive, since P is nonnegative, and we estimate the other one recalling that P is bounded because P is Lipschitz continuous. By applying the Hölder inequality, and recalling (20), we have for every δ > 0 that We notice that the function s → S 2 (s) C, τ + µ 2 (s) A, ρ 2 belongs to L ∞ (0, T), thanks to the regularity (37) and (40) of µ 2 and S 2 , respectively. Therefore, by choosing δ > 0 small enough, and applying the Gronwall lemma, we obtain the estimate (47) with a constant K M whose dependence on the data agrees with that specified in the statement. We now come back to uniqueness. As both u i ∈ L 2 (0, T; H) and the corresponding solutions (µ i , ϕ i , S i ) are arbitrary in the above argument (since no restriction on M is made), we conclude that (µ 1 , ϕ 1 , S 1 ) = (µ 2 , ϕ 2 , S 2 ) if u 1 = u 2 , which shows the uniqueness for the solution to problem (32)- (35). With this, the proof is complete.
Finally, we can state our well-posedness and stability result. Here, and later on in this paper, we use the notation where R is a positive real parameter.
Moreover, for every R > 0, there exist a constant K 1 (R) and a compact interval [a R , b R ] ⊂ (a, b), which depend only on the structure of the system, the initial data, T and R, such that both the estimate (49) and the separation property hold true for every u ∈ B R and the corresponding solution (µ, ϕ, S). Finally, if R > 0, u i ∈ B R , i = 1, 2, and (µ i , ϕ i , S i ) are the corresponding solutions, then the estimate holds true with a constant K 2 (R) that depends only on the structure of the system, the initial data, T, and R.
Proof. Uniqueness follows from Theorem 2. Let us come to the existence of a solution and to the estimates (49) and (51). However, we do not give a complete proof. Indeed, one can adapt the arguments of [1] on account of Remark 3, and we briefly explain the reason for this. The procedure used there is based on the Yosida regularization of the nonlinearity f , a time discretization of the regularized system, and the derivation of suitable a priori estimates, and the same line of argumentation can be followed in our situation. Here, is the main remark: in [1], some estimates for µ have been derived from estimates of ∂ t µ, and this term is missing in (32), in contrast to (1). In the present case, an estimate of the norm of µ, e.g., in L 2 (0, T; H), can be deduced from an estimate of A ρ µ in the same space as shown in Remark 2, by using the assumption (17) on the first eigenvalue λ 1 of A. Hence, we do not repeat the arguments of [1] with the corresponding modifications. However, for the reader's convenience, we sketch the formal proofs of the estimates that would be obtained step by step in the rigorous procedure in order to prove the existence of a solution. We assume u ∈ B R from the very beginning, so that these estimates eventually lead to (49) as well. In order to simplify notation, we use the same symbol c without any subscript for possibly different constants that depend only on the structure of our system, the initial data and T, but neither on u nor on R. Moreover, the symbol c R stands for (possibly different) constants that depend on the constant R, in addition, but still not on u. Thus, it is understood that the actual values of such constants may vary from line to line and even in the same chain of inequalities. Notice that the notations used for constants we want to refer to (like, e.g., those used in (25)) are different.
First a priori estimate. We test (43), (44) and (45), written at the time s, by µ(s), ∂ t ϕ(s), and S(s), respectively, in the scalar product of H. Then, we sum up and integrate over (0, t), where t ∈ (0, T) is arbitrary, noting that the terms involving the product µ ∂ t ϕ cancel each other. By also adding |Ω|C 2 to both sides (see (25)), we obtain the identity Recalling the consequence (22) of (17), taking (25) into account and applying the Gronwall lemma, we conclude that Consequence. We can also infer that Indeed, the estimate is right for the first term, since P is bounded by (27); for the bound of the second term, we can argue as in (Section 4.4 [1]).
Second a priori estimate. We would like to test (43) by ∂ t µ even though ∂ t µ does not appear in the equation (in contrast to (1)). In fact, the estimate we derive here by a formal procedure should be performed rigorously at the level of the discete scheme, which can contain that time derivative multiplied by a viscosity coefficient that tends to zero at some point of the procedure. Thus, we test (32) and (34) formally by ∂ t µ and ∂ t S, respectively. At the same time, we formally differentiate (33) with respect to time and test the resulting equality by ∂ t ϕ. Then, we sum up and integrate with respect to time, as usual. Since the terms involving the product ∂ t µ ∂ t ϕ cancel each other, we obtain the identity The integral containing u can be handled using Young's inequality, and the one involving f is easily treated using the second inequality in (25) and (52). We postpone the estimate of the last integral and first deal with the initial values appearing on the right-hand side of (53). By using the initial conditions for ϕ and S, we write (32) and (33) at the time t = 0 in the following way: Since P is nonnegative, by multiplying the first identity in (54) by µ(0), we have that (see (22)) and, on account of (36), we also deduce from the second identity in (54) that Finally, we deal with the last integral on the right-hand side of (53), which we term I for brevity. We perform an integration by parts in time, recall that P is bounded by (27), and invoke (55). By recalling that · p denotes the norm in L p (Ω), we then have 2 4 ds .
Finally, we notice that (52) and some of the embeddings in (18). This allows us to apply Gronwall's lemma, whence we conclude that Third a priori estimate. By taking v = µ(t) in (43) and recalling that P is nonnegative and bounded, we obtain the following inequality for a.a. t ∈ (0, T) By accounting for Remark 2, we deduce that In particular, the norm of µ in L ∞ (0, T; H) is bounded by some constant c R . Since the same holds for S due to (52), we infer that Consequence. By comparison in (32), we deduce that Combining this with (58) and (56), we conclude that Then, the first embedding in (18) yields that µ is bounded (as claimed in the statement) and that By comparison in (34), we also deduce that No better estimate for S is available, since u ∈ L 2 (0, T; H) only.
Fourth a priori estimate. We recall that Remark 3 provides a splitting of f as f 1 + f 2 , with f 1 monotone and vanishing at the origin and f 2 Lipschitz continuous. Thus, we can write (33) for a.a. t ∈ (0, T) in the form and test this identity by f 1 (ϕ(t)). More precisely, in the correct argument, f 1 is replaced by its Yosida regularization, which is monotone and Lipschitz continuous and vanishes at the origin, and the equation itself is replaced by a scheme, which is obtained by discretizing time differentiation and for which the analogue of ϕ(t) belongs to V 2σ B . Hence, assumption (19) can actually be applied. Here, we formally apply it to the above identity with v = ϕ(t) and ψ = f 1 . We obtain that On account of the previous estimates, and by a comparison in (33), we conclude that Conclusion. This concludes the formal proof of the existence part of Theorem 3 and of estimate (49). As already said, in the rigorous argument, the above bounds are established for the solution to an approximating problem, and one has to perform some limiting procedure. The estimates provide convergence of weak and weak-star type. However, even strong convergence in L 2 (0, T; H) for the approximations of ϕ and S is obtained. Indeed, the embeddings V σ B ⊂ H and V τ C ⊂ H are compact due to (7), so that one can apply the Aubin-Lions lemma (see, e.g., (Thm. 5.1, p. 58 [60])). Therefore, the nonlinear terms can be correctly managed.
Separation. Let us come to estimate (50). This is a trivial consequence of the above estimate and Theorem 1. Indeed, this theorem can be applied with M := c R + 1, where c R is the constant that appears in (60). The corresponding compact interval [a R , b R ] of the statement is nothing but the interval [a M , b M ] considered in (46), which depends only on the structure of the system, the initial data, T, and R.
Continuous dependence. In addition, (51) is a trivial consequence of a fact already proved, namely, of Theorem 2. Indeed, if u i ∈ B R , i = 1, 2, then the L ∞ bound for the corresponding µ i is ensured by (49), and Theorem 2 can be applied with M = K 1 (R) + 1. Hence, we can take as K 2 (R) the constant K M that appears in (47). In addition, this constant depends only on the structure of the system, the initial data, T and R.

Remark 4.
The existence part of Theorem 2.6 is closely connected to the existence result proved in (Theorem 3.4 [38]), where, however, no statement concerning separation or uniqueness was proved. For purposes of control theory, however, it is indispensable to have uniqueness, since otherwise no control-to-state operator can be defined, and this seems to be available only under the assumptions made here.

The Control-to-State Mapping
The results of the previous section ensure that we can correctly define a control-to-state mapping to be used in the control problem under investigation. Taking into account that the cost functional to be minimized depends only on the components ϕ and S of the solution corresponding to a given u, we set More precisely, we need to consider the restriction of these maps to B R for any given radius R > 0. The choice of the space Y mainly is due to the following fact: the inequality (51) implies that A very important consequence of the separation property (50) and of the regularity of f ensured by (24) is the following global boundedness condition: where K 3 (R) depends only on the structure of the system, the initial data, T, and R. The Fréchet differentiability of the maps S is strictly related to the properties of the linearized problem we introduce now. To this end, we fix u ∈ L 2 (0, T; H). The linearized system associated with u and the variation h ∈ L 2 (0, T; H) is the following: where µ := S 1 (u), ϕ := S 2 (u), and S := S 3 (u). We also write the weak formulation of (67)-(69): the identities have to hold true for every v ∈ V ρ A , v ∈ V σ B , and v ∈ V τ C , respectively, and for a.a. t ∈ (0, T). We have the following results.
Theorem 4. Suppose that the assumptions (36) on the initial data of problem (32)- (35) are fulfilled, and let u ∈ L 2 (0, T; H) and h ∈ L 2 (0, T; H). Then, the linearized problem (67)-(70) has a unique solution (η, ξ, ζ) satisfying the regularity requirements Moreover, if R > 0 and u ∈ B R , then this solution satisfies the estimate where the constant K 4 (R) depends only on the structure of the system (32)-(34), the initial data ϕ 0 and S 0 , T, and R.
Proof. We notice that the coefficients P(ϕ), P (ϕ), f (ϕ), as well as µ, are bounded functions. Moreover, if u belongs to some B R , then the L ∞ bounds are uniform, i.e., they just depend on R and not on u.
On the contrary, S might be unbounded. However, as stated in Theorem 3, it is smooth. Thus, the linear system is not worse than the nonlinear one and can be solved by the same argument (which we do not repeat here) based on time discretization that has been used in [1] (see also [31] for the linearized system associated with the Cahn-Hilliard equations). This first leads to a solution to the variational problem (71)-(73) and then to the strong formulation (67)-(69). However, we perform at least some formal estimates that can justify both the regularity asserted in the statement and the validity of estimate (77). To this end, we fix R > 0 and assume that u ∈ B R at once. In addition, in this section, i.e., in this proof and later on, we adopt a convention on the constants similar to the one used in the previous section: c stands for possibly different constants depending only on the structure, the data, and T, while the notation c R indicates an additional dependence on R.
First a priori estimate. We formally test (71)-(73) by η, ∂ t ξ, and ζ, respectively, sum up and integrate over (0, t). Moreover, we add the same quantities (1/2) ξ(t) 2 = t 0 (ξ(s), ∂ t ξ(s)) ds and t 0 ζ(s) 2 to both sides of the resulting identity, in order to recover the full norms in V σ B and V τ C on the left-hand side. In addition, in this case, a cancellation occurs, and we have that The first integral on the right-hand side is nonpositive, while the next one, which we term I, needs some treatment. By using the Hölder inequality, two of the inequalities (20), Remark 2, and the Young inequality, we obtain that where we notice that the function s → S(s) 2 C, τ + µ(s) 2 A, ρ belongs to L ∞ (0, T) and that its norm is bounded by a constant like in (49), due to Theorem 3 applied to u. By treating the last terms of (78) using the Schwarz and Young inequalities, and applying Gronwall's lemma, we conclude that Second a priori estimate. We estimate the right-hand side of (67). On account of the embeddings (18), we have a.e. in (0, T) that On account of (79), we conclude that Since (79) also yields an estimate for ∂ t ξ, a (formal) comparison in (67) allows us to conclude that Third a priori estimate. We test (69) by ∂ t ζ and integrate in time, as usual. On account of (80) and the Young inequality, we obtain that We thus deduce that Now that ∂ t ζ is estimated, a comparison in (69) provides a bound for C 2τ ζ. Hence, we conclude that This ends the list of the formal estimates and formally leads to a strong solution satisfying (77). Even though uniqueness formally follows by taking h = 0, we remark that it can be proved rigorously. Indeed, by assuming the regularity (74)-(76), the procedure used to obtain the above estimates is justified. Proof. Fix any u ∈ L 2 (0, T; H), and let (µ, ϕ, S) be the corresponding state. For every h ∈ L 2 (0, T; H), let (µ h , ϕ h , S h ) be the state corresponding to u + h. Finally, let (η, ξ, ζ) be the solution to the linearized problem (67)-(70) associated with u and h. We set, for convenience, According to the definitions of differentiability and derivative in the sense of Fréchet, we have to prove that the (linear) map h → (ξ, ζ) is continuous from L 2 (0, T; H) into Y and that there exist a real number h > 0 and a function Λ : (0, h) → R satisfying The first fact is ensured by (77) once R is chosen larger than u L ∞ (0,T;H) . Hence, we fix R > u L 2 (0,T;H) once and for all. As for the construction of Λ, we set h := R − u L 2 (0,T;H) , and we assume that h L 2 (0,T;H) < h. This implies that u and u + h belong to B R , so that Theorem 3 can be applied to both of them. We thus derive uniform estimates for the corresponding states, hence for the coefficients of the corresponding linearized systems. This entails uniform estimates for the corresponding solutions. In order to establish (85), we observe that (η h , ξ h , ζ h ) satisfies the regularity properties and solves the problem (in a strong form, i.e., the equations are satisfied a.e. in Q, since all the contributions are L 2 functions) where Q h 1 and Q h 2 are defined by It is convenient to rewrite the functions Q h i by accounting for the Taylor expansions of P and f . Using the formula with integral remainder, it is immediately checked that while it is more complicated to find a convenient representation of Q h 1 . However, simple algebraic manipulations show that where Notice that, by (66), both R h 1 and R h 2 are bounded uniformly with respect to h: After this preparation, we start estimating. To this end, we test (86), (87), and (88), by η h , ∂ t ξ h , and ζ h , respectively. Then, we sum up and integrate in time. There is a usual cancellation. By adding the same contributions to both sides similarly as in the previous proof, we obtain We have to estimate only the integrals involving Q h 1 and Q h 2 . The first term produces four integrals, termed I j for j = 1, . . . , 4 for brevity, which correspond to the four summands, in that order, of (91). Clearly, I 1 is nonpositive. As for I 2 , we use the Hölder inequality and the embeddings (18) as well as the estimate (51) applied with u 1 = u + h and u 2 = u. Hence, by omitting the integration variable s to shorten the lines, we have, for every δ > 0, Next, by also accounting for (49) applied to µ and S, we similarly have that and for the fourth contribution, thanks to (92), we obtain that Finally, we estimate the term involving Q h 2 by accounting for (90) and (92) in this way: By treating the last two terms of (93) in a trivial way, recalling all the inequalities derived above, choosing δ > 0 small enough, and applying the Gronwall lemma, we conclude that If we term C R the value of the constant c R of the last inequality, then we obtain (85) with Λ defined on (0, h) by Λ(s) := C R s 2 . This completes the proof.

The Control Problem
As announced in the Introduction, the main aim of this paper is the discussion of a control problem for the state system studied in the previous sections. For this problem, we assume that u min , u max ∈ L ∞ (Q), and u min ≤ u max a.e. in Q.
Then, the cost functional J and the set U ad of the admissible controls are defined by For the above problem, we prove the existence of an optimal control, and we derive the first-order necessary conditions for optimality. This involves an adjoint problem for which we prove a well-posedness result. We recall that the control-to-state mapping S is defined in (64) and state our first result. Theorem 6. Under the assumptions (94)-(95), the control problem has at least one solution, that is, there is some u ∈ U ad satisfying the following condition: for every v ∈ U ad we have that J (u, ϕ, S) ≤ J (v, ϕ, S), where (ϕ, S) = S(u) and (ϕ, S) = S(v).
Proof. Since U ad is nonempty, the infimum of J under the constraints given in (98) is a well-defined real number d ≥ 0, and we can pick a minimizing sequence {u n } ⊂ U ad . Hence, denoting by (µ n , ϕ n , S n ) the state corresponding to u n for n ∈ N, we have that J (u n , ϕ n , S n ) → d as n → ∞.
Since U ad is bounded and closed in L 2 (0, T; H) (in fact, it is even bounded and closed in L ∞ (Q)), we can assume that u n → u weakly in L 2 (0, T; H) for some u ∈ U ad . Moreover, we can choose some R > 0 such that U ad ⊂ B R . Therefore, we can apply Theorem 3 to u n and deduce that (µ n , ϕ n , S n ) satisfies the estimate (49), as well as the separation and global boundedness properties (50) and (66), for all n ∈ N. Hence, for a subsequence indexed again by n, we have that It follows that the initial conditions (35) are satisfied by the limiting pair (ϕ, S). Moreover, thanks to the compact embedding V σ B ⊂ H ensured by (7), and consequently of H 1 (0, T; V σ B ) into L 2 (0, T; H), we deduce that ϕ n → ϕ strongly in L 2 (0, T; H).
Since f and P are Lipschitz continuous in [a R , b R ], we also infer that f (ϕ n ) → f (ϕ) and P(ϕ n ) → P(ϕ), strongly in L 2 (0, T; H).
It follows that (µ, ϕ) solves (33). From the above strong convergence and the weak convergence of {µ n } and {S n } at least in L 2 (0, T; H), we deduce that {P(ϕ n )(S n − µ n )} converges to P(ϕ)(S − µ) weakly in L 1 (Q). Hence, the limiting triplet (µ, ϕ, S) satisfies Equations (32) and (34) as well, i.e., (µ, ϕ, S) is the state corresponding to the control u. On the other hand, we have that J (u, ϕ, S) ≤ lim inf n ∞ J (u n , ϕ n , S n ) = d by semicontinuity. We conclude that u is an optimal control. The rest of the section is devoted to the derivation of the first-order necessary conditions for optimality. Hence, we fix an optimal control u ∈ U ad and the corresponding (µ, ϕ, S) once and for all. If we introduce the so-called reduced cost functionalJ by setting J (u) := J (u, S 2 (u), S 3 (u)) for u ∈ L 2 (0, T; H), we immediately find from the convexity of U ad that the Fréchet derivative (DJ )(u) ∈ L(L 2 (0, T; H); R) must satisfy provided that it exists. However, this is the case due to the obvious differentiability of the quadratic functional J and the differentiability of the operator S, which takes its values in Y ⊂ (C 0 ([0, T]; H)) 2 . Hence, by accounting for the full statement of Theorem 5, we can even apply the chain rule and rewrite the above inequality as where ξ and ζ are the components of the solution (η, ξ, ζ) to the linearized system (67)-(70) associated with u and h = v − u.
As usual in control problems, a condition of this sort is not satisfactory, since it requires solving the linearized problem for infinitely many choices of h ∈ L 2 (0, T; H) because v is arbitrary in U ad . Therefore, we have to eliminate ξ and ζ from (103), which can be done by introducing and solving a proper adjoint problem. This is a backward-in-time problem for the adjoint state variables (q, p, r) that formally reads as follows: However, in order to give this system a proper meaning according to the regularity that we will prove, we need some preliminaries. First, due to the density of V σ B in H, we can identify H with a subspace of the dual space in such a way that v, w = (v, w) for every v ∈ H and w ∈ V σ B , where · , · denotes the duality pairing between V −σ B and V σ B . Now, thanks to the obvious formula (B 2σ v, w) = (B σ v, B σ w), which holds for every v ∈ V 2σ B and w ∈ V σ B and, owing to the above identification, can also be read in the form B 2σ v, w = (B σ v, B σ w), one can extend the operator B 2σ : V 2σ B → H to a continuous linear operator, still termed B 2σ , from V σ B to V −σ B by means of the above formula, namely, At this point, it is meaningful to postulate the following regularity for the adjoint variables: Indeed, then all of the equations, as well as the final conditions, have a precise meaning, by also accounting for the properties of the other ingredients which we recall for the reader's convenience: P(ϕ), P (ϕ), f (ϕ), and µ are bounded, and S ∈ H 1 (0, T; H) ∩ L ∞ (0, T; V σ B ), whence, in particular, S ∈ L ∞ (0, T; L 4 (Ω)). However, we also consider a variational formulation of the adjoint system, which makes sense in a much weaker regularity setting for (q, p, r), namely, We require that where we have introduced the abbreviating notation In addition, for brevity, and in order to shorten the exposition, we have omitted the integration time variable termed s. We will do the same in the following.
Then, the function t → (y(t), z(t)) H is absolutely continuous on [0, T], and its derivative is given by where ( · · ) H and V * · · V denote the inner product in H and the dual pairing between V * and V, respectively. Proof. We first notice that (113) implies the pointwise variational inequality A and a.e. in (0, T).
On the other hand, the conditions w ∈ V ρ A , g ∈ H, and (A ρ w, A and A 2ρ w = g, as one immediately sees by using the spectral representation. Hence, we obtain (109) and (104). The same argument can be used to deduce that r belongs to L 2 (0, T; V 2τ C ) and solves (106), since even ∂ t r belongs to L 2 (0, T; H) by assumption. The last condition r ∈ L ∞ (0, T; V τ C ) in (111) then follows from interpolation. Much more work has to be done for the second equations. First, for the same test functions v as in (114), we deduce that where, for brevity, we have set We immediately infer that In particular, we have for some constant c > 0 that This exactly means that ∂ t (q + p) ∈ (L 2 (0, T; V σ B )) * = L 2 (0, T; V −σ B ). Thus, we can replace the . The variational equation we obtain is just (105) understood in the sense of V −σ B . It remains to derive the first of the final conditions (107). To this end, we also assume that v(0) = 0 and exploit (120) once more. Moreover, we can apply Lemma 1 with . Finally, we account for the already proved Equation (105). We then obtain that Therefore, we have that (q + p)(T), v(T) = g 2 , v(T) for every v with the required properties, and the desired final condition obviously follows.
Thus, we can choose between the strong form (104)-(107) and the weak formulation (113)-(116), according to our convenience, in proving a well-posedness result, which is our next goal. We prepare the existence part by introducing a Faedo-Galerkin scheme with viscosity that looks like an approximation of (104)-(107). We recall (8)  Then, we look for a triplet (q n , p n , r n ) satisfying and solving the system − 1 n ∂ t q n + A 2ρ q n − p n + P(ϕ)(q n − r n ), v = 0 for every v ∈ V (n) A and a.e. in (0, T), −∂ t (q n + p n ) + B 2σ p n + f (ϕ)p n − P (ϕ)(S − µ)(q n − r n ), v = (g 1 , v) for every v ∈ V (n) B and a.e. in (0, T), −∂ t r n + C 2τ r n − P(ϕ)(q n − r n ), v = (g 3 , v) for every v ∈ V (n) C and a.e. in (0, T), as well as the final conditions The following result holds true.
Proof. The requirements (121) mean that q n (t) = n ∑ j=1 q n j (t) e j p n (t) = n ∑ j=1 p n j (t) e j and r n (t) = n ∑ j=1 r n j (t) e j for a.a. t ∈ (0, T) and some functions q n j , p n j , r n j ∈ H 1 (0, T). Moreover, an equivalent system is obtained by taking for i = 1, . . . , n just v = e i , v = e i , and v = e i , in the three variational equations, respectively. Hence, (122)-(124) become an ODE system having the column vectors q n := (q n j ), p n := (p n j ), and r n := (r n j ), as unknowns. This system reads as follows: −(e j , e i ) d dt r n j + (λ j ) 2τ (e j , e i )r n j − P(ϕ)e j , e i q n j + P(ϕ)e j , e i r n j = (g 3 , e i ), where the index i runs over {1, . . . , n} in all of the equations, which are understood to hold a.e. in (0, T). Thus, thanks to the orthogonality conditions in (8), it takes the form − 1 n q n + M 1 q n + M 2 p n + M 3 r n = 0 M 4 q n − p n + M 5 q n + M 6 p n + M 7 r n = b n − r n + M 8 q n + M 9 r n = b n for some (possibly time dependent, but bounded) (n × n) matrices M k , k = 1, . . . , 9, and column vectors b n b n ∈ L 2 (0, T; R n ). Therefore, one can solve the first equation for q n and replace q n in the second one by the resulting expression. At the same time, one multiplies the first equation by n and keeps the third one as it is. This procedure leads to an equivalent system of the form −y + My = b for some matrix M ∈ L ∞ (0, T; R 3n×3n ) and some vector b ∈ L 2 (0, T; R 3n ) in the unknown y ∈ H 1 (0, T; R 3n ) obtained by rearranging the triplet (q n p n r n ) as a 3n-column vector.
On the other hand, the final conditions (125) provide a final condition for y. Hence, standard results for ODEs show the unique solvability.
At this point, we are ready to solve the adjoint problem. We need, however, the following additional compatibility condition: Remark 5. The compatibility condition (126) is satisfied if either κ 4 = 0 or S Ω ∈ V τ C . Obviously, κ 4 = 0 means that we do not have a tracking of the solution variable S at the final time T; while this is not desirable, it is not too much of a restriction, since one is rather interested in monitoring the final tumor fraction ϕ(T) than S(T). On the other hand, the assumption S Ω ∈ V τ C is not overly restrictive in view of the fact that S ∈ H 1 (0, T; H) ∩ L 2 (0, T; V 2τ C ), whence it follows that S ∈ C 0 ([0, T]; V τ C ) by continuous embedding, and thus S(T) ∈ V τ C . To assume the same regularity for S Ω is certainly not unreasonable.
We have the following result. Proof. In order to prove the existence of a solution, we start from the finite-dimensional problem (122)-(125), perform an a priori estimate, and let n tend to infinity. In addition, in this section, we simplify the notation as far as constants are concerned and use the same symbol c for different constants that can depend only on the structure, the data, T, the optimal control u, and the corresponding state (µ, ϕ, S). A priori estimate. We write the Equations (122)-(124) at the time s and test them by −∂ t q n (s), p n (s), and −∂ t r n (s), respectively. Then, we sum up, integrate over (t, T) with respect to s, and notice that the terms involving the product p n ∂ t q n cancel each other. Moreover, we add the same quantities whence we infer that C τ r n (T) ≤ C τ g 4 . In conclusion, we have shown the estimate r n (T) C,τ ≤ g 4 C,τ for all n ∈ N.
Next, we consider the first two terms on the right-hand side, which we denote by Y 1 and Y 2 . We only need to estimate these terms, since the remaining other ones can easily be handled using Young's inequality and, eventually, Gronwall's lemma. As for Y 1 , we first integrate by parts, and one of the important terms we obtain is nonpositive. Then, we account for the Hölder and Young inequalities, the equivalence of norms in V ρ A related to (22), and the embeddings (18) as follows: Concerning Y 2 , we have that By treating the remaining terms on the right-hand side of (127) as announced before, and applying the Gronwall lemma, we conclude that Existence. The above estimate ensures that, for a subsequence again indexed by n, 1 n ∂ t q n → 0 strongly in L 2 (0, T; H), q n → q weakly star in L ∞ (0, T; V ρ A ), p n → p weakly star in L ∞ (0, T; H) ∩ L 2 (0, T; V σ B ), r n → r weakly star in H 1 (0, T; H) ∩ L ∞ (0, T; V τ C ).
We aim at proving that (q, p, r) is the desired solution to the weak form (113)-(116) of the adjoint problem. Clearly, (116) is satisfied, and we have to prove that the variational equations are satisfied as well. We confine ourselves to the second equation, which is the most complicated one. To this end, we write an integrated version of (123). We fix any integer m > 1, take any v ∈ H 1 (0, T; H) ∩ L 2 (0, T; V A , so that v(s) is admissible in (123) written at the time s, and we can test the equation in the inner product of H. Moreover, we replace (B 2σ p n (s), v(s)) by (B σ p n (s), B σ v(s)). Then, we integrate over (0, T) with respect to s. Now, we observe that v(T) is admissible in the second identity of (125). Thus, by an integration by parts, we obtain that T 0 (q n + p n , ∂ t v) + (B σ p n , B σ v) + f (ϕ) p n − P (ϕ)(S − µ)(q n − r n ), v ds = T 0 (g 1 , v) ds + g 2 , v(T) .
As the other terms are easier, we conclude that (114)  ) and vanishes at t = 0. Hence, we can use it in the equality just obtained. As m is arbitrary, we can take the limit as m → ∞. By noting that v m converges even strongly to v in H 1 (0, T; H) ∩ L 2 (0, T; V σ B ), we conclude that (114) is satisfied for such a v. By similarly reasoning for the other equations, we can conclude. Hence, the existence part of the statement is proved.
Uniqueness. By linearity, we can assume that all the right-hand sides of the strong formulation (104)-(107) vanish, so that the problem becomes (q + p)(T) = 0 and r(T) = 0.
Next, we add the identities just obtained to each other. Several cancellations occur, and what remains is just the following identity: At this point, we observe that the left-hand side equals (g 4 , ζ(T)) by (116), so that the above identity becomes T 0 (g 1 , ξ) ds + g 2 , ξ(T) + Hence, by recalling the notation (117), and comparing with (103), we obtain (140).

Conclusions
In this paper, we analyzed a system of three evolution equations involving fractional operators and nonlinearities. This system turns out to be a generalization of a Cahn-Hilliard type phase field system modeling tumor growth. Actually, we considered the distributed control problem in terms of the control in the third equation, by dealing with a tracking-type cost functional. We discussed separation properties, continuous dependence with respect to the control, and well-posedness results for the state system (cf. Theorems 1-3). Then, we introduced the linearized problem and showed its solvability by Theorem 4. The Fréchet differentiability of the control-to-state operator, with the Fréchet derivative mapping the increment into the solution of the linearized problem, has been established by Theorem 5. The existence of solutions to the optimal control problem has been the subject of Theorem 6. Finally, the associated adjoint system has been introduced and its well-posedness has been proved (see Theorem-7), then the first-order necessary conditions of optimality have been derived in a simple and clear form, as stated in Theorem 8.
Funding: This research received no external funding.