Abstract
In this paper, we study a distributed optimal control problem for a diffuse interface model for tumor growth. The model consists of a CahnHilliard type equation for the phase field variable coupled to a reaction diffusion equation for the nutrient and a Brinkman type equation for the velocity. The system is equipped with homogeneous Neumann boundary conditions for the tumor variable, the chemical potential and the nutrient as well as a ”nofriction” boundary condition for the velocity. The control acts as a medication by cytotoxic drugs and enters the phase field equation. The cost functional is of standard tracking type and is designed to track the phase field variable during the evolution and at some fixed final time. We prove that the model satisfies the basics for calculus of variations and we establish firstorder and secondorder conditions for local optimality. Moreover, we present a globality condition for critical controls and we show that the optimal control is unique on small time intervals.
Keywords: Optimal control with PDEs, calculus of variations, tumor growth, CahnHilliard equation, Brinkman equation, firstorder necessary optimality conditions, secondorder sufficient optimality conditions, uniqueness of globally optimal solutions.
MSC Classification: 35K61, 76D07, 49J20, 92C50.
Optimal control theory and advanced optimality conditions
for a diffuse interface model of tumor growth
Matthias Ebenbeck and Patrik Knopf
[1mm] Department of Mathematics, University of Regensburg, 93053 Regensburg, Germany
[1mm] Matthias.Ebenbeck@ur.de, Patrik.Knopf@ur.de
[5mm] Please cite as:
1 Introduction
The evolution of cancer cells is influenced by a variety of biological mechanisms like, e.g., cellcell adhesion or mechanical stresses, see [4]. Although cancer is one of the most common causes of death, the knowledge about underlying processes is still at an unsatisfying level. Due to the complexity of the evolution, experimental and medical research may not be sufficient to predict growth and to establish general treatment strategies. Thus, it is of high importance to develop biologically realistic and predictive mathematical models to identify the influence of different growth factors and mechanisms and to set up individual treatment.
In the recent past, diffuse interface models gained much interest (e.g., [21, 26, 29, 30]) and some of them seem to compare well with clinical data, see [1, 4, 19]. These models are typically based on mass and momentum balance equations and mass/momentum exchange between the different phases and can be closed by appropriate constitutive laws. Several biological mechanisms like chemotaxis, mitosis, angiogenesis or necrosis can be incorporated or effects due to stress, plasticity or viscoelasticity can be included, see [13, 25, 36].
In general, living biological tissues behave like viscoelastic fluids, see [9]. Since relaxation times of elastic materials are rather short (see [22]), it is reasonable to consider Stokes flow as an approximation of those tissues, see e.g., [17]. Indeed, many authors used Stokes flow to describe the tumour as a viscous fluid, see [18, 20]. In classical tumor growth models, velocities are modelled with the help of Darcy’s law. In these models the velocity is assumed to be proportional to the pressure gradient caused by the birth of new cells and by the deformation of the tissue, see [7, 28]. Brinkman’s law is an interpolation between the viscous fluid and the Darcytype models, see e.g., [38, 46].
Introduction of the model.
In this paper, we consider the following model: Let with be a bounded domain. For a fixed final time , we write . By we denote the outer unit normal on and denotes the outward normal derivative of the function on . Our state system is given by
(1a)  
(1b)  
(1c)  
(1d)  
(1e)  
(1f)  
(1g)  
(1h) 
where the viscous stress tensor is defined by
(2) 
and the symmetric velocity gradient is given by
In (1), denotes the volumeaveraged velocity of the mixture, denotes the pressure and denotes the concentration of an unknown species acting as a nutrient. The function denotes the difference of the local relative concentrations of tumor tissue and healthy tissue where represents the region of pure tumor tissue, stands for the surrounding pure healthy tissue and represents the transition between these pure phases. Furthermore, denotes the chemical potential associated with . The positive constant represents the mobility for the phase variable . The thickness of the diffuse interface is modelled by a small parameter and the constant is related to the fluid permeability. Moreover, the constants and are nonnegative and represent the shear and the bulk viscosity, respectively. The proliferation rate , the apoptosis rate and the chemotaxis parameter are nonnegative constants. By , we denote the nutrient concentration in a preexisting vasculature and is a positive constant. Hence, the term models the nutrient supply from the blood vessels if and the nutrient transport away from the domain for . The term in (1c) models the elimination of tumor cells by cytotoxic drugs and the function will act as our control.
Since it does not play any role in the analysis, we set .
We investigate the following distributed optimal control problem:
\linenomath
subject to the control constraint \linenomath
(3) 
for boxrestrictions and the state system (CHB). Here, and are nonnegative constants. The optimal control problem can be interpreted as the search for a strategy how to supply a medication based on cytotoxic drugs such that

