{'fraction_non_alphanumeric': 0.07038866333945333, 'fraction_numerical': 0.024462775723410306, 'mean_word_length': 3.5689981096408316, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 48, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal control of treatment time in a diffuse interface model of tumor growth

We consider an optimal control problem for a diffuse interface model of tumor growth. The state equations couples a Cahn-Hilliard equation and a reaction-diffusion equation, which models the growth of a tumor in the presence of a nutrient and surrounded by host tissue. The introduction of cytotoxic drugs into the system serves to eliminate the tumor cells and in this setting the concentration of the cytotoxic drugs will act as the control variable. Furthermore, we allow the objective functional to depend on a free time variable, which represents the unknown treatment time to be optimized. As a result, we obtain first order necessary optimality conditions for both the cytotoxic concentration and the treatment time.


Introduction
There has been a recent surge in the development of phase field models for tumor growth. These models aim to describe the evolution of a tumor colony surrounded by healthy tissues which experience biological mechanisms such as proliferation via nutrient consumption, apoptosis, chemotaxis and active transport of specific chemical species. For the case of a young tumor, before the development of quiescent cells, the phase field models often consist of a Cahn-Hilliard equation coupled with a reaction-diffusion equation for the nutrient [13,22,25,26,36]. One may also treat the tumor cells and the healthy cells as inertia-less fluids, and include the effects of fluid flow into the evolution of the tumor, leading to the development of a Cahn-Hilliard-Darcy system [14,15,22,44].
Current treatments for cancer include surgery, immunotherapy (strengthening the immune system), radiotherapy (using radiation to kill cancer cells) and chemotherapy (using drugs to kill cancer cells). The latter three treatments are typically conducted in cycles. A cycle is a period of treatment followed by a (longer) period of rest, so that the patient's body can build new healthy cells. The goal of these therapeutic treatments is to shrink the tumor into a more manageable size for which surgery can be applied. Further therapeutic treatments may be necessary in order to destroy the cancer cells that may remain after the surgery.
In this work, we consider an optimal control problem involving a cancer treatment with cytotoxic drugs. It is well-known that while cytotoxic drugs mainly target and damage rapidly dividing cells such as tumor cells, the drugs can also accumulate in the body and cause adverse side-effects to the immune system and various vital organs such as the kidneys and the liver. In a worst case scenario, too much cytotoxic drugs may allow tumor cells to mutate and become resistant to the treatment. Thus, from the viewpoint of the patient, the shortest treatment time in which the objectives of the chemotherapy are achieved is the most ideal. Therefore, the optimal control problem we study involves finding the optimal drug distribution and the optimal treatment time.
For T > 0, in a bounded domain Ω ⊂ R 3 with C 3 -boundary Γ, we consider the following Cahn-Hilliard model for tumor growth, Here, α is a positive constant, ϕ denotes the difference in volume fraction, where ϕ = 1 represents the tumor phase and ϕ = −1 represents the healthy tissue phase. The function µ is a chemical potential associated to ϕ, Ψ ′ (ϕ) is the derivative of a potential Ψ(ϕ) with equal minima at ϕ = ±1, σ denotes the concentration of an unspecified chemical species acting as nutrient for the tumor cells, while u denotes the concentration of cytotoxic drugs. The function h(ϕ) is an interpolation function such that h(−1) = 0 and h(1) = 1, and the parameters P, A, C, and B denote the constant proliferation rate, apoptosis rate, nutrient consumption rate, and nutrient supply rate, respectively. The positive constants A and B are related to the thickness of the interfacial layer and the surface tension, while ∂ ν f = ∇f ⋅ ν denotes the normal derivative of f where ν is the unit outward normal of Γ.
The term h(ϕ)Pσ models the proliferation of tumor cells which is proportional to the concentration of the nutrient, the term h(ϕ)A models the apoptosis of tumor cells, and Ch(ϕ)σ models the consumption of the nutrient only by the tumor cells. The term αuh(ϕ) models the elimination of the tumor cells by the cytotoxic drugs at a constant rate α. Meanwhile, σ S denotes the nutrient concentration in a pre-existing vasculature, and B(σ S − σ) models the supply of nutrient from the blood vessels if σ S > σ and the transport of nutrient away from the domain Ω if σ S < σ.
In comparison with the models of [22], we have neglected the effects of chemotaxis and active transport, but the new feature of (1.1) is the inclusion of the effects of cytotoxic drugs via the term αuh(ϕ), and in this work the function u will act as our control. For realistic applications the control u ∶ [0, T ] → [0, 1] should be spatially constant, where u = 1 represents a full dosage and u = 0 represents no dosage. However, in the subsequent analysis, we allow for spatial dependence (see Assumption 2.1 below).
For positive constants r, β u and β T , and nonnegative constants β Q , β Ω , and β S , we consider the objective functional J r given as (1.2) In particular, (1.2) can be seen as the relaxation of the following more natural objective functional Here, τ ∈ (0, T ] represents the treatment time, ϕ Q represents a desired evolution for the tumor cells while ϕ Ω represents a desired final distribution. The first two terms of J are of standard tracking type as often considered in the literature of parabolic optimal control, and the third term of J measures the size of the tumor at the end of the treatment. The fourth term penalizes large concentrations of the cytotoxic drugs, and the fifth term of J penalizes long treatment times. Let us make the following comments: 1. A large value of ϕ − ϕ Q 2 would mean that the patient suffers from the growth of the tumor, and a large value of u 2 would mean that the patient suffers from high toxicity of the drug.
2. The function ϕ Ω can be a stable configuration of the system, so that the tumor does not grow again once the treatment is completed. One can also choose ϕ Ω as a configuration which is suitable for surgery.
3. The variable τ can be regarded as the treatment time of one cycle, i.e., the amount of time the drug is applied to the patient before the period of rest, or the treatment time before surgery. 4. It is possible to replace β T τ by a more general function f (τ ) where f ∶ R ≥0 → R ≥0 is continuously differentiable and increasing. 5. We consider T ∈ (0, ∞) as a fixed maximal time in which the patient is allowed to undergo a treatment obtained from this optimal control problem.
For technical reasons highlighted below, we consider an optimal control problem with the relaxed objective functional (1.2) and the state equations (1.1). We denote the space of admissible controls as U ad (see Assumption 2.1 below) and the optimal control problem we study in this work can be expressed as follows, minimise J r (ϕ, u, τ ) subject to (1.1), u ∈ U ad , τ ∈ (0, T ). (P) The optimal control problem (P) is a problem involving a free terminal time, and we say that (u * , τ * ) is a minimizer of (P) if where the infimum is taken over triplets (φ, w, s) such that w ∈ U ad , s ∈ [0, T ] and φ solves (1.1) with datum w. In ODE constrained optimal control where the cost functional depends on the free terminal time, the necessary optimality condition can be derived with the help of the corresponding Hamiltonian function, see for instance [34,Chapter 20] and [24,32,37]. One may use the notion of Hamiltonian functional to derive the optimality condition for the free terminal time when the state equations are partial differential equations, see in particular [2,38,39] for semilinear parabolic state equations.
Below we illustrate with an example the optimality conditions obtained with the Hamiltonian from ODE theory and with the Lagrangian method for PDE-constrainted optimization, see for instance [43, §2.10] and [31, §1.6.4]. Suppose the objective functional is of the form and ϕ satisfies for example Let u * denote an optimal control with corresponding state ϕ * . The Hamiltonian H is defined as where p act as the adjoint variable to ϕ * . From the works of [2,38,39] and also from the theory of ODE-constraint optimal control, the optimality condition for the optimal time τ * is Now, let us define the Lagrangian then one obtains from formally differentiating L with respect to τ the optimality condition for τ * , which is The adjoint equation for p is a terminal time boundary value problem: Using the terminal condition for p in the expression for ∂L ∂τ (τ * , ϕ * , u * ) we see that and thus ∂L ∂τ (τ * , ϕ * , u * ) = 0 is equivalent to (1.4). That is, the optimality conditions for the free terminal time obtain from the Hamiltonian formulation and the Lagrangian formulation coincide.
Let us briefly explain the issues with the objective functional (1.3). Formally differentiating (1.3) with respect to τ , we obtain 5) and in order for the terms in (1.5) to be well-defined, we need that Furthermore, to rigorously establish the Fréchet differentiability of J with respect to τ , it turns out that we require ∂ tt ϕ * ∈ L 2 (0, T ; L 2 (Ω)), and u * , ϕ Q ∈ H 1 (0, T ; L 2 (Ω)).
Thus, the main mathematical difficulties arise from establishing high temporal regularity for the state variables. A preliminary analysis shows that it is possible to derive such regularity but only under rather strong assumptions such as ϕ 0 ∈ H 5 (Ω), σ 0 ∈ H 3 (Ω) and ∂ t u L 2 (0,T ;L 2 (Ω)) ≤ K for some fixed K > 0. The assumption on the a priori boundedness of ∂ t u is not meaningful as in applications it will be hard to verify this condition. Furthermore, using the Lagrangian method, one can compute that the terminal condition for the adjoint variable p to ϕ * is p(τ * ) = β Ω (ϕ * (τ * ) − ϕ Ω ) + β S 2 , and so we can write (1.5) more compactly as But this would mean that we require the weak formulation for the equation of ϕ * to be satisfied pointwise in [0, T ], that is, This in turn implies that we need These difficulties motivates the current study with the relaxed objective functional J r (1.2). There have been many recent contributions regarding the well-posedness and asymptotic behaviour for phase field type tumor models, see for example [7,8,9,17,20,21] for the Cahn-Hilliard variant, and [3,15,19,33,35] for the Cahn-Hilliard-Darcy variant. From the aspect of optimal control, we mention the works of [5,6,11,12,28,46] for the Cahn-Hilliard equation, [40,45,47,48] for the convective Cahn-Hilliard equation and [18,27,29,30] for the Cahn-Hilliard-Navier-Stokes system. In the context of PDE constraint optimal control for diffuse interface tumor models, we have the recent work of [10], where the objective functional (1.3) with β S = β T = 0 and no dependence of on τ is studied with state equations given by the model of [26] and the control enters the nutrient equation as a source term, similar to the term B(σ S − σ) in (1.1c). With this work we aim to provide a contribution to the theory of free terminal time optimal control in the context of diffuse interface tumor models.
Let us provide some future directions of research motivated by this study: 1. An optimal control u that is periodic in time, reflecting the cyclic nature of therapeutic treatments.

