Mathematical analysis of a phase-field model of brain cancers with chemotherapy and antiangiogenic therapy effects

1.   Introduction

There has been recently a strong interest in the development, mathematical study and numerical analysis of phase-field models for tumor growth. Such models describe the evolution of a tumor surrounded by healthy tissues and take into account mechanisms such as proliferation of cells via nutrient consumption, therapies, clustering effects in brain tumors, etc. We refer the interested reader to, e.g., ^{[1,3,4,5,8,9,10,11,12,18,21,22,23]} for more details.

In this paper we address a phase-field model of tumor growth driven by a vital nutrient and subject to medical treatment. The model takes into account both the effects of a cytotoxic drug inhibiting the tumor proliferation rate, and those of an antiangiogenic therapy which reduces the nutrient supply. Indeed, since one of the hallmarks of cancer is microvascularization and it is more pronounced in certain tumors (such as gliomas), chemotherapy is often supported by an attendant antiangiogenic drug (usually bevacizumab, see e.g. ^[7]). This interesting approach has been introduced in the mathematical literature only recently by Colli et al. in ^[3,4] in the context of prostatic cancer.

Therein, the authors propose a model of PCa growth and chemotherapy based on ^[26], and then perform a complete mathematical analysis of the resulting system. This is based on three parabolic equations, each ruling the evolution of a characteristic variable. The first is a phase-field variable $u$ that identifies the spatial location and geometry of the tumor. In particular, the healthy tissue corresponds to tumor absence, that is $u = 0$, while $u = 1$ in the tumoral case. The second variable in play is the concentration of a vital nutrient denoted by $\sigma$, while the third variable $p$ identifies the prostate-specific antigen released by cancerous prostatic cells. In the present paper, aiming to model different tumors, we will not consider the last equation. Instead, we will focus on the equation for $\sigma$ and we will modify it in a form that seems more suitable to describe the evolution of oxygen - one of the main nutrients for brain cancers like glioma (see ^[13]).

Here is a description of the model. We assume that the evolution of $u$ is governed by a nonconserved phase-field equation, as justified by Xu, Vilanova and Gomez in ^[26] by using the concept of tumor free energy and gradient dynamics. It reads as follows

where $\lambda$ is the diffusion coefficient of tumor cells, and $\Psi = \Psi({u}, \sigma, {{\textsf{c}}})$ is a double-well potential having local minima at $u = 0$ and $u = 1$. Note that the so-called chemical free energy $\Psi$ here also depends on the administered cytotoxic drug ${{\textsf{c}}}$. We will detail later the choice of the functional $\Psi$, explaining in particular how nutrient and chemotherapy influence tumor proliferation rate.

Concerning the nutrient dynamics, the starting point is the reaction diffusion equation proposed in ^[3,4]

where $\gamma_c$ and $\gamma _h$ are the nutrient uptake rate in cancerous and healthy tissue, respectively. Analogously $S_c$ is the nutrient supply rate in the cancerous tissue, while $S_h$ refers to nutrient supply in the healthy tissue. Besides, the model incorporates the action of an antiangiogenic treatment, via the control function ${{\textsf{s}}}$ providing a reduction of the intratumoral nutrient supply rate.

Indeed, since in this paper we take oxygen as nutrient, we modify the equation by considering a nonlinear term of the form

that accounts for the oxygen uptake by cells, assuming Michaelis–Menten kinetics (and setting some biological constants equal to $1$), see e.g. ^[13].

Accordingly, we propose the following equation for the nutrient dynamics

i.e.,

having set $\gamma_{ch} = \gamma _c-\gamma _h$.

We stress that the nonlinear term $g(\sigma)$ in the equation of the nutrient is in particular relevant in brain cancers and, more precisely, in gliomas (see ^[13]). Note that another important nutrient for brain cancers development is lactate (see ^[2,17]; see also ^[19,20]); in that case, one also has Michaelis–Menten kinetics, leading to a similar equation for lactate. We will address this model in a forthcoming paper.

Our aim in this work is twofold.

● First task: Well posedness of the model. Assuming that the controls $({{\textsf{c}}}, {{\textsf{s}}})$ are given functions as above, in the first part of the paper we shall provide an existence and uniqueness result for the proposed model, see (1.2) below, after setting the proper mathematical framework.

This will be addressed in Sections 2-4; Section 1 is devoted to the precise mathematical setting.

● Second task: Optimal control. We study the problem of finding the controls $({{\textsf{c}}}, {{\textsf{s}}})$ that provide the optimal cytotoxic and antiangiogenic effects to treat a certain glioma described by (1.2).

From the mathematical viewpoint, the problem consists in

minimizing a certain cost function $J({{{c}}}, {{{s}}})$ subject to system $(1.2)$, in a prescribed class of admissible controls $({{{c}}}, {{{s}}})\in \mathcal U_{ad}$.

This will be accomplished in Sections 5-6, exploiting classical tools of control theory suitably formulated in Appendix 7.

1.1.   The mathematical model

Let $\Omega$ be a bounded and regular domain of ${\mathbb{R}}^N$ with $N = 1, 2, 3$ being the spatial dimension. Let $T > 0$ be a finite time. By defining the space/time sets $Q_T : = \Omega \times (0, T)$ and $\Sigma_T = \partial\Omega\times (0, T)$ the model can be formulated as