a desired evolution of the tumor cells is realized as good as possible;

a therapeutic target (expressed by the final distribution ) is achieved in the best possible way;

the amount of supplied drugs does not cause harm to the patient (expressed by both the control constraint and the last term in the cost functional).
The ratio between the parameters , and can be adjusted according to the importance of the individual therapeutic targets.
In the case when , the term models the elimination of tumor cells by a supply of cytotoxic drugs represented by the control . This specific control term has been investigated in [16] and also in [24] where a simpler model was studied in which the influence of the velocity is neglected. However, in some situations it may be more reasonable to control, for instance, the evolution at the interface and one has to use a different form for , see Remark 1 below. Therefore, we allow to be rather general.
Summary of our main results.
In Section 3, we prove the existence of a controltostate operator that maps any admissible control onto a corresponding unique strong solution of the state equation (CHB). Furthermore, we show that this controltostate operator is Lipschitzcontinuous, Fréchet differentiable and satisfies a weak compactness property. In particular, we establish the fundamental requirements for calculus of variations.
In Section 4, we investigate the adjoint system. Its solution, that is called the adjoint state or the costate, is an important tool in optimal control theory as it provides a better description of optimality conditions. We prove the existence of a controltocostate operator which maps any admissible control onto its corresponding adjoint state. Then, we show that this controltocostate operator is Lipschitz continuous and Fréchet differentiable.
Eventually, in Section 5, we investigate the above optimal control problem. First, we show that there exists at least one globally optimal solution. After that, we establish firstorder necessary conditions for local optimality. These conditions are of great importance for possible numerical implementations as they provide the foundation for many computational optimization methods. We also present a secondorder sufficient condition for strict local optimality, a globality criterion for critical controls and a uniqueness result for the optimal control on small time intervals.
Comparison with the results established in [16].
We want to outline the main features and novelties of our work, especially compared to the results in [16]. The main difference lies in the equation for the nutrient concentration. In [16] the equation
(4) 
was used where denotes some positive boundary permeability constant. However, in this paper, we use equation (1e) endowed with a homogeneous Neumann boundary condition (as studied in [24]) to describe the nutrient distribution.
Apart from the fact that (1e) might be more reasonable in some situations from a modeling point of view, it provides some advantages with regard to analysis and optimal control theory. In particular, Lipschitz continuity and Fréchet differentiability of both the controltostate operator and the controltocostate operator can be established in much better function spaces. As a consequence, we can prove that the cost functional is twice continuously differentiable which was not possible in [16]. This enables us to establish secondorder conditions for (strict) local optimality.
Moreover, as mentioned above, we present a globality criterion for critical controls as well as a uniqueness result for the optimal control on small time intervals. These new conditions have neither been established in [16] nor for related models in literature. However, we are convinced that analogous conditions could also be proved for the optimal control problem in [16].
Related results in literature.
Finally, we want to mention further works where optimal control problems for tumor models are studied.
Results on optimal control problems for tumor models based on ODEs are investigated in [35, 37, 39, 42]. In the context of PDEbased control problems we refer to [5] where a tumor growth model of advectionreactiondiffusion type is considered.
There are various papers analyzing optimal control problems for CahnHilliard equations (e.g., [10, 31]). Furthermore, control problems for the convective CahnHilliard equation where the control acts as a velocity were investigated in [12, 27, 44, 45] whereas in [6] the control enters in the momentum equation of a CahnHilliardNavierStokes system. As far as control problems for CahnHilliardbased models for tumor growth are considered, there are only a few contributions where an equation for the nutrient is included in the system. In [11], they investigated an optimal control problem consisting of a CahnHilliardtype equation coupled to a timedependent reactiondiffusion equation for the nutrient, where the control acts as a righthand side in this nutrient equation. The model they considered was firstly proposed in [29] and later wellposedness and existence of strong solutions were established in [21]. However, effects due to velocity are not included in their model and mass conservation holds for the sum of tumor and nutrient concentrations. Lastly, we want to cite the paper [24] about an optimal control problem of treatment time where the control represents a medication of cytotoxic drugs and enters the phase field equation in the same way as ours. Although their nutrient equation is nonstationary, some of the major difficulties do not occur since the velocity is assumed to be negligible ().
2 Preliminaries
We first want to fix some notation: For any (real) Banach space , its corresponding norm is denoted by . denotes the dual space of and stands for the duality pairing between and . If is an inner product space, its inner product is denoted by . We define the scalar product of two matrices by
For the standard Lebesgue and Sobolev spaces with , , we use the notation and with norms and respectively. In the case we write and the norm . Sometimes, we will also write instead of . By , and , we denote the corresponding spaces of vector or matrix valued functions. For Bochner spaces, we use the notation for a Banach space and . For the dual space of a Banach space , we introduce the (generalized) mean value by
Moreover, we introduce the function spaces
Then, the NeumannLaplace operator is positive definite and selfadjoint. In particular, by the LaxMilgram theorem and the Poincaré inequality, the inverse operator is welldefined, and we set for if and
We have dense and continuous embeddings and the identifications , for all and .
Furthermore, we define the function spaces
\linenomath
endowed with the standard norms.
Moreover, we make the following assumptions which we will use in the rest of this paper:
Assumptions.