2.
A feedback mechanism taking into account the patient's response to the therapy, and the tumor's resistance to the drug.
3. Analysis and identification of stable equilibria for diffuse interface models of tumor growth.
Plan of the paper. The paper is organised as follows. In Section 2 the general assumptions are outlined and the main results are stated. The well-posedness of the state equations (1.1) is established in Section 3. The existence of a minimizer to (P) is proved in Section 4, while the unique solvability of the linearized state equations and the Fréchet differentiability of the control-to-state mapping and of the functional J r are contained in Section 5. In Section 6, the unique solvability of the adjoint equations is studied and the first order necessary optimality conditions are derived.

General assumptions and main results
Notation. For convenience, we will often use the notation L p ∶= L p (Ω) and W k,p ∶= W k,p (Ω) for any p ∈ [1, ∞], k > 0 to denote the standard Lebesgue spaces and Sobolev spaces equipped with the norms ⋅ L p and ⋅ W k,p . In the case p = 2 we use H k ∶= W k,2 and the norm ⋅ H k . Moreover, the dual space of a Banach space X will be denoted by X * , and the duality pairing between X and X * is denoted by ⟨⋅, ⋅⟩ X . The space-time cylinder Ω × (0, T ) will be denoted by Q, and we use the notation L p (Q) to denote the spaces L p (Ω × (0, T )) for 1 ≤ p ≤ ∞. Using Fubini's theorem we have the isometric isomorphism then it holds that The following Taylor's theorem with integral remainder will be used to show the Fréchet differentiability of the control-to-state mapping. For f ∈ C 2 (R) and a, x ∈ R, it holds that The Gagliardo-Nirenberg interpolation inequality in dimension d is also useful (see [16,Theorem 10.1,p. 27]): Let Ω be a bounded domain with C m boundary, and f ∈ W m,r (Ω) ∩ L q (Ω), 1 ≤ q, r ≤ ∞. For any integer j, 0 ≤ j < m, suppose there is an α ∈ R such that There exists a positive constant C depending only on Ω, m, j, q, r, and α such that We consider the following assumptions.
The initial conditions satisfy ϕ 0 ∈ H 3 with compatibility condition ∂ ν ϕ 0 = 0 on Γ, in Ω, while the target functions satisfy ϕ Q , ϕ Ω ∈ L 2 (Q). The vasculature nutrient concentration σ S satisfies 0 ≤ σ S ≤ 1 a.e. in Q. (A2) The interpolation function h ∶ R → [0, 1] is twice continuously differentiable and Lipschitz continuous (with Lipschitz constant L h ). The parameters P, A, C, and B are nonnegative constants, and α is a positive constant. (A3) The space of admissible controls is given as (A4) The potential Ψ ∶ R → R ≥0 is three times continuously differentiable and satisfies for some positive constants {k j } 5 j=0 , for all s, t ∈ R.
We point out that as Ω is a bounded domain, there exists an open set U ⊂ L 2 (Q) such that U ad ⊂ U. In this work, we consider quartic potentials Ψ for the state equations, for which the classical double-well potential Ψ(s) = 1 4 (1 − s 2 ) 2 is one example. The wellposedness of the state equations with higher polynomial growth for Ψ is also possible, see for instance the procedure in [17, Proof of Theorem 1], but we restrict our current analysis to that of quartic potentials to simplify the computations.
such that ϕ(0) = ϕ 0 , σ(0) = σ 0 , and for a.e. t ∈ (0, T ) and for all ζ ∈ H 1 , for some positive constant C not depending on (ϕ, µ, σ, u). Let (ϕ i , µ i , σ i ) i=1,2 denote two weak solutions to (1.1) satisfying (2.8) corresponding to {u i } i=1,2 with the same initial data ϕ 0 and σ 0 . Then, there exists a positive constant C cts depending only on ϕ i L ∞ (0,T ;L ∞ ) , A, B, P, A, C, α, k 5 , T and the Lipschitz constant L h , such that for all s ∈ (0, T ], (2.9) The existence of solutions to the state equations (1.1) is proved via a fixed point argument, see also [15] for a similar argument applied to a multispecies tumor model. One may also use a Galerkin approximation, which has been applied to similar systems in [7,17,20,21,19,33,35]. The key difference here are that we establish boundedness of the nutrient concentration σ, which comes from the application of a weak comparison principle. Here we also point out that the gradient ∇ϕ is continuous on the boundary up to initial time by the embedding ϕ ∈ L ∞ (0, T ; H 2 ) ∩ H 1 (0, T ; L 2 ) ⊂⊂ C 0 ([0, T ]; H β ) for β < 2 and the trace theorem. Hence, the initial condition ϕ 0 needs to fulfill the boundary conditions.
The unique solvability of the state equations (1.1) allows us to define a solution operator S given as where the triplet (ϕ, µ, σ) is the unique weak solution to (1.1) with data (ϕ 0 , σ 0 , u) over the time interval [0, T ]. We use the notation ϕ = S 1 (u) for the first component of S(u). Then, we deduce the existence of a minimizer to the (P).
Theorem 2.2 (Existence of minimizer). Under Assumption 2.1, there exists at least one minimizer (ϕ * , u * , τ * ) to (P). That is, ϕ * = S 1 (u * ) with Note that we cannot exclude the trivial cases where τ * = 0 or τ * = T . To establish the Fréchet differentiability of the solution operator with respect to the control u, we first investigate the linearized state equations. For arbitrary but fixed u ∈ U ad , let (ϕ, µ, σ) = S(u) denote the unique solution triplet to (1.1) from Theorem 2.1. For w ∈ L 2 (Q), we consider the following linearized state equations, The unique solvability of (2.10) is obtained via a Galerkin procedure.
such that for a.e t ∈ (0, T ) and for all ζ ∈ H 1 , The expectation is as follows. Let u,û ∈ U ad ⊂ U be arbitrary, with (ϕ, µ, σ) = S(u) and (φ,μ,σ) = S(û) denote the unique solution triplets to (1.1) corresponding to u and u, respectively. Denote by w ∶=û − u ∈ L 2 (Q) and let (Φ w , Ξ w , Σ w ) denote the unique solution to the linearized state equations (2.10) associated to w. We define the remainders to be then, it holds that the solution operator S is Fréchet differentiable at u, the Fréchet derivative with respect to the control u, denoted as D u S, belongs to L(L 2 (Q), Y) and satisfies With the unique solvability of the linearized state equations, we have the following result on the Fréchet differentiability of the solution operator.
Theorem 2.4 (Fréchet differentiability with respect to the control). Under Assumption 2.1, the solution operator S is Fréchet differentiable in U as a mapping from L 2 (Q) to the product Banach space That is, for anyû, u ∈ U ad ⊂ U with w =û − u ∈ L 2 (Q), there exists a positive constant C diff,u not depending onû, u and w such that
We now define a reduced functional For any u ∈ U ad ⊂ U, set w = u − u * ∈ L 2 (Q) and let (Φ w , Ξ w , Σ w ) be the unique solution to (2.10) corresponding to w, the optimal control u * and the corresponding state variables (ϕ * , µ * , σ * ). By Theorem 2.4, J is Fréchet differentiable with respect to the control with (2.14) Next, we make use of the following adjoint system to eliminate the presence of the linearized state variable Φ w in (2.14), Note that the adjoint system is supplemented with terminal conditions at the optimal treatment time τ * . We now state the unique solvability result.
The first order necessary optimality conditions for the minimizer (u * , τ * ) of Theorem 2.2 also requires the Fréchet derivative of J with respect to τ , and for this we make the additional assumption on the target functions ϕ Q and ϕ Ω .
Furthermore, we extend ϕ to negative times using the initial condition, i.e., ϕ(t) = ϕ 0 for t < 0. Theorem 2.6 (Fréchet differentiability of the reduced functional with respect to time). Under Assumptions 2.1 and 2.2, let u ∈ U ad be arbitrary with corresponding state variables S(u) = (ϕ, µ, σ). The reduced functional J (u, τ ) is Fréchet differentiable with respect to τ with The first order necessary optimality conditions to (P) for the minimizer (u * , τ * ) are given as follows.