where ${u}_0$ and $\sigma_0$ are sufficiently smooth given functions defined in $\Omega$. We recall that $\gamma_h, \, S_h, \, S_c$ are positive constants related with the biological mechanisms in the glioma, while $\gamma_{ch} = \gamma _c-\gamma$. Besides,

The given functions $({{\textsf{c}}}, {{\textsf{s}}})$ correspond respectively to a cytotoxic drug administered as chemotherapy and to an antiangiogenic therapy; they are supposed to be positive and bounded. We further suppose that

Let us now detail the choice of the nonlinearities $F, m, h$ appearing in the first equation.

1.2.   The chemical potential $\Psi$

As derived in ^[26], the phenomenological equation for tumor phase-field relies on the choice of the tumor chemical free energy functional. This is defined as

where $F({u}) = {u}^2(1-{u})^2$ is the prototype double-well potential, the wells being ${u} = 0$ and ${u} = 1$. It is perturbed by the term $m(\sigma) h({u})$, where

is an interpolation function with the property $h(0) = 0$, $h(1) = 1$ and $h'(0) = h'(1) = 0$, while the tilting function $m$ accounts for the effects of hypoxia. Indeed, a possible choice of $m$ is

where $\sigma^{h-v}$ is the threshold oxygen concentration for hypoxia. According to ^[26], one should assume that $| m(\sigma)| < 1/3$ so that the perturbation does not destroy the double-well structure of $F$. As a consequence, the free energy functional $\Psi$ has a local maximum in $(0, 1)$ and preserves two local minima in $u = 0$ and $u = 1$ for any oxygen concentration $\sigma > 0$.

To understand how hypoxia influences proliferation, observe that, when $\sigma$ is below the threshold $\sigma^{h-v}$, the function $m(\sigma)$ is negative. Thus, since $\Psi(0, \sigma) = 0 < \Psi(1, \sigma) = - m(\sigma)$, the preferable energy level corresponds to healthy tissue. In turn, without hypoxia, $\sigma > \sigma^{h-v}$, we have $\Psi(0, \sigma) = 0 > \Psi(1, \sigma) = -m(\sigma)$, so the phase-field equation will favor tumor growth.

At this point, aiming at incorporating in the equation the tumor-inhibiting effect of a cytotoxic drug ${{\textsf{c}}}$, we define

Accordingly, the phase-field equation (1.1) reads

that can also be written as

We refer the reader to ^[6,13] for alternative models of hypoxia effects on tumor growth.

Remark 1.1. For further use, we observe that, when $\sigma\geq 0$, the Michaelis-Menten nonlinearity $g$ satisfies $0\leq g(\sigma) < 1$ and $0 < g'(\sigma) < 1$. Besides,

Furthermore, as tilting term $m$ we consider any function satisfying

Without loss of generality, we set $\lambda = 1$.

1.3.   The cost functional

In view of the optimal control problem, we consider a cost function that is based on prescribed target functions for the tumor volume and the oxygen, respectively, on $Q_T$ and on $\Omega$ at the final time $T$ of the pharmacological treatment. Accordingly, for assigned ${u}_Q, \sigma_Q\in L^2(Q_T)$, ${u}_\Omega, \sigma_\Omega\in L^2(\Omega)$, we define

where $(u, \sigma)$ is the (unique) solution to (1.2) originated by any observed initial state $(u_0, \sigma_0)$ of the system. Besides, the set of admissible controls will be

where the given quantities $U_{\rm max} > 0$ and $0 < S_{\rm max}\leq S_c$ are two threshold positive values. We shall prove the existence of an optimal control and we shall provide a necessary condition for a control to be optimal that, in particular, allows its identification via numerical simulations.

1.4.   Functional framework

We will use the classical Lebesgue spaces $L^{p}(\Omega)$ ($p\geq 1$), denoting their norms by $\|\cdot\|_{L^p}$, and the Sobolev spaces $H^{k}(\Omega)$ of functions in $L^{2}(\Omega)$ with distributional derivative of order less than or equal to $k$ in $L^{2}(\Omega)$. As customary, we set $H = L^{2}(\Omega)$ with inner product denoted by $(\cdot, \cdot)$ and corresponding norm $\Vert \cdot \Vert$. We also set $V = H^{1}(\Omega)$ equipped with the norm

and by $V'$ its dual space, the symbol ${\langle} \cdot, \cdot{\rangle}$ standing for the corresponding duality pairing.

Finally, we set

We will also make use of spaces of functions that depend on time with values in a Banach space. Hence, given a generic Banach space $B$ with norm $\|\cdot\|_B$ and an interval $I\subseteq [0, \infty)$, $L^p (I; B)$ is the set of measurable functions $f : I\to B$ such that $t\mapsto \|f(t)\|_{B}$ belongs to $L^p(I)$. Recall that $L^2(I; H)$ is isomorphic to $L^2(\Omega\times I)$. With the symbol $W^{1, p}(I; B)$ we will denote functions $f : I\to B$ such that both $f$ and its (weak) derivative ${\partial_t} f$ belong to $L^p (I; B)$. The family of continuous functions $f : I \to B$ is denoted by $\mathcal{C}(I, B)$.

