1 Introduction
Abstract

In this paper, we study a distributed optimal control problem for a diffuse interface model for tumor growth. The model consists of a Cahn-Hilliard 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 ”no-friction” 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 first-order and second-order 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, Cahn-Hilliard equation, Brinkman equation, first-order necessary optimality conditions, second-order sufficient optimality conditions, uniqueness of globally optimal solutions.

MSC Classification: 35K61, 76D07, 49J20, 92C50.

\pdfstringdefDisableCommands\pdfstringdefDisableCommands

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., cell-cell 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 Darcy-type 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 volume-averaged 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 non-negative and represent the shear and the bulk viscosity, respectively. The proliferation rate , the apoptosis rate and the chemotaxis parameter are non-negative constants. By , we denote the nutrient concentration in a pre-existing 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

\endlinenomath

subject to the control constraint \linenomath

(3)
\endlinenomath

for box-restrictions 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 control-to-state operator that maps any admissible control onto a corresponding unique strong solution of the state equation (CHB). Furthermore, we show that this control-to-state operator is Lipschitz-continuous, 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 control-to-costate operator which maps any admissible control onto its corresponding adjoint state. Then, we show that this control-to-costate 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 first-order 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 second-order 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 control-to-state operator and the control-to-costate 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 second-order 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 PDE-based control problems we refer to [5] where a tumor growth model of advection-reaction-diffusion type is considered. There are various papers analyzing optimal control problems for Cahn-Hilliard equations (e.g., [10, 31]). Furthermore, control problems for the convective Cahn-Hilliard 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 Cahn-Hilliard-Navier-Stokes system. As far as control problems for Cahn-Hilliard-based 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 Cahn-Hilliard-type equation coupled to a time-dependent reaction-diffusion equation for the nutrient, where the control acts as a right-hand side in this nutrient equation. The model they considered was firstly proposed in [29] and later well-posedness 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 non-stationary, 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 Neumann-Laplace operator is positive definite and self-adjoint. In particular, by the Lax-Milgram theorem and the Poincaré inequality, the inverse operator is well-defined, 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

\endlinenomath

endowed with the standard norms.

Moreover, we make the following assumptions which we will use in the rest of this paper:
Assumptions.

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

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

  3. The non-negative function belongs to , i.e., is bounded, three times continuously differentiable and its first, second and third-order derivatives are bounded. Without loss of generality, we assume that .

  4. is the smooth double-well potential, i.e., for all .

Remark 1.

  1. In principal, it would be possible to consider more general potentials . However, since the double-well potential is the classical choice for Cahn-Hilliard-type equations (apart from singular potentials like the logarithmic or double-obstacle potential) and to avoid being too technical, we focus on the above choice for in this work.

  2. 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 control-to-state 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 so-called state) of the system (CHB). Therefore, we can define the control-to-state-operator 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 Lipschitz-continuous, Fréchet-differentiable 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)
\endlinenomath

is referred to as the set of admissible controls. Its elements are called admissible controls.

Note that this box-restricted set of admissible controls is a non-empty, bounded subset of the Hilbert space since for all , \linenomath

(6)
\endlinenomath

This means that \linenomath

(7)
\endlinenomath

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 non-negativity 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)
\endlinenomath

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)
\endlinenomath

Due to Hölder’s and Young’s inequalities and (10)-(12), the Sobolev embedding and elliptic estimates, it follows that \linenomath

(14)
\endlinenomath

Now, we observe that

Using Hölder’s, Young’s and Gagliardo-Nirenberg’s inequalities, the assumptions on , elliptic regularity theory and (10), this implies \linenomath

(15)
\endlinenomath

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 Gagliardo-Nirenberg’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

\endlinenomath

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

\endlinenomath

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)
\endlinenomath

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)
\endlinenomath

where we used that

Using Proposition 3, (24)-(25) together with Hölder’s and Young’s inequalities, it follows that \linenomath

(27)
\endlinenomath

Using the continuous embedding , Proposition 3, (24), the assumptions on and the elliptic estimate

we obtain \linenomath

(28)
\endlinenomath

Now, we observe that \linenomath

\endlinenomath

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