The domain is bounded with boundary , the initial datum and are prescribed.

The constants , , , , are positive and the constants , , , are nonnegative.

The nonnegative function belongs to , i.e., is bounded, three times continuously differentiable and its first, second and thirdorder derivatives are bounded. Without loss of generality, we assume that .

is the smooth doublewell potential, i.e., for all .
Remark 1.

In principal, it would be possible to consider more general potentials . However, since the doublewell potential is the classical choice for CahnHilliardtype equations (apart from singular potentials like the logarithmic or doubleobstacle potential) and to avoid being too technical, we focus on the above choice for in this work.

For the function , there are two choices which are quite popular in the literature. In, e.g., [23, 26], the choice for is given by
satisfying , . Other authors preferred to assume that is only active on the interface, i.e., for values of between and , which motivates functions of the form
see, e.g., [30, 32]. Surely, we would have to use regularized versions of these choices to fulfill (A3).
3 The controltostate operator and its properties
In this section, we consider the equation (CHB) as presented in the introduction. First, we define a certain set of admissible controls that are suitable for the later approach. We will see that each of these admissible controls induces a unique strong solution (the socalled state) of the system (CHB). Therefore, we can define the controltostateoperator mapping any admissible control onto its corresponding state. We show that this operator has some important properties that are essential for calculus of variations: It is Lipschitzcontinuous, Fréchetdifferentiable and satisfies a weak compactness property.
3.1 The set of admissible controls
The set of admissible controls is defined as follows:
Definition 2.
Let be arbitrary fixed functions with almost everywhere in . Then the set \linenomath
(5) 
is referred to as the set of admissible controls. Its elements are called admissible controls.
Note that this boxrestricted set of admissible controls is a nonempty, bounded subset of the Hilbert space since for all , \linenomath
(6) 
This means that \linenomath
(7) 
Obviously, the set is also convex and closed in . Therefore, it is weakly sequentially compact (see [43, Thm. 2.11]).
3.2 Strong solutions and uniform bounds
We can show that the system (CHB) has a unique strong solution for every control :
Proposition 3.
Let be arbitrary. Then, there exists a strong solution quintuple of (CHB). Moreover, every strong solution quintuple satisfies the following bounds that are uniform in :
(8) 
where is a constant that depends only on the system parameters and on , and .
Proof.
The assertion follows with slight modifications in the proof of [16, Thm. 4]. We will sketch the main differences in the following:
Step 1:
Testing (1e) with , using (1f), the nonnegativity of and Hölder’s and Young’s inequalities, we obtain
meaning
(9) 
Then, with exactly the same arguments as in [16, Proof of Thm. 4] it follows that \linenomath
(10) 
with a constant depending only on the system paramters and on , and .
Step 2:
Now, we want to establish higher order estimates. Using elliptic regularity theory, the assumptions on and , (1e)(1g), (9) and (10), it is easy to check that
(11) 
Together with the boundedness of and the Sobolev embedding , this implies
(12) 
Testing (1c) with , (1c) with , integrating by parts and summing the resulting identities, we obtain \linenomath
(13) 
Due to Hölder’s and Young’s inequalities and (10)(12), the Sobolev embedding and elliptic estimates, it follows that \linenomath
(14) 
Now, we observe that
Using Hölder’s, Young’s and GagliardoNirenberg’s inequalities, the assumptions on , elliptic regularity theory and (10), this implies \linenomath
(15) 
Using (14)(15) in (13), recalling (10) and using elliptic regularity theory, a Gronwall argument yields
(16) 
Then, a comparison argument in (1d) yields
(17) 
A further comparison argument in (1c) yields
(18) 
Using (10)(12), (16)(17), the assumptions on and GagliardoNirenberg’s inequality, it is easy to check that
Due to [14, Lemma 1.5], this implies
(19) 
which completes the proof. ∎
From interpolation theory, we can conclude that the component of a strong solution quintuple has a representative that is continuous on .
Corollary 4.
Let and be arbitrary and let denote the strong solution of the system (CHB). Then, has the following additional properties: \linenomath
for some constant that depends only on the system parameters and on , , and .
Proof.
The assertion can be established by an interpolation argument. The proof proceeds completely analogously to the proof of [16, Cor. 5]. ∎
Furthermore, we can show that any control induces a unique strong solution of the system (CHB):
Theorem 5.
Let and be arbitrary and let denote the corresponding strong solution as given by Proposition 3. Then, this strong solution is unique.
Proof.
Let be arbitrary and let denote a generic nonnegative constant that depends only on , , and and may change its value from line to line. For brevity, we set \linenomath
where and are strong solutions of (CHB) to the controls and . In particular, this means that both strong solutions satisfy the initial condition (1h), i.e., holds almost everywhere in .
Then, the following equations are satisfied:
(20a)  
(20b)  
(20c)  
(20d)  
(20e)  
(20f)  
(20g)  
(20h) 
Now, we will show that if . The argumentation is split into two steps:
Step 1:
In [16], it has been shown that the following inequalities hold: For any and all , \linenomath
(21)  
(22) 
Now, multiplying (20e) with , integrating by parts and using (20f), it follows that
Using the assumptions on , Proposition 3 and Hölder’s and Young’s inequalities, it is therefore easy to check that
(23) 
Then, we can follow the arguments in [15, Sec. 5] to deduce that
(24) 
Step 2:
We now prove higher order estimates. Using elliptic regularity theory, Proposition 3, (23)(24) and the assumptions on , it is easy to check that
(25) 
Multiplying (20c) with and inserting the expression for given by (20d), we obtain \linenomath
(26) 
where we used that
Using Proposition 3, (24)(25) together with Hölder’s and Young’s inequalities, it follows that \linenomath
(27) 
Using the continuous embedding , Proposition 3, (24), the assumptions on and the elliptic estimate
we obtain \linenomath
(28) 
Now, we observe that \linenomath
Due to the assumptions on and because of Proposition 3, it is straightforward to check that
where we used the continuous embedding and elliptic regularity theory. With similar argument, using the Sobolev embedding and the assumptions on , we obtain
From the last two inequalities, we obtain
(29) 
Therefore, we have
(30) 
Plugging in (27)(30) into (26), we obtain \linenomath