Throughout the paper, by $C > 0$ we shall denote a constant that may change from line to line, depending on the problem parameters, the final time $T$, the norms of the initial data, and possibly on the norms of ${{\textsf{c}}}$ and ${{\textsf{s}}}$.

2.   Well-posedness: statement of the results

For a fixed $T > 0$, we set $Q_T = \Omega\times (0, T)$ and we introduce the phase space

Definition 2.1. Let ${{\textsf{c}}}\in L^\infty(Q_T)$, ${{\textsf{s}}}\in L^\infty(Q_T)$ be given and take $(u_0, \sigma_0)\in H\times H$. A (weak) solution on $[0, T]$ to the initial value problem (1.2) endowed with Neumann boundary conditions is a pair $(u, \sigma)$ with

satisfying

for almost every $t\in (0, T)$. Moreover, ${{\partial_{\bf n}}} {u} = {{\partial_{\bf n}}} \sigma = 0$ almost everywhere on $\Sigma_T$ and $({u}(0), \sigma(0)) = ({u}_0, \sigma_0)$ almost everywhere in $\Omega$.

Remark 2.2. By a classical result (see, e.g., ^[24]), the regularity $f\in L^2(0, T; V)$, $\partial_t f\in L^2(0, T; V')$ ensures that $f\in{{\mathcal C}}([0, T], H)$. Besides, $t\mapsto \|u(t)\|^2$ is absolutely continuous and $\frac{{{{\rm{d}}}}}{{{{\rm{d}}}} t}\|f\|^2 = {\langle} \partial_t f, f{\rangle}$.

Theorem 2.3. Let ${{\textsf{c}}}\in L^\infty(Q_T)$, ${{\textsf{s}}}\in L^\infty(Q_T)$ with ${{\textsf{s}}}\leq S_c$ be given, and $({u}_0, \sigma_0)\in H\times H$ be such that

Then, system $(1.2)$ has a unique weak solution $({u}, \sigma)\in X\times X$ such that

Besides, the following uniform estimate holds:

Furthermore, if $\sigma_0\in L^\infty(\Omega)$, then $\sigma\in L^\infty(Q_T)$ and

In particular, our theorem tells that the biologically relevant region

is invariant for the differential system (1.2), namely: if we consider any biologically meaningful initial datum $z_0 = (u_0, \sigma_0)\in{\mathcal S}$, then any weak solution of (1.2) departing from $z_0$ remains in ${\mathcal S}$ for every time.

The next Section 3 and Section 4 are devoted to the proof of Theorem 2.3 via a number of steps. The first consists in the introduction of an auxiliary problem where we suitably modify some of the nonlinearities involved in system (1.2).

3.   An auxiliary problem

In this and the next section, according to the assumptions of Theorem 2.3, we assume $({{\textsf{c}}}, {{\textsf{s}}})\in L^\infty(Q_T)\times L^\infty(Q_T)$, with ${{\textsf{s}}}\leq S_c$, fixed. We introduce the cut–off function

Note that $k$ is globally bounded and Lipschitz on $\mathbb{R}$. Then, defining

we consider the auxiliary system

Due to the special form of the nonlinearities, it is easy to show that any solution of system (3.1) originated from an initial data $(u_0, \sigma_0)\in {\mathcal S}$ belongs to ${\mathcal S}$ a.e. in $Q_T$. This is done in the sext section.

3.1.   Maximum principles

Let $(u(t), \sigma(t))\in H\times H$ be any weak solution on $[0, T]$ to the auxiliary Cauchy problem (3.1).

Lemma 3.1. If $0\leq u_0\leq 1$ a.e. in $\Omega$, then $0\leq u(t)\leq 1$ a.e. in $Q_T$.

Proof. Testing the first equation of (3.1) by $-u^-\in V$, where $u^- = \max (0, -u)$, we obtain

Indeed, the rhs is identically zero since $k(u) = 0$ whenever $u\leq 0$. As a result,

and the Gronwall lemma yields

This means that $u\geq 0$ a.e. in $Q_T$. Now we consider $w = u-1$. Note that $w$ solves the equation

hence, testing by $w^+$ we get

Again the rhs is identically zero since by construction $k(u) = 0$ whenever $u\geq 1$, namely, on the support of $w^+$. Reasoning as above, we reach the conclusion

Since $w(0) = u_0-1\leq 0$, then $w^+(0) = 0$: as a consequence $w^+(t) = 0$ in $\Omega\times [0, T]$, meaning that $u\leq 1$ a.e., as claimed.

Lemma 3.2. If $0\leq u_0\leq 1$ and $\sigma_0\geq 0$ a.e. in $\Omega$, then $\sigma (t)\geq 0$ a.e. in $Q_T$.

Proof. Testing the second equation of (3.1) by $-\sigma^-\in V$, and taking into account that $0\leq u\leq 1$ a.e., we get

Integrating over $[0, t]$ the final differential inequality $\frac{{{{\rm{d}}}}}{{{{\rm{d}}}} t}\|\sigma^-\|^2\leq 2|\gamma_{ch}|\|\sigma^-\|^2$, we obtain

for all times. Since by assumption $\sigma_0\geq 0$, it turns out that $\sigma^-(0) = 0$, yielding the thesis.

3.2.   Local in time existence

Let $(u_0, \sigma_0)\in H\times H$ be arbitrarily given. In this section we prove that the Cauchy problem (3.1) admits (at least) a local solution, which is defined in a maximal time interval $[0, \tau)$, for some $\tau > 0$.

3.2.1.   Galerkin approximation

Let $\{e_j\}_{j = 1}^\infty$ be a smooth orthonormal basis in $H$ which is also orthogonal in $V$. Then define $V_n = {\rm Span}\{e_1, \dots, e_n\}$ and denote by $\mathbb P_n$ the corresponding projection. Now, for any fixed $n\in {\mathbb{N}}$, we consider the following finite dimensional problem: Find $t_n > 0$ and functions $a_j, b_j\in {{\mathcal C}}^1([0, t_n))$ such that

satisfy, for almost every $t\in (0, t_n)$,

and

for every test function $v\in V_n$ and $w\in V_n$, along with the initial conditions

Indeed, choosing $v = w = e_j$ for any $j\in\{1, \ldots, n\}$, everything boils down to a system of $2n$ nonlinear ordinary differential equations with locally Lipschitz nonlinearities. Hence, by classical results in ODE's theory, the local existence (and uniqueness) of a solution $(u_n, \sigma_n)$ is guaranteed on a certain maximal interval $[0, t_n)$. Besides, the solution satisfies

We now wish to find estimates that are independent of $n$.

3.2.2.   Energy estimates

Along the proof, $C > 0$ will stand for a generic constant independent of $n$. Test the first equation by $u_n$ and the second one by $\sigma_n$ to find

We now estimate the rhs. Recalling that $|k(r)|\leq c(1+r^2)$ and that $\|\tilde{f}(\sigma_n, {{\textsf{c}}})\|_{L^\infty}\leq C$ by construction, we have

Besides, we easily obtain

and

It is then apparent that we end up with the inequality

There we can control the $L^3$-norm of $u_n$ via the Young inequality with exponents $\frac43, 4$ as follows

Thus we arrive at

having set

Therefore, by Gronwall's lemma,

so that there exists $\tau > 0$ independent on $n$ such that

Since there exists $C > 0$ such that

then, going back to (3.2) and integrating in time over $[0, \tau]$, we further learn that

and

Hence, by comparison,

3.2.3.   Passage to the limit

Due to uniform bounds above, there exists $(u, \sigma)$ such that

By the uniform control on ${\partial_t} u_n$ and ${\partial_t}\sigma_n$, we also learn that

so, in particular,

This allows to pass to the limit in the weak formulation to prove that $(u, \sigma)$ is a weak solution of (3.1) on $[0, \tau]$. All the convergences are straightforward but those involving the nonlinear terms. We start by proving that for any $w\in V$, for any $\varphi\in C^\infty_0(0, t)$ with $t\leq \tau$,

This is easily seen by noticing that $f_n = u_n^3-u^3\to 0$ a.e. in $Q_\tau$ and $\|f_n\|_{L^{4/3}(Q_\tau)}\leq C$, hence $f_n\to 0$ weakly in $L^{4/3}(Q_\tau)$ (Lebesgue convergence, weak form). Let us now prove that

To this end, observe that, since $\tilde{f}$ is globally Lipschitz and $k$ is a bounded function,

by the strong convergence $\sigma_n\to \sigma$ in $L^2(0, \tau; H)$. Analogously, exploting the fact that $\tilde{f}$ is bounded and $k$ globally Lipschitz

by the strong convergence $u_n\to u$ in $L^2(0, \tau; H)$.

In the second equation, we have to show that

Indeed, we rewrite the difference as

and, noticing that $\frac{|s|}{1+|s|}\leq 1$, we obtain

where

On account of the Lipschitz continuity of $\frac{s}{1+|s|}$, we find

Hence, the last term in the second equation can be handled as

3.3.   Global existence

In this section we show that any solution to the auxiliary initial value problem (3.1) originated from $(u_0, \sigma_0)\in {\mathcal S}$ is defined for all positive times.

Theorem 3.3. Let $T > 0$ be given and $z_0 = (u_0, \sigma_0)\in {\mathcal S}$. Then, any weak solution $(u, \sigma)$ to $(3.1)$ departing from $z_0$ is global in time on $[0, T]$.

Proof. Let us define

We know by the previous section that there exists a solution $(u, \sigma)\in {\mathcal S}$ defined on $[0, \tau]$, hence $\overline{t}\geq \tau > 0$ and

Besides, inequality (3.2) holds for $(u, \sigma)$ in place of $(u_n, \sigma_n)$ since all the involved constants are independent of $n$. Then Gronwall's lemma yields

for some $c > 0$ independent of $\overline{t}$. This implies that $\overline{t} = T$: indeed, the uniform bounds of the $H$-norms tell that $\lim_{t\to \overline{t}}u(t)$ and $\lim_{t\to \overline{t}}\sigma(t)$ exist in $H$ (at least for a subsequence). Now we can consider a solution to the Cauchy problem with initial datum $(u(\overline{t}), \sigma(\overline{t}))$, which is defined on an interval $[\overline{t}, \overline{t}+\delta])$, for some $\delta > 0$ (see also the extension theorem [14, Lemma 3.1, p. 13]). In this way we contradict the definition of $\overline{t}$.