17)
and which allows for the interpretation that the optimal control u * is the L 2 (Q)-projection of

Results on the state equations
We show the existence of weak solutions to the state equations (1.1) by means of a fixed point argument. The idea is to consider the following two auxiliary problems. Let φ be given, we define the solution mapping with homogeneous Neumann boundary condition and initial condition σ 0 . Then, we define the solution mapping M by ϕ = M(φ), where ϕ is the unique solution to with homogeneous Neumann boundary conditions and initial condition ϕ 0 . Ifφ is a fixed point for M, withσ = M 1 (φ) andμ = AΨ ′ (φ) − B∆φ, then the triplet (φ,μ,σ) is a solution to (1.1).

Auxiliary problems
to (AP1) such that σ(0) = σ 0 and 0 ≤ σ ≤ 1 a.e. in Q. Furthermore there exists a positive constant C AP1 not depending on φ such that Proof. As (AP1) is a linear parabolic equation in σ, the existence of weak solutions can be shown using a Galerkin approximation, and we will only present the derivation of a priori estimates here. The weak formulation of (AP1) is for a.e. t ∈ (0, T ) and for all ζ ∈ H 1 .
Neglecting the nonnegative term Ch(φ) σ 2 + B σ 2 , and the application of the Gronwall inequality leads to where we used the boundedness of h. Integrating in time and using that σ 0 ∈ H 1 , σ S ∈ L 2 (Q) we have Third estimate. Note that (3.2) can be seen as a weak formulation for the following elliptic problem, As the right-hand side of (3.5) belongs to L 2 for a.e. t ∈ (0, T ), elliptic regularity theory [23, Theorem 2.4.2.7] yields that σ(t) ∈ H 2 for a.e. t ∈ (0, T ), with the estimate where C is a positive constant not depending on σ and φ. Integrating in time gives The a priori estimates (3.3), (3.4) and (3.6) are sufficient to deduce the existence of a weak solution σ satisfying (3.2). The initial condition is attained by the use of the continuous embedding L 2 (0, T ; , and using weak lower semicontinuity of the norms, we obtain (3.1). We now establish the boundedness property and continuous dependence on the data φ.
As the integrand is nonnegative, we neglect the second term on the left-hand side and upon integrating yields where we used that σ 0 ≥ 0 a.e. in Ω, and so σ − (0) = 0 a.e. in Ω. Thus σ ≥ 0 a.e. in Q. On the other hand, consider ζ = (σ − 1) Using that h is nonnegative and σ S ≤ 1 a.e. in Q, so that (1 − σ S )(σ − 1) + is nonnegative, we find that the integrand is nonnegative. Thus, after integrating from 0 to t, we obtain where we used that σ 0 ≤ 1 a.e. in Ω. This implies that σ ≤ 1 a.e. in Q.
, respectively, and with the same initial condition σ 0 and nutrient supply σ S . Then the difference σ ∶= σ 1 − σ 2 satisfies for a.e. t ∈ (0, T ) and for all ζ ∈ H 1 . Substituting ζ = σ in (3.9), neglecting the nonnegative term Ch(φ 2 ) σ 2 + B σ 2 , and integrate over [0, s] for s ∈ (0, T ], we obtain where we have used the boundedness of σ 1 and the Lipschitz property of h. Applying Gronwall's inequality (2.1) yields where we used that Next, substituting ζ = ∂ t σ in (3.9) and integrate over [0, s] leads to From (3.10) we have (3.11) and so this yields Remark 3.1. The main reason we do not employ a Galerkin approximation for the state equation (1.1) is that the computations in the weak comparison principle seems not to apply to the Galerkin solutions, in particular we cannot show that the Galerkin solutions to σ is nonnegative and bounded above by 1. Indeed, for Galerkin solutions ϕ n and σ n belong to some finite dimensional subspace W n of H 1 satisfying is nonnegative. There is also a similar issue with the nonnegativity of as it is not guaranteed that the projection of a nonnegative function is nonnegative.
Due to the estimate (3.11), we can define a continuous mapping  That is, C AP2 does not depend on φ.
Proof. Let {w i } i∈N denote the eigenfunctions of the Neumann-Laplacian with corresponding eigenvalues {λ i } i∈N : Then, it is well-known that {w i } i∈N forms an orthonormal basis of L 2 and an orthogonal basis of H 1 . As constant functions are eigenfunctions, we take w 1 = 1 with λ 1 = 0. Let n ∈ N be fixed and we define W n ∶= span{w 1 , . . . , w n } as the finite dimensional space spanned by the first n eigenfunctions, with the corresponding projection operator Π n . We consider sequences {φ n } n∈N , {u n } n∈N ⊂ C 0 (0, T ; L 2 ) such that φ n → φ and u n → u strongly in L 2 (0, T ; L 2 ) and look for functions of the form where the coefficients a n ∶= {a n,i } n i=1 and b n ∶= {b n,i } n i=1 satisfy the following initial-value problem with prime denoting the time derivative and for 1 ≤ i, j ≤ n, Without loss of generality, we assume that 0 ≤ u n ≤ 1 a.e. in Q for all n ∈ N and from Lemma 3.1 it holds that 0 ≤ M 1 (φ n ) ≤ 1 a.e. in Q and M 1 (φ n ) ∈ C 0 ([0, T ]; L 2 ) for all n ∈ N. Substituting (3.15b) into (3.15a) leads to a system of ODEs in a n with right-hand side depending continuously on t and a n . By the Cauchy-Peano theorem [4, Chapter 1, Theorem 1.2], there exists a t n ∈ (0, T ] such that (3.15) has a local solution a n on [0, t n ) with a n ∈ C 1 ([0, t n ); R n ). Then, b n can be defined by the relation (3.15b), and we obtain functions ϕ n , µ n ∈ C 1 ([0, t n ); W n ) satisfying In the following we will derive a series of a priori estimates leading to the uniform boundedness (in n) of (ϕ n , µ n ) in the following Bochner spaces: 3. µ n ∈ L ∞ (0, T ; L 2 ) ∩ L 2 (0, T ; H 2 ), ϕ n ∈ L ∞ (0, T ; H 2 ), ∂ t ϕ n ∈ L 2 (0, T ; L 2 ).
In particular for the third estimate, we have to differentiate (3.17) in time to obtain a system of ODEs involving ∂ t µ n . Thus, we prescribe additional initial conditions, namely we set Note that by Assumption 2.1, there exists a positive constant C ini , not depending on φ and n, such that Furthermore, to approximate ϕ 0 by a linear combination of eigenfunctions of the Neumann-Laplacian in H 2 , we require that ϕ 0 satisfies zero Neumann boundary conditions.
First estimate. Multiplying (3.17a) with µ n and (3.17b) with ∂ t ϕ n , integrate over Ω and integrate by parts, upon adding and using the boundedness of h, M 1 (φ n ) and u n , we obtain Let C u ∶= P + A + α, then by the Poincaré inequality in L 1 (with constant C p > 0 depending only on Ω), Hölder's inequality and Young's inequality, the right-hand side of (3.18) can be estimated as follows, From integrating (3.17b) over Ω, and (2.4), we find that Hence, we obtain the following differential inequality By the Sobolev embedding H 1 ⊂ L 6 and the growth assumption (2.6), it holds that Applying the Gronwall inequality (2.1) gives where the constant C depends only on T , C u , C p , k 0 , k 1 , A, B, Ω , and ϕ 0 H 1 . From (3.20) and (3.22), the mean of µ n is bounded in L ∞ (0, T ), and the Poincaré inequality gives that µ n is bounded in L 2 (0, T ; L 2 ). Meanwhile, by (2.5) we see that and thus by (3.22), the mean of ϕ n is bounded in L ∞ (0, T ), and by the Poincaré inequality we obtain that ϕ n is also bounded in L ∞ (0, T ; L 2 ). Thus, there exists a positive constant C, not depending on φ n and n such that and as a result, this guarantees that the Galerkin solutions (ϕ n , µ n ) can be extended to the interval [0, T ], and thus t n = T for each n ∈ N.
Second estimate. From (2.6) and the Sobolev embedding H 1 ⊂ L 6 , we have that where C is a positive constant depending only on k 4 and Ω. Since Π n (Ψ ′ (ϕ n )) L 2 ≤ Ψ ′ (ϕ n ) L 2 , applying elliptic regularity to (3.17b) yields that ϕ n (t) ∈ H 2 for a.e. t ∈ (0, T ) and satisfies Then, by (2.6) we have , and so Ψ ′ (ϕ n ) ∈ L 2 (0, T ; H 1 ). Application of elliptic regularity yields that ϕ n (t) ∈ H 3 for a.e. t ∈ (0, T ) and Multiplying (3.24) with µ n and (3.17a) with B∂ t ϕ n , and then integrating over Ω, we obtain upon summing (3.25) From (2.6), we find that and so Ψ ′′ (ϕ n ) is bounded in L ∞ (0, T ; L 3 ). Applying Hölder's inequality on the right-hand side of (3.25) yields where we recall C u = P+A+α and C Sob is the positive constant from the Sobolev embedding H 1 ⊂ L 6 depending only on Ω. Then, integrating in time and using that µ n is bounded in L 2 (0, T ; H 1 ), and µ n (0) 2 Together with the improved regularity µ n ∈ L ∞ (0, T ; L 2 ), when we revisit the elliptic equation (3.17b) we find that with a positive constant C not depending on φ n and n. Similarly, viewing (3.17a) as an elliptic problem for µ n , and as ∂ t ϕ n ∈ L 2 (0, T ; L 2 ), we have by elliptic regularity where the positive constant C does not depend on φ n or n.
Compactness. From the above a priori estimates, we obtain for a relabelled subsequence, where d is the space dimension, we find that H 2 (Ω) is compactly embedded into C 0 (Ω). Hence, by [42, §8, Corollary 4] we have the following strong convergences for any 1 ≤ r < 6. The initial condition ϕ 0 is attained from that fact that ϕ ∈ C 0 ([0, T ]; H 1 ). It follows from standard arguments that the pair (ϕ, µ) satisfies (3.13), see for instance [20,21]. Furthermore, by weak lower semicontinuity of the norms, we obtain (3.14).
Continuous dependence. Let {(ϕ i , µ i )} i=1,2 denote two solution pairs satisfying (3.13) with the same initial condition ϕ 0 and corresponding data {(φ i , u i )} i=1,2 , respectively. Then, it holds that the difference ϕ ∶= ϕ 1 − ϕ 2 and µ ∶= µ 1 − µ 2 satisfy Substituting ζ = Bϕ in (3.27a) and ζ = µ in (3.27b), integrating over [0, s] for s ∈ (0, T ] and upon adding we obtain (3.28) By the boundedness of M 1 (φ 1 ) and u 1 , the Lipschitz continuity of h, we obtain while by Hölder's inequality and Young's inequality, and the boundedness of h, we have Using (2.7) and the fact that ϕ i ∈ C 0 (Q), we find that  Then, substituting the above three estimates into (3.28) we obtain for s ∈ (0, T ], is a positive constant. Applying (2.1) yields for any s ∈ (0, T ], where we used that W (t) ∶= u 2 L 2 (0,t;L 2 ) + M 1 2 L 2 (0,t;L 2 ) is a nondecreasing function of t, and thus We point out that although the source term in (1.1a) closely resembles that of [20], we obtain a priori estimates for potentials Ψ with quartic growth (see (2.6)), which is in contrast to the quadratic potentials considered in [20]. The main difference is that here we have the boundedness of the nutrient, and thus we only require a bound on the mean of µ in (3.19). But in [20], the presence of the active transport mechanism (modeled by the term div (n(ϕ)χ∇ϕ) in the nutrient equation) prevents us from applying a weak comparison principle to deduce the boundedness of the nutrient. Without the boundedness of the nutrient, we have to control the square of the mean of µ in order to estimate the source term in (3.19).

Existence by Schauder's fixed point theorem
Note that if {φ n } n∈N is a bounded sequence in L 2 (Q), by Lemma 3.1 the corresponding sequence {σ n ∶= M 1 (φ n )} n∈N satisfies 0 ≤ σ n ≤ 1 a.e. in Q, and by Lemma 3.2 we have that the corresponding solution pair {ϕ n , µ n } n∈N is bounded uniformly in which yields a strongly convergent (relabelled) subsequence {ϕ n } n∈N in L 2 (Q), due to the compact embedding The problem ϕ = λM(ϕ) translates to By Lemma 3.1 we have that 0 ≤ σ ≤ 1 a.e. in Q for all λ ∈ [0, 1], and thus we can choose M to be the constant C AP2 in (3.14) which does not depend on ϕ and λ ∈ [0, 1]. Thus Schauder's fixed point theorem yields the existence of a weak solution (ϕ, µ, σ) to the state equations (1.1) with 0 ≤ σ ≤ 1 a.e. in Q and for some positive constant C not depending on (ϕ, µ, σ, u).

Substituting this into (3.30) leads to
we obtain from (2.1) that Combining with (3.10) and (3.12), we find that there exists a positive constant C 1 , depending only on B, Q, C, L h , T such that for s ∈ (0, T ]. Next, we find using (2.7) and the fact that ϕ i ∈ C 0 (Q) for i = 1, 2, , and so viewing (3.27b) as an elliptic problem for ϕ, we obtain by elliptic regularity where C 2 is a positive constant depending only on Ω, A, k 5 , ϕ i L ∞ (Q) , T and C 1 .

Existence of a minimizer
From (3.31) it holds that where C is a positive constant independent of (ϕ, µ, σ, u). Hence, we obtain that As J r is bounded from below, we can consider a minimising sequence (u n , τ n ) n∈N with u n ∈ U ad , τ n ∈ (0, T ) and corresponding weak solutions (ϕ n , µ n , σ n ) n∈N on the interval [0, T ] with ϕ n (0) = ϕ 0 and σ n (0) = σ 0 for all n ∈ N, such that In particular, u n ∈ U ad implies that 0 ≤ u n ≤ 1 a.e. in Q for all n ∈ N. As {τ n } n∈N is a bounded sequence, there exists a relabelled subsequence such that where (ϕ * , µ * , σ * , u * ) satisfy (2.8) with 0 ≤ u * , σ * ≤ 1 a.e. in Q. Note that by the dominated convergence theorem, for all p ∈ [1, ∞), Then, by the strong convergence of ϕ n −ϕ Q to ϕ * −ϕ Q in L 2 (Q) and the strong convergence  Next, using the weak lower semicontinuity of the L 2 (Q)-norm, we have and so, passing to the limit n → ∞ in the above inequality and using the convergence τ n → τ * , we see that Then, by passing to the limit n → ∞ in J r (ϕ n , u n , τ n ) and using (4.1), (4.2), and (4.3), we have which implies that (u * , τ * ) is a minimizer of (P).

Unique solavability of the linearized state equations
Recalling the set {w i } i∈N of eigenfunctions of the Neumann-Laplacian from the proof of Lemma 3.2, we look for functions of the form for all v ∈ W n . Substituting v = w j leads to γ ′ n = −Sδ n − M h n Pη n + J n − K n γ n , (5.2a) where the matrix S has been defined in (3.16), and for 1 ≤ i, j ≤ n, Taking an approximating sequence in C 0 ([0, T ]; L 2 ) for u, which we will abuse notation and reuse the variable u, and then supplementing (5.2) with the initial conditions γ n (0) = 0 and η n (0) = 0 leads to a system of ODEs with right-hand sides depending continuously on (t, γ n , η n ). Thus, by the Cauchy-Peano theorem, there exists t n ∈ (0, T ] such that (5.2) has a local solution (γ n , δ n , η n ) on [0, t n ] with γ n , δ n , η n ∈ C 1 ([0, t n ); R n ). Then, we obtain functions Φ n , Ξ n , Σ n ∈ C 1 ([0, t n ); W n ) satisfying (5.1).

First estimate.
Substituting v = Φ n in (5.1a), v = ∆Φ n in (5.1b) and v = Σ n in (5.1c), integrating over [0, t] for t ∈ (0, T ], and integrating by parts, we obtain after summation where we used that Σ n (0) = Φ n (0) = 0 and have neglected the nonnegative term (B + Ch(ϕ)) Σ n 2 . From Theorem 2.1, we have ϕ ∈ C 0 (Q), and as Ψ ′′ , Ψ ′′′ , h ′ and h ′′ are continuous with respect to their arguments, it holds that there exists a constant C * > 0 such that Then, applying Hölder's inequality and Young's inequality we obtain Using the estimates for I 1 , I 2 , I 3 and I 4 , we obtain where C 5 , C 6 > 0 are positive constants depending only on C * , A, B, P, A, C, and α.
Applying the integral form of Gronwall's inequality we obtain that for some constant D 1 not depending on n, which in turn implies that {Φ n } n∈N is bounded uniformly in L ∞ (0, T ; L 2 ), Second estimate. Substituting v = ∂ t Σ n in (5.1c), we obtain Applying Gronwall's inequality yields that where D 2 is a positive constant not depending on n. Hence, {Σ n } n∈N is bounded uniformly in L ∞ (0, T ; H 1 ) ∩ H 1 (0, T ; L 2 ).

Fréchet differentiability with respect to the control
In this section, we use the notation ϕ w =φ, µ w =μ, σ w =σ. The remainders (θ w , ρ w , ξ w ) from (2.12) satisfy for a.e. t ∈ (0, T ) and for all ζ ∈ H 1 with Using the Taylor's theorem with integral remainder (2.2) we see that Thanks to the fact that ϕ, ϕ w ∈ C 0 (Q) and the continuity of Ψ ′′′ and h ′′ , we see that there exists a constant C * * > 0 such that Furthermore, we can express Let X w ∶= Pσ w − A − α(u + w) and X ∶= Pσ − A − αu. Then, it holds similarly that and thus, we see that (θ w , ρ w , ξ w ) satisfy for a.e. t ∈ (0, T ) and for all ζ ∈ H 1 .
First estimate. Let us first compute the following preliminary estimates, using the continuous dependence estimate (2.9), the Lipschitz continuity of h, Hölder's inequality, Young's inequality and the embedding L 2 (0, T ; where C 10 is a positive constant depending only on C cts , C Sob and L h . Meanwhile, using the boundedness of σ, h ′ (ϕ) and R w 2 in Q, we see that Thus, when we substitute ζ = ξ w in (5.15c), integrating over [0, s] for s ∈ (0, T ], and neglecting the nonnegative term B ξ w 2 + Ch(ϕ) ξ w 2 , we obtain Next, substituting ζ = θ w in (5.15a), ζ = θ w in (5.15b) and ζ = 1 B ρ w in (5.15b), integrating by parts and integrating over [0, s] for s ∈ (0, T ], and upon adding leads to (5.17) Using (2.9), Hölder's inequality, Young's inequality, the boundedness of Ψ ′′ (ϕ) and R w 1 in Q, we have where C 11 is a positive constant depending only on C * , C * * , C cts , A and B. Meanwhile, by the Lipschitz continuity of h, and the fact that we see that for some positive constant C 12 depending only on B, L h , P and α. Furthermore, using the boundedness of X, h ′ (ϕ), R w 2 and h(ϕ) in Q, we have K 3 ≤ (P + A + α)C * θ w 2 L 2 (0,s;L 2 ) + P ξ w L 2 (0,s;L 2 ) θ w L 2 (0,s;L 2 ) where we recall C u = (P + A + α). Substituting the above estimates into (5.17) we obtain where the positive constants C 13 , C 14 depend only on C, C cts , C * , C * * , C 10 , C 11 , C 12 , P, A, α, and B. Applying Gronwall's inequality to (5.19) we have that Second estimate. Substituting ζ = ∂ t ξ w in (5.15c), integrating over [0, s] for s ∈ (0, T ] leads to Using the Lipschitz continuity of h, the boundedness of σ, h(ϕ), h ′ (ϕ), and R w 2 in Q, Hölder's inequality, Young's inequality and (5.20), we obtain for some positive constant C 16 depending only on T , L h , C, C cts , C * , C * * and C 15 .
6 First order necessary optimality conditions 6.1 Unique solvability of the adjoint system We apply a Galerkin approximation and consider a basis {w i } i∈N of H 2 that is orthonormal in L 2 , and we look for functions of the form for all v ∈ W n ∶= span{w 1 , . . . , w n }. Substituting v = w j leads to where S is defined in (3.16), and and we supplement the above backward-in-time system of ODEs with the condition r n (τ * ) = 0, p n (τ * ) = 0.
Once again, we consider approximating sequences in C 0 ([0, T ]; L 2 ) for u, ϕ Q and ϕ Ω and use the same variables to denote the approximating functions. Note that the right-hand side of (6.2) depends continuously on (P n , Q n , R n ) but due to the term χ (τ * −r,τ * ) (t)G n in the equation for P ′ n , we cannot apply the Cauchy-Peano theorem directly. But we can consider first solving (6.2) on the interval (τ * − r, τ * ], that is, P n and R n satisfy for t ∈ (τ * − r, τ * ], which would yield, via the Cauchy-Peano theorem, the existence of t n ∈ [τ * − r, τ * ) and a local solution pair (P n , R n ) ∈ C 1 ((t n , τ * ]; R n ) 2 to (6.3). The a priori estimates derived below will allow us to deduce that (P n , R n ) can be extended to τ * − r, that is, t n = τ * − r for all n ∈ N. Then, we then extend the solutions by solving the system with terminal conditions at time τ * − r. Overall, this procedure yields functions p n , q n , r n ∈ C 1 ((t n , τ * ]; W n ) satisfying (6.1) for some t n ∈ [0, τ * ). We now derive the a priori estimates.
Second estimate. Viewing (6.1b) as the weak formulation of an elliptic problem for p n , and using that q n is bounded uniformly in L 2 (0, τ * ; L 2 ), we have by elliptic regularity that p n L 2 (0,τ * ;H 2 ) ≤ C q n L 2 (0,τ * ;L 2 ) + p n L 2 (0,τ * ;L 2 ) , for some positive constant C not depending on n.