4.   Well-posedness: proof of Theorem 2.3

Let $T > 0$ and take an initial datum $(u_0, \sigma_0)\in {\mathcal S}.$ In light of the above analysis, let $(u, \sigma)$ be any global weak solution to the auxiliary system (3.1) originated by $(u_0, \sigma_0)$. We actually proved in Section 3.1 that $(u, \sigma)\in {\mathcal S}$ a.e. on $Q_T$, implying that

It turns out that the pair $(u, \sigma)$ actually solves the original problem (1.2) on $[0, T]$. This proves the first part of Theorem 2.3, namely the global existence of solutions to the original model and the invariance of the set ${\mathcal S}$. Let us now prove uniform energy estimates and a continuous dependence result, where we highlight the role of the controls for further use.

Theorem 4.1. Let $(u, \sigma)$ be a solution to $(1.2)$ originated from $(u_0, \sigma_0)\in {\mathcal S}$. Then, the following uniform estimate holds

Proof. Along the line of Section 3.2.2, we multiply the first equation of (1.2) by $u$ and the second one by $\sigma$ to find

At this point we estimate the rhs, recalling that $0\leq u\leq 1$, but now paying attention to the role of ${{\textsf{c}}}$ and ${{\textsf{s}}}$. Since

and $m$ is bounded, we have

Besides, we easily get

and

Hence, calling

we obtain the differential inequality

for some $\omega > 0$. An application of the Gronwall lemma on $[0, T]$ yields

for every $t\in [0, T]$, saying that

for every $t$. Going back to the differential inequality and integrating in time over $[0, T]$, we obtain the desired control for $u$ and $\sigma$ in $L^2(0, T;V)$ and $\|u \|_{L^4(Q_\tau)}$. Finally, by comparison in the system we get an analogous estimate for $\|({\partial_t} u, {\partial_t} \sigma)\|_{[L^2(0, \tau; V')]^2}\leq C,$ completing the proof.

Theorem 4.2. Let $(u_i, \sigma_i)$ be two solutions to to $(1.2)$ corresponding to controls $({{\textsf{c}}}_i, {{\textsf{s}}}_i)\in L^\infty(Q_T)\times L^\infty(Q_T)$, with ${{\textsf{s}}}_i\leq S_c$, and initial data $z_i = (u_{0, i}, \sigma_{0, i})\in {\mathcal S}$, $i = 1, 2$. Then, the following continuous dependence estimate holds:

for all $t\in [0, T]$.

Proof. Observe that, by Section 3, we know that $0\leq u_i\leq 1$ and $\sigma_i\geq 0$ a.e. in $Q_T$, for $i = 1, 2$. Let us denote by $(u, \sigma) = (u_1-u_2, \sigma_1-\sigma_2)$ and $({{\textsf{c}}}, {{\textsf{s}}}) = ({{\textsf{c}}}_1-{{\textsf{c}}}_2, {{\textsf{s}}}_1-{{\textsf{s}}}_2)$. Then

Recalling that $\tilde f(\sigma, {{\textsf{c}}}) = [1-3(m(\sigma)-{{\textsf{c}}})],$ we rewrite

Then, multiplying the first equation by $u$, taking into account that $u_1+u_2\leq 2$, we get

Since $m$ is a globally Lipschitz function and $k(u)$ is bounded, we immediately get

Besides, exploiting the global Lipschitz continuity of $k$ and the boundedness of $\tilde f$,

We thus end up with

As a second step we consider the second equation in the differential system solved by $(u, \sigma)$. We observe that

hence a multiplication by $\sigma$ yields

Since $|g(\sigma_1)-g(\sigma_2)|\leq |\sigma|$ and $0\leq u_1\leq 1$, the first term on the rhs is easily estimated as

We proceed noticing that, since $0\leq g(\sigma) < 1$,

Finally,

Collecting everything, we end up with the differential inequality

for some $\omega > 0$. Let now $T > 0$ be fixed. An application of the Gronwall lemma on $[0, T]$ yields

where $u(0) = u_{0, 1}-u_{0, 2}$ and $\sigma(0) = \sigma_{0, 1}-\sigma_{0, 2}$, proving the claimed continuous dependence estimate.

As an immediate consequence of Theorem 4.1, we see that the (global) weak solution to (1.2) departing from any $(u_0, \sigma_0)\in\mathcal S$ is unique. Indeed, let us denote by $(u_i, \sigma_i)$, $i = 1, 2$ two solutions, corresponding to a fixed pair of controls $({{\textsf{c}}}, {{\textsf{s}}})\in L^\infty(Q_T)\times L^\infty(Q_T)$ with ${{\textsf{s}}}\leq S_c$, departing from the same initial datum $z_0 = (u_{0}, \sigma_{0})\in {\mathcal S}$. Then, setting ${{\textsf{c}}}_1 = {{\textsf{c}}}_2 = {{\textsf{c}}}$, ${{\textsf{s}}}_1 = {{\textsf{s}}}_2 = {{\textsf{s}}}$ and $z_1 = z_2 = z_0$ in the continuous dependence estimate, we get

saying that $(u, \sigma)\equiv (0, 0)$ hence $(u_1, \sigma_1)\equiv(u_2, \sigma_2)$ in $[0, T]$.

As a last step, we are left to prove that, if the initial datum $\sigma_0$ is bounded, then the solution remains bounded for all times. To this aim, we rewrite the equation for $\sigma$ as

noticing that the right-hand side belongs to $L^\infty(0, T;H)$. If we consider $\sigma_0\in L^\infty(\Omega)$, then, by a classical result in the theory of linear parabolic PDEs (see e.g. Theorem 7.1 in ^[16]), we immediately find the desired conclusion $\sigma\in L^\infty(Q_T),$ together with the uniform estimate (2.1). The proof of Theorem 2.3 is now completed.

5.   The control-to-state mapping

From now on, let

be fixed. In light of the existence result Theorem 2.3, we can define the control-to-state mapping as

where $(u, \sigma)$ is the unique weak solution to (1.2) corresponding to $({{\textsf{c}}}, {{\textsf{s}}})$ with initial datum $(u_0, \sigma_0)$ as above. Here we set

Besides, by Theorem 2.3 we also know that

Observe that the mapping $G$ is Lipschitz continuous (having endowed $L^\infty(Q_T)$ with the $L^2$- topology); indeed, owing to Theorem 4.2 we have

for all $({{\textsf{c}}}_i, {{\textsf{s}}}_i)\in \mathcal U$ and associated states $(u_i, \sigma_i) = G({{\textsf{c}}}_i, {{\textsf{s}}}_i)$.

Let us now show that $G$ possesses certain directional derivatives at any point in $\mathcal U$. To this aim, let $({{\textsf{c}}}^*, {{\textsf{s}}}^*)\in \mathcal U$ be fixed and denote by $({u}^*, \sigma^*) = G({{\textsf{c}}}^*, {{\textsf{s}}}^*)$ the corresponding state. Then, for $({{\textsf{c}}}, {{\textsf{s}}})\in \mathcal U$, we introduce the linearized system at $({u}^*, \sigma^*)$, defined as

where the coefficients are defined as follows

and $S_{ch} = S_c-S_h.$ Notice that the four coefficients, as well as the source terms $-{{\textsf{c}}} h'({u}^*)$ and $-{{\textsf{s}}}{u}^*$, are in $L^\infty(Q_T)$, due to the assumptions on the nonlinearities and the fact that $0\leq u^*\leq 1$ and $\sigma^*\geq 0$ a.e. in $Q_T$. By the theory of linear parabolic equations (see the subsequent Theorem 7.1), there exists a unique strong solution to (5.3) with

At this point, we consider any $({\bar {{\textsf{c}}}}, {\bar {{\textsf{s}}}})\in {{\mathcal U}}$ and notice that

for any $\mu \in (0, 1)$. Therefore, we can consider the corresponding state $({u^\mu}, {\sigma^\mu}) = G({{{\textsf{c}}}^\mu}, {{{\textsf{s}}}^\mu})$ satisfying all the results proven in the previous sections. Note that, letting $\mu\to 0$, by construction ${{{\textsf{c}}}^\mu}\to {{\textsf{c}}}^*$ and ${{{\textsf{s}}}^\mu}\to{{\textsf{s}}}^*$ in $L^2(Q_T)$ As a consequence, since $G$ is Lipschitz continuous by (5.2), we have

Lemma 5.1. In the limit $\mu \to 0^+$ we have

where $(Y, Z)$ is the solution to the linearized system $(5.3)$ with $({{\textsf{c}}}, {{\textsf{s}}}) = (\bar {{\textsf{c}}}-{{\textsf{c}}}^*, \bar {{\textsf{s}}}-{{\textsf{s}}}^*)$.

Proof. We set

Accordingly, we have to prove that $Y^\mu\to 0$ and $Z^\mu\to 0$ in ${\mathcal H}$. The first step consists in writing in a suitable form a differential system for $(Y^\mu, Z^\mu)$. After some computations, it is not difficult to check that the following holds:

having defined

and

Here, $x^\mu$, $x_\mu$ and $s^\mu$, $s_\mu$ are measurable functions arising from the application of an extension of Lagrange Theorem (see [4, Appendix]) as follows:

We recall that $x^\mu$ and $x_\mu$ attain intermediate values between the ones of $u^\mu$ and $u^*$, while $s^\mu$ and $s_\mu$ are in between $\sigma^\mu$ and $\sigma^*$.

As a second step, we test system (5.6) with the pair $(Y^\mu, Z^\mu)$, so obtaining the differential equality

We proceed by estimating the rhs. Since $A_i, B_j\in L^\infty(Q_T)$ in light of (5.1) and thanks to the regularity of the involved nonlinearities, the first three terms in the rhs are easily controlled by

Besides, the last two terms can be estimated exploiting the Lipschitz continuity of $h'$ as

Finally, exploiting the fact that $Y, Z\in L^\infty(0, T;V)$ by (5.4),

Integrating on $[0, t]$, observing that $\|Y^\mu(0)\|^2+\|Z^\mu(0)\|^2 = 0$, we get

having set

We claim that

Indeed, we have the convergences (5.5), implying in turn that

Now the conclusion follows by invoking the Lipschitz continuity of $F", h', h"$ and of $m', g, g'$. As a final step we apply the Gronwall lemma on $[0, T]$ that yields

Letting $\mu\to 0$ the proof is done.

6.   Optimal control problem

Our optimal control problem consists in finding the control functions ${{\textsf{c}}}^*$ and ${{\textsf{s}}}^*$ (if any) that provide the optimal cytotoxic and antiangiogenic effects to treat a certain glioma whose evolution is modeled by (1.2).

In order to state the problem, we first fix the desired targets for the tumor phase and for the oxygen in $Q_T$ and in $\Omega$ at the final time $T$, respectively given by

Then, for any $(u, \sigma)\in [{{\mathcal C}}([0, T], H)]^2$ and any $({{\textsf{c}}}, {{\textsf{s}}})\in [L^2(0, T;H)]^2$, we introduce the functional

where $k_i$ are given nonnegative constants, with at least one strictly positive.

Next, we define the set of all the admissible controls $({{\textsf{c}}}, {{\textsf{s}}})$. Given two positive thresholds $U_{\rm max} > 0$ and $0 < S_{\rm max}\leq S_c$, we define

and we set

Finally, given the initial state

we consider the control-to-state map defined in Section 5. Accordingly, for any pair $({{\textsf{c}}}, {{\textsf{s}}})\in {{\mathcal U}}_{ad}$ we set $(u, \sigma) = G({{\textsf{c}}}, {{\textsf{s}}})$ as the weak solution to (1.2) corresponding to $({{\textsf{c}}}, {{\textsf{s}}})$ with initial datum $(u_0, \sigma_0)$, and we define the (reduced) cost functional

Our control problem can be stated as follows: find, if possible, an optimal control $({{\textsf{c}}}^*, {{\textsf{s}}}^*)\in {{\mathcal U}}_{ad}$ such that

We start the analysis by proving an existence result.

6.1.   Existence of the optimal control

Theorem 6.1. Under assumption $(6.1)$, for any fixed $(u_0, \sigma_0)\in {\mathcal S}$ with $\sigma_0\in L^\infty(\Omega)$ there exists at least a solution $({{\textsf{c}}}^*, {{\textsf{s}}}^*)\in {{\mathcal U}}_{ad}$ to $(6.2)$ with corresponding optimal state $(u^*, \sigma^*)$.

Proof. Indeed, since $J\geq 0$, it is immediate to see that $\inf_{({{\textsf{c}}}, {{\textsf{s}}})\in {{\mathcal U}}_{ad}} J({{\textsf{c}}}, {{\textsf{s}}}) = \delta\geq 0$. We can consider then a minimizing sequence $({{\textsf{c}}}_n, {{\textsf{s}}}_n)\in {{\mathcal U}}_{ad}$ such that

and, according to Theorem 2.3, the corresponding state $(u_n, \sigma_n)$: this, in particular is uniformly bounded in $X\times X$ with $0\leq u_n\leq 1$ a.e. in $Q_T$ and $\sigma_n\geq 0$ a.e. in $Q_T$ satisfies (2.1). By the boundedness of ${{\mathcal U}}_{ad}$ and Theorem 4.1, we can select subsequences (that we still denote as) $({{\textsf{c}}}_n, {{\textsf{s}}}_n)$ and $(u_n, \sigma_n)$ such that

It is worth noticing that so far no relation connects $({{\textsf{c}}}^*, {{\textsf{s}}}^*)$ and $(u^*, \sigma^*)$. Our aim will be to prove that $(u^*, \sigma^*) = G({{\textsf{c}}}^*, {{\textsf{s}}}^*)$, namely, that $(u^*, \sigma^*)$ is the state corresponding to the control, and that $J({{\textsf{c}}}^*, {{\textsf{s}}}^*) = \delta$. First of all, by compactness,

and, owing to the Ascoli-Arzelá Theorem,

Therefore, it follows that $(u^*(0), \sigma^*(0)) = (u_0, \sigma_0)$. Besides, $(u^*, \sigma^*)\in {\mathcal S}$ and $\sigma^*$ satisfies (2.1). Furthermore, due to the boundedness and Lipschitz continuity of all the involved nonlinear functions, these convergences allow to pass to the limit in the problem solved by $(u_n, \sigma_n)$, proving that $(u^*, \sigma^*)$ solves the initial boundary value problem correspondig to $({{\textsf{c}}}^*, {{\textsf{s}}}^*)$, that is, $(u^*, \sigma^*) = G({{\textsf{c}}}^*, {{\textsf{s}}}^*)$.

To accomplish our second task, we decompose the functional $J$ in three parts, namely,

where

Now, convergence (6.3) immediately gives

By the uniform boundedness of $\|u_n(T) \|$ following from Theorem 4.1, we infer that, up to a subsequence,

so that

The last functional $J_3$ is weakly lower semicontinuous thus

Collecting all our computations, we conclude

showing that indeed $J$ realizes its minimum value at $({{\textsf{c}}}^*, {{\textsf{s}}}^*)$.

Once the existence of an optimal control is established, the next goal is devising a necessary condition for a control to be optimal that, in particular, allows its identification by numerical simulations.

6.2.   Optimality conditions

Let $({{\textsf{c}}}^*, {{\textsf{s}}}^*)\in {{\mathcal U}}_{ad}$ be an optimal control and denote by $({u}^*, \sigma^*)$ the corresponding optimal state. Then, for any $({\bar {{\textsf{c}}}}, {\bar {{\textsf{s}}}})\in {{\mathcal U}}_{ad}$, we notice that

for any $\mu \in (0, 1)$. Therefore, we can consider the corresponding state $({u^\mu}, {\sigma^\mu})$ and observe that

Now, owing to Lemma 5.1, we can pass to the limit as $\mu \to 0^+$ in (6.5), saying that the derivative of $J$ at $({{\textsf{c}}}^*, {{\textsf{s}}}^*)$ in the direction of $({\bar {{\textsf{c}}}}-{{\textsf{c}}}^*, {\bar {{\textsf{s}}}}-{{\textsf{s}}}^*)$ is nonnegative. Invoking Lemma 5.1 once again, we easily obtain

where $(Y, Z)$ solves the linearized problem (5.3). The above inequality is the so-called first order optimality condition, although from its expression it is really difficult to identify the optimal control even by numerical simulations.

Therefore, as it is done in the classical control theory (see the subsequent Theorem 7.1), we introduce the so-called adjoint problem, here defined as

where $k_i$, $i = 1, \dots, 5$, and ${u}_Q, \sigma_Q, {u}_\Omega, \sigma_\Omega$ are exactly the constants and the target functions appearing in the cost functional.

By Theorem 7.1, we learn that there exists a unique weak solution $(w, z)\in X\times X$ to (6.7) and that the analogous of (7.1) holds true, namely,

As a consequence, inequality (6.6) turns into the much simpler form

that we write as

Notice that this is equivalent to

The geometric meaning of these inequalities is clear: indeed, leaning on the elementary theory of projections in Hilbert spaces (see e.g. Remark 4.6 in ^[4]), we have obtained the first order optimality conditions

Theorem 6.2 (First order optimality conditions). Let $({{{c}}}^*, {{{s}}}^*)\in {{\mathcal U}}_{ad}$ be an optimal control, with corresponding state $({u}^*, \sigma^*)$. Then

where $(w, z)\in X\times X$ is the solution to the adjoint system $(6.7)$.

Let us conclude our analysis by expressing the optimal control in the easiest possible form, recalling that the projection of $v\in L^2(Q_T)$ into

is

As a result, the optimal control is characterized by the following two formulas

and

7.   Appendix. A useful result for distributed controls

In classical control theory (see e.g. ^[25]), the more feasible expression of the first order optimality condition relies on the solutions to a suitable linear problem and its adjoint. We briefly describe the main tools in a suitable form to treat the model under study in this paper. We consider the two linear systems

where the unknowns are the vectors $y = (Y, Z)^T$ and $p = (w, z)^T$ so that, in particular, $\Delta y = (\Delta Y, \Delta Z)^T$. Besides, $c_0$ is a $2\times 2$ matrix whose transpose is $c^T_0$, while all the other given quantities $b, v, a_Q, a_\Omega$ are vector functions. The second system $(L^*)$ is called the adjoint to system $(L)$. The link between the two systems, that turns out to be quite useful in order to identify the optimal control, is expressed by (7.1) in the next result.

Theorem 7.1. Provided that the entries of the matrix $c_0$ and of the vectors $b$ and $v$ belong to $L^\infty(Q_T)$, if $a_Q\in [L^2(Q_T)]^2$ and $a_\Omega\in H\times H$, then there exists a unique strong solution $y\in [{{\mathcal C}}([0, T], V)\cap L^2(0, T;W)]^2$ to problem $(L)$ such that

Moreover, there exists a unique weak solution $p\in X\times X$ to the adjoint problem $(L^*)$ such that

Finally, the following equality holds true

Proof. First of all, we see that problem $(L)$ is well posed: indeed, under our assumptions, the coefficients belong to $L^\infty(Q_T)$, the source terms to $L^2(Q_T)$ and the null initial data are in particular in $V\times V$, hence by classical results on parabolic systems (see, e.g., [15, Theorem 1.1]), there exists a unique strong solution $y\in [{{\mathcal C}}([0, T], V)\cap L^2(0, T;W)]^2$ satisfying (7.1). Reversing time by the change of variable $t\mapsto T-t$, problem $(L^*)$ turns into a forward system with $L^\infty(Q_T)$ coefficients, $L^2(Q_T)$ source terms and initial data in $H\times H$. Thus, the aforementioned theorem on linear parabolic system applies, yielding the existence of a unique weak solution $p\in X\times X$ to $(L^*)$ satisfying (7.1). Since $y$ is a strong solution to problem $(L)$ then it is also a weak solution and, by definition,

for every $\phi\in X\times X$. Choosing $\phi = p$ as weak solution to problem $(L^*)$ and integrating by parts, we obtain

Besides,

We thus find

Collecting our computations, we obtain (7.1).

Conflict of interest

All authors declare no conflicts of interest in this paper.