Optimal control of a chemotaxis equation arising in angiogenesis

1.   Introduction

Taxis is understood as the motion of an organism towards or away from an external stimulus. In particular when the stimulus is a chemical, it is called chemotaxis. It seems natural to address the problem of driving the motion of the organism by modifying or applying a chemical gradient. In this paper we have in mind a process in which chemotaxis takes place, the tumor angiogenesis. However, we could extend the same idea for other biological processes where the chemotaxis is involved.

Tumor angiogenesis starts when, as a response to nutrient deprivation, cancer cells secrete a chemical factor known as Tumor Angiogenic Factors (TAF)s. (TAF)s diffuse in the extracellular matrix and activate endothelial cells which migrate, via chemotaxis, towards the source of TAF i.e., the tumor. When endothelial cells reach the tumor more nutrients are supplied to the tumor which grows further. See for instance [16] for additional details.

It is clear that angiogenesis is a crucial step in the development of the tumor. Therefore, if we act in the (TAF)s molecules or if we modify the response of the endothelial cells to the (TAF)s molecules we may either reduce angiogenesis or avoid it.

In this paper we consider two variables $u$, the density of endothelial cells and $z$ the concentration of anti-angiogenic drug that modifies the sensitivity of a chemical gradient $v$ (TAF) and the growth of the endothelial cells. We also assume that $u$, $z$ and $v$ are defined in a bounded domain $\Omega \subset {{{\rm{I}}}\hskip -0.85mm{{\rm{R}}}}^3$. Here we are interested in the minimum of the functional defined by

where $a$ and $b$ are positive parameters and $(u, z)$ is a solution of the nonlinear differential equation

and $z \in \mathcal{U}_{ad}$ where

The functional has two terms; the first term refers to the amount of endothelial cells at the final stage and the second one can be seen as the cost of the drug over the time.

Minimizing a functional (1.1), subject to a differential problem (1.2) and a convex constraint, $z\in \mathcal{U}_{ad}$, is a problem of the optimal control theory. It is classical in this theory to call $u$ state variable and to call $z$ control variable, because it controls the state through the equation. The nonlinearity of (1.2) leads that the minimal problem is out of the convex framework. We will apply a generalization of the Lagrange multipliers theorem to obtain the optimality system which provides the necessary condition for the solution. This method, called Dubovitskii-Milyutin formalism (see [7]) has been mostly applied to ordinary differential equations (see [13]). In [6] it is applied to a linear partial differential problem with the feature that it is non well-posed. The optimal problem we study in this paper have the difficulty of the nonlinearity of the partial differential equation and the control constraint $z\in \mathcal{U}_{ad}$ which has empty interior in $L^2(Q)$.

The formalism gives an optimality system constituted by two partial differential coupled equations, one is the state equation given in (1.2) and the other one is a linear equation for the adjoint variable, and a condition for the optimal control given by a projection operator. We will prove the uniqueness of solution of the optimality system for $T$ small enough and this turns out the uniqueness of solution of the optimal problem.

The chemotaxis is described in (1.2) by the term ${\rm div} (V(u) \alpha (z) \nabla v)$. A similar description of the chemotaxis was introduced in [12] where there is an additional equation for $v$. Since, the minimal system proposed in [12] could generate singularities in finite time (see [11]) then to avoid the singularities it is proposed a model with a bounded drift term in [9]. In the case of angiogenesis, one of the first continuous models related to angiogenesis is introduced in [1]. A similar model to the one in the paper, without the variable $z$ but with an additional equation for chemoattractant is given in [4]. The case of a therapy $z$ satisfying an additional parabolic equation is considered in [17] and [5].

The structure of the paper is as follows. In Section 2 we prove the existence of a unique weak solution of (1.2), global in time and positive, fixed $z$. The optimal control problem, the existence of the optimal control, its characterization and the uniqueness when $T$ is small enough are studied in Section 3. In Section 4 we show some numerical simulations to illustrate the theoretical results.

2.   Global existence and uniqueness for the equation

Let ${\Omega}\subset {{{\rm{I}}}\hskip -0.85mm{{\rm{R}}}}^N$ be a regular bounded domain and $T > 0$ a fixed number. We will denote $Q = {\Omega}\times (0, T)$ and ${\Gamma} = {\partial}{\Omega}\times (0, T)$. In this paper, we consider the problem

Here, $v\in W^{2, \infty}({\Omega})$ is a known function, $z\in L^2(Q)$, and $u_0\in\, \, L^2({\Omega})$. The function $V:\, [0, +\infty)\rightarrow [0, +\infty)$ verify $V(0) = 0$ and $V\in C^2([0, \infty))\cap L^\infty([0, \infty))$. If it is needed we will extend this function by zero for negative values. The functions ${\alpha}:\, [0, +\infty)\rightarrow [0, +\infty)$ is in $W^{1, \infty} ([0, +\infty))$ and ${\beta}:\, [0, +\infty) \times [0, +\infty) \rightarrow [0, +\infty)$ is in $C^2([0, \infty))\cap L^\infty([0, \infty))$. We would like to point out here that the main difficulty of the problem is that $z$ is not defined on ${\partial}{\Omega}$.

Let $T > 0$ and $X$ is a Banach space, we define the space $L^p(0, T; X)$ of equivalence classes of measurable functions $u:I\rightarrow X$ such that $t\in (0, T)\rightarrow \Vert u\Vert_X$ belongs to $L^p(I)$ which is a Banach space for the norm

For instance, we will use

The definition of a weak solution for the problem is

Definition 2.1. We say that $u$ is a weak solution of (1.2) if the following conditions are verified

$1)$ $u\in W(0, T)$ where

$2)$ $\forall w\in C^1([0, T]\times\bar{\Omega}):\, \, w(T, x) = 0, \, \forall x\in\bar {\Omega}$, it holds

Theorem 2.2. If $u_0 \in L^\infty (\Omega)$, $u_0\ge 0$, then there is a unique positive global weak solution to the problem (1.2) that also satisfies

Proof. We denote by

and we define the mapping $S:B_R\rightarrow L^2(0, T; L^2({\Omega}))$ such that for each $\bar u\in B_R$, $S(\bar u) = u$ is the weak solution of the linear problem

where

We will prove the existence of a fixed point for $S$ which is a weak solution of a truncated nonlinear problem. After that we will justify that the weak solution to the truncated problem is a solution of (1.2). At the end we will show the uniqueness of weak solution of (1.2).

$\underline {{\rm{Step 1}}}$ First of all we show that the operator $S$ is well defined. The existence and the uniqueness of the weak solution for the linear problem (2.3) follow from [14]. More precisely, by [14], Theorem 6.39] the problem (see [14], p. 136])

(where $f^i = -V(\bar u){\alpha}(z)\partial_iv$ and $c^0(x, t) = {\beta} (z, v)-T_k(\bar u)$ in (2.3) and $\nu = (\nu_1, .., \nu_N)$ is the outward normal to ${\Omega}$), has a unique weak solution i.e. $u\in L^2(0, T; H^1({\Omega}))\cap L^\infty(0, T; L^2({\Omega}))$ such that

for each $w\in C^1([0, T]\times\bar{\Omega})$ such that $w(T, x) = 0, \, \, \forall x\in {\Omega}$. Note that, if $u_t\in L^2(0, T; (H^1({\Omega}))')$ then the previous definition it is exactly the one given in Definition (2.1). Therefore we should show that $u_t\in L^2(0, T; (H^1({\Omega})')$.

where $\langle \cdot, \cdot \rangle$ stands for the duality pairing between $(H^1 (\Omega))'$ and $H^1 (\Omega)$. Therefore

and

i.e., $u_t\in L^2(0, T; (H^1({\Omega}))')$. Hence, $S$ is well defined.

In what follows we will apply the Schauder fixed point theorem to get the existence of a fixed point for $S$ in $B_R$.

$\underline {{\rm{Step 2}}}$ We claim that there exists $R > 0$ and $T > 0$ such that $S(B_R)\subset B_R$. In fact, multiplying the equation of (2.3) by $u$ and integrating in ${\Omega}$, it results

Thus for $y(t) = { }{\int_{\Omega} u^2}$ we obtain the problem

for some positive constants $a, \, \, b$. Consequently,

We can choose any $R > y(0) = \Vert u_0\Vert_2$ and determine $T > 0$ such that

and

if we take $R > 1$ and $T < 1$. So, for $R\geq \Vert u_0\Vert_2+1$ there exists $T < 1$ such that $S(B_R)\subset B_R$.

On the other hand, taking into account that

we can choose $R$, determine $T$ and integrate the inequality on the interval $(0, T)$ to obtain

Consequently, $\Vert u\Vert_{L^2(H^1)}$ is bounded.

$\underline {{\rm{Step}}\; 3}$ We claim that $S$ is a continuous mapping. We will prove that

with $\Phi(k, s)\to 0$ when $s\to 0$ for each fixed $k$, and $T$ as in Step 2. If we denote $u_1 = S(\bar u_1)$ and $u_2 = S(\bar u_2)$, it holds

and taking $w = u_1-u_2 = S(\bar u_1)-S(\bar u_2)$, we have

On multiplying the previous inequality by $w$ and integrating on ${\Omega}$ we obtain

Next, we provide a bound to the terms in the right hand side of (2.9)

The Sobolev embedding $H^1({\Omega})\hookrightarrow L^6({\Omega})$ together with the estimate $u_2\in L^2(0, T; H^1({\Omega}))$ implies that $u_2 \in L^2 (0, T; L^6 (\Omega))$. On the other hand if we apply the Hölder inequality with exponents $3$ and $3/2$ to the last term of the above inequality we obtain

Using the interpolation inequality,

and, thus

Then (2.9) becomes

We denote by $y(t) = \int_{\Omega} w(x, t)^2\, dx$. Since $y(0) = 0$ it follows from (2.10) that

By the Hölder inequality it holds

Since

it results from (2.11),

Integrating on $(0, T)$, we obtain

which proves the desired continuity.

$\underline {{\rm{Step}}\; 4}$ We claim that $S$ is a compact mapping in $L^2 (Q)$. We know that for each $\overline u \in B_R$ $S(\overline u) = u$ is bounded in $L^2 (0, T; H^1 (\Omega))$ and $S(\overline u)_t = u_t$ is bounded in $L^2 (0, T; (H^1(\Omega))')$ (see (2.7) and (2.5)) then by the Lions-Aubin Lemma (see for instance [18]) $S(B_R)$ is embedded compactly in $L^2 (0, T; L^2 (\Omega))$.

$\underline {{\rm{Step}}\; 5}$ By the Schauder Theorem we have a solution to the problem

We will prove that any solution of (2.13) is nonnegative. Let $u^- = \min (u, 0)$. Taking $u^-$ as a test function in the above equation we obtain

Integrating on the space variable infer

Therefore, $u^- (t) = 0$ a.e. in $\Omega$.

$\underline {{\rm{Step}}\; 6}$ We prove that a solution to (2.3) satisfies $u\in L^\infty(0, T; L^\infty({\Omega}))$ independently of the $k$, as a consequence a solution of (2.13) is a solution to (1.2) for $k$ sufficiently large. We multiply (2.3) by $pT_k(u)^{p-1}$ for $p\geq 2$, a cut off function that approximate to $pu^{p-1}$ (see for instance [3], p. 1196]) to get

Since $V(0) = 0$ and $V$ is a Lipschitz function then,

Adding on both sides of the above inequality $\int_\Omega T_k(u)^p$ we deduce that

At this point we will provide a bound for $\int_\Omega u^p$. We claim that for any $\epsilon > 0$

Next lines are devoted to the proof of the above inequality. Let $\overline z = \frac{1}{|\Omega|} \int_\Omega z$. We know that

On the other hand, by the Hölder inequality and the Poincare-Wintinger inequality we have

Therefore if we take $z = T_k(u)^{p/2}$ the claim easily follows. By (2.15) we get from (2.14) that

For any $0\leq t\leq T < T_{max}$, where $T_{max}$ stand for the maximal existence time, the above inequality asserts

Let us define

From (2.16), by a recursive procedure, we have

In particular, for $p = 2^j$ we deduce

where $s_j = \sum_{i = 0}^j 2^{-i}$. Taking into account that

Therefore we can take $k\rightarrow +\infty$ to get

Since

we can take $j \rightarrow \infty$ to conclude that

$\underline {{\rm{Step}}\; 7}$ Uniqueness of solution for (1.2). Let be two solutions of (1.2), $u_1, \, \, u_2$, and we take $w = u_1 -u_2$. We can argue as in Step 3 to obtain

for some $a' > 0$. By the boundedness of $u_1 +u_2$ we easily get the result by the Gronwall Lemma.

3.   The optimal control problem

In this section we are interested in the minimum of the functional

defined by

$a$ and $b$ are positive parameters and $(u, z)$ has to be a solution of the nonlinear differential equation

and besides, we claim $z$ that

The function $u$ represents the endothelial cell (EC) density, the initial data $u_0\in L^{\infty}(\Omega)$, $u_0\ge 0$ and not identically zero, means that the angiogenesis process has already begun, the TAF concentration is given by $v$ and it is known in the equation and the concentration of the anti-angiogenic drug is represented by $z$ which is the control function. Minimizing $J$ means to find the appropriated therapy $z$ such that the total amount of EC at the end of the treatment, this is $T$, is the smallest as possible applying the least amount of drug possible. The term $\int_{\Omega}u(T)^2$ represents the total amount of EC in $\Omega$ at the final time $T$ and the term $\int_Qz^2$ plays the total amount of drug in $\Omega$ for the time interval $(0, T)$.

3.1.   The existence of solution of the optimal problem

The proof of the existence of solution of the optimal problem (3.1)–(3.3) is standard by a minimizing sequence. We will do the proof for the reader's convenience.

Theorem 3.1. The optimal control problem (3.1)–(3.3) has a solution $(\hat{u}, \hat{z})$ in $W(0, T)\times L^2(Q)$.

Proof. Let be $\{u_n, z_n\}$ a minimizing sequence such that $\{J(u_n, z_n)\}$ is decreasing to $\hat{\beta}$, which is the infimum of $J(u, z)$ subject to (3.2) and (3.3). Then, $\{u_n(T)\}$ is bounded in $L^2(\Omega)$ and $\{z_n\}$ is bounded in $L^2(Q)$. So, there exist $\hat {z}\in L^2(Q)$ and a subsequence of $\{z_n\}$ which converges to $\hat {z}$ in $L^2(Q)$-weakly. Multiplying the differential equation by $u_n$ and integrating by parts, we obtain

Using that $\alpha, V\in L^{\infty}({{{\rm{I}}}\hskip -0.85mm{{\rm{R}}}})$, $\beta\in L^{\infty}({{{\rm{I}}}\hskip -0.85mm{{\rm{R}}}}^2)$ and $u_n\ge 0$, we have

and so,

Applying Gronwall's Lemma we deduce that $\{u_n\}$ is bounded in $L^{\infty}(0, T; L^2(\Omega))$. Using this fact and returning to (3.4) we obtain that $\{u_n\}$ is bounded in $L^2(0, T; H^1(\Omega))$. This boundness gives that $\{(u_n)_t\}$ is bounded in $L^2(0, T; (H^1(\Omega))')$. Effectively, taking any $w\in H^1(\Omega)$ and applying the duality product $\langle H^1(\Omega)', H^1(\Omega)\rangle$ in the equation of (3.2) we have

Since $\{u_n\}$ is bounded in $L^{\infty}(0, T; L^2(\Omega))$, $\{u_n\}$ is bounded in $L^2(0, T; H^1(\Omega))$, applying the Hölder inequality, we get

Then, we obtain

So,

which proves that $\{(u_n)_t\}$ is bounded in $L^2(0, T; (H^1(\Omega))')$.

And this gives that $\{u_n\}$ is bounded in $W(0, T)$. By the Aubin-Lions Lemma (see for instance [15]) we can assure that there exist a subsequence of $\{u_n\}$ and a function $\hat{u}$ in $W(0, T)$ such that (still denoting by $u_n$ the subsequence) $\{u_n\}$ converges to $\hat{u}$ in $L^2(Q)$ strongly.

The variational formulation of the problem (3.2) for the functions $(u_n, z_n)$ is

The functions $\alpha$ and $\beta$ are convex and continuous (for $\beta$ it is only necessary to be convex and continuous in the first variable), $V$ is a continuous function in ${{{\rm{I}}}\hskip -0.85mm{{\rm{R}}}}$. Then, using the convergence of the sequences, $\{u_n\}$ and $\{z_n\}$, it is possible to pass to the limit in (3.5) and so, we obtain that $(\hat{u}, \hat{z})$ solves the problem (3.2). Besides, as $\{z_n\}$ weakly converges to $\hat{z}$ in $L^2(Q)$ and $z_n\ge 0$, we have that $\hat {z}\ge 0$. This proves that $(\hat{u}, \hat{z})$ belongs to the admissible set of the optimal control problem. The functional $J$ is continuous and convex, so

which gives that $J(\hat{u}, \hat{z}) = \hat{\beta}$ and then, $(\hat{u}, \hat{z})$ is a solution of the problem (3.1)–(3.3).

Remark 1. The solution of the optimal problem might not be unique.

3.2.   The optimality system

In order to obtain the first order necessary conditions of the optimal control problem, we enunciate the Duvobitskii-Milyutin theorem (see [7]), which is a generalization of the Lagrange's Multipliers theorem.

Theorem 3.2. Let $X$ be a normed space. Assume that the functional $J:X\to {{{\rm{I}}}\hskip -0.85mm{{\rm{R}}}}$ has a local minimum with constraints $Z = \bigcap\limits_{i = 1}^{n+1}Z_i \subset X$ at a point $x_0\in Z$. Assume that $J$ is regularly decreasing at $x_0$, with decreasing (and convex) cone $DC_0$; the inequality constraints $Z_i$, $1\le i\le n$, are regular at $x_0$, with feasible (and convex) cones $FC_i$, $1\le i\le n$; the equality constraint $Z_{n+1}$ is also regular at $x_0$, with tangent (and convex) cone $TC_{n+1}$. Then, there exist continuous linear functionals $f_0\in DC_0^{\ast}$, $f_i\in FC_i^{\ast}$ for $1\le i\le n$ and $f_{n+1}\in TC_{n+1}^{\ast}$ (we denote by $^{\ast}$ the corresponding dual cone), not all identically zero, such that they satisfy the Euler-Lagrange equation:

The constraint (3.3) must be considered as an equality constraint, together with (3.2), because the convex $U = \{z\in L^2(Q):\ z\ge 0\}$ is a set with empty interior in $L^2(Q)$. So, applying the Duvobitskii-Milyutin theorem to (3.1)–(3.3) there exist two functionals, $f_0\in (DC)^{\ast}$, $f\in (TC)^{\ast}$ such that

and they are not simultaneously identically zero.

The identification of the cones. Following [2], we are going to define the cones and their dual ones in a point. The functional $f_0\in (DC)^{\ast}$ where $(DC)^{\ast}$ is the dual decreasing cone in $(\hat{u}, \hat{z})$, is given by

The computation of the tangent cone and its dual one is technically more complicated because there are two constraints, (3.2) and (3.3), considered like equalities generating a tangent cone. Each of these constraints generates, by itself, a tangent cone. We denote $Z_1$ the constraint set given by (3.2) and $TC_1$ the tangent cone associated to $Z_1$ in $(\hat{u}, \hat{z})$. By Lyusternik's theorem (see [2]), this cone is the set

And $(TC_1)^{\ast}$ is

If we consider the convex constraint (3.3) and we call $Z_2$ to this set, the dual tangent cone to $Z_2$ in the point $(\hat{u}, \hat{z})$, denoted by $(TC_2)^{\ast}$, is

The question is if we could assure that the dual cone to the cone associated to the set $Z_1\bigcap Z_2$ in $(\hat{u}, \hat{z})$, denoted by $(TC)^{\ast}$, would be equal to $(TC_1)^{\ast}+(TC_2)^{\ast}$. We can always assure that $(TC_1)^{\ast}+(TC_2)^{\ast}$ is included in $(TC)^{\ast}$, but the converse inclusion will be true if we prove that $(TC_1)^{\ast}$ and $(TC_2)^{\ast}$ are a system of the same sense cones (see [19], Definition 2.2.1]). The definition of a system of the same sense cones (SSS) is the following (see [19]):

Definition 3.3. Let $\{C_i\}_{i = 1}^n$ a system of cones in a normed space $X$. It is a (SSS) if for every $M > 0$ there exist $M_1, ...M_n$ such that, for each $x = \sum_{i = 1}^n x_i$, $x_i\in C_i$, the inequality $\Vert x\Vert \le M$ implies the inequalities $\Vert x_i\Vert \le M_i, \ \ i = 1, ..., n$.

This definition is not handle to determine if a system of cones is or not a (SSS). We enunciate a theorem that characterizes a (SSS) constituted by two cones when one of them is given by a linear continuous operator (see [19]).

Theorem 3.4. Let $E = X\times Y$ be the product space of two normed spaces and $C_1, C_2$ the following cones:

where $A$ is a linear continuous operator from $Y$ to $X$, and

and $\tilde{C_2}$ is a cone in $Y$. If we denote by $(C_i)^{\ast}$ the dual cone to $C_i$, then

and

Besides, the system $\{(C_1)^{\ast}, (C_2)^{\ast}\}$ is a (SSS).

In our case, $C_1 = TC_1$ and $C_2 = TC_2$. The tangent cone associated to (3.2), $TC_1$, is given by a linear continuous operator, the operator of the differential equation (3.6). By this theorem, we have that $(TC)^{\ast} = (TC_1)^{\ast}+(TC_2)^{\ast}$.

In the following theorem we obtain the optimality system of the optimal control problem (3.1)–(3.3).

Theorem 3.5. If $(\hat{u}, \hat{z})$ is a solution of the optimal control problem (3.1)–(3.3), then there exists $p\in W(0, T)$ such that $(\hat{u}, \hat{z}, p)$ satisfies

Proof. The Euler-Lagrange equation

The Euler-Lagrange equation, (E-L), is

with $f_0\in (DC)^{\ast}$, $f_1\in (TC_1)^{\ast}$ and $f_2\in (TC_2)^{\ast}$. We know that

We will see that $\lambda$ cannot be zero.

Suppose that $\lambda = 0$. Then, the (E-L) equation is

Applying this equation to any $(u, z)\in TC_1$ we have $\langle f_2, (u, z)\rangle = 0$ and so, $\langle \tilde{f}_2, z \rangle = 0$ for every $z$ such that there exists $u$ verifying $(u, z)\in TC_1$. Since for every $z\in L^2(Q)$ there exists this $u$, we have that $\tilde{f}_2\equiv 0$ and then, $f_2\equiv 0$ and $f_1\equiv 0$, which is impossible. So, $\lambda \ne 0$, we can rescale $\lambda = 1$ and the (E-L) equation is

We have that

this is

In the following, we describe a standard procedure in order to rewrite the first term in (3.8) using the optimal control $\hat{z}$. This is done defining an adjoint problem. We multiply the differential equation of (3.6) by a function $p$, which will be the solution of the adjoint problem, and we integrate by parts. Taking account, we obtain

We define the adjoint problem

Then

and replacing this equality in (3.8) we obtain

for every $z\in L^2(Q)$ because, as we have already said, for every $z\in L^2(Q)$ there is a function $u$ verifying $(u, z)\in TC_1$. So, we have obtained the functional $\tilde{f}_2$:

Since $(\tilde{f}_2, z-\hat{z})\ge 0$ for every $z\ge 0$, we have that

and this is equivalent to

where ${\mathcal{U}}$ is the convex set $\{ z\in L^2(Q):\ \ z\ge 0\}$ and $P$ is the projection operator of $L^2(Q)$ on ${\mathcal{U}}$.

In this case, because of the particular convex set we have, the operator $P$ is well-known. It is

hence,

We have just obtained the optimality system, given by the theorem.

Next, we are going to prove the uniqueness of solution of (3.7) for $T$ small enough. To this aim it will be useful to obtain an uniform estimate in time for $\nabla p$ on the interval $[0, T]$. In order to get this estimate we will suppose additional restrictions on $\alpha$ and $v$. In particular, we require that $\alpha = \alpha_0$ is a constant function and $\partial_n v = 0$. In this case the optimality system is given by the following equations:

Lemma 3.6. We have that $\|u(T)\|_{W^{1, 3} (\Omega)} \leq C$.

Proof. Using the variations of constants formula in (3.11)–(3.13) we get

where $A = -\Delta + I$ with domain $\{ \psi \in W^{2, 3} (\Omega) : \partial_n u = 0 \}$ and

Let $X^\beta$, $\beta \in (0, 1)$ the fractional powers of $X$. By [8], Theorem 1.6.1] we have that

for any $\beta \in (1/2, 1)$. Taking the norm $W^{1, 3}$ in (3.18) and applying the above embedding in the right hand side we infer

By [8], Theorem 1.4.3] we deduce

Let us notice that by the regularity of $\beta$, $V$ and $v$ and the embedding $W^{1, 3} (\Omega)\hookrightarrow L^6(\Omega)$ we have

As a consequence, we can conclude by the singular Gronwall Lemma (see [8], Section 1.2.1]).

In what follows we will show the regularity for the backward linear parabolic equation (3.14)–(3.16). The backward equation can be written in a forward way by the change of variable $\overline p (t) = p(T-t)$. Therefore, the problem is

Lemma 3.7. We have that $\max\limits_{t \in [0, T]} \| p (t) \|_{W^{1, 3} (\Omega)} \leq C\, .$

Proof. As in the previous Lemma we apply the variation of constants formula for $\overline p$ to get

where

We take $W^{1, 3} (\Omega)$ in the variations of constants formula for $p$ and we apply [8], Theorem 1.4.3, Theorem 1.6.1] and ([10], p. 59]) to infer

with $\beta \in (1/2, 1)$. Since $\|f(u, v, \overline p)\|_{L^3} \leq C \| \overline p\|_{W^{1, 3} (\Omega)}$ then by the Gronwall Lemma (see [20]) we can conclude the result.

Theorem 3.8. Let be $\alpha = \alpha_0$ a constant function, the TAF concentration $v$ satisfies that $\partial_n v = 0$ on $\Gamma$, the function $\beta$ is twice differentiable with respect to the first variable and it verifies that $\partial ^2_{11}\beta$ is bounded and its norm in $L^{\infty}$ is small. Then, the optimality system (3.7) has a unique solution when the time $T$ is small enough.

Proof. Using the hypothesis $\alpha = \alpha_0$, the optimality system is given by these equations:

We call ${\mathcal T}$ to the following operator:

where $L^2(L^2)$ is the space $L^2(0, T; L^2(\Omega))$. The function $u$ is obtained by the data $z$:

Having $z$ and $u$, the function $p$ can be solved:

With this notation, the operator $T$ is

The existence of an optimal control is the existence of a fixed point of ${\mathcal T}$. This is guaranteed by the Theorem 3.1. We want to prove that ${\mathcal T}$ is contractive. Then, there will be a unique fixed point, $\hat{z}$, and we can say that the optimality system has a unique solution and so, the optimal control problem has a unique optimal control.

Let be $z_1, z_2\in L^2(L^2)$.

where we have called $u_i = {\mathcal T}_1(z_i)$ and $p_i = {\mathcal T}_2(z_i, {\mathcal T}_1(z_i))$, $i = 1, 2$. Adding and subtracting the appropriate terms and using the triangular inequality we have

By Theorem 2.2, it is known that $\Vert u\Vert_{L^{\infty}(L^{\infty})}$, where $u$ is a solution of (3.11), (3.12) and (3.13), and $\Vert p\Vert_{L^{\infty}(L^{\infty})}$, $p$ is a solution of (3.14), (3.15) and (3.16), are bounded.

Let $w = u_1-u_2$. Writing the equation for $w$, multiplying by $w$ and integrating in $\Omega$ we obtain the following estimate:

Applying the medium value theorem, $V'$, $\beta$, $\partial _1 \beta$ and $\Vert u_2\Vert_{L^{\infty}(Q)}$ are bounded and the Young's inequality, we have

Then,

Therefore,

We do a similar reasoning with $p_1-p_2$. Let $\eta = p_1-p_2$. Writing the equation for $\eta$, multiplying it by $\eta$ and integrating in $\Omega$, we get

We multiply the previous equality by $-1$ and we apply the Young inequality and the Hölder inequality to obtain

Hence, the Sobolev embedding and (3.21) entails

Next, we solve the differential equation to get

We apply (3.21) for $t = T$ to obtain

Therefore, after integration on the interval $[0, T]$ we have

After all of these inequalities we get the following equation

where $C(T)$ goes to zero as $T$ goes to zero. If we assume that $\Vert \partial ^2_{11}\beta\Vert_{L^{\infty}}$ is small enough we can say that the operator ${\mathcal T}$ is contractive if $T$ is small enough. Therefore, there exists a unique optimal control if $T$ is small.

4.   Numerical simulations

We are going to solve the optimality system in some particular cases. The program has been done in FreeFEM, using P1-Lagrange finit element method and Euler method to solve the problem of $u$ and $p$. In the case of $u$, because of the nonlinearity of the equation, we have used the Newton's method, and for the problem of $p$ we have to solve a backward equation. The optimality system is a recursive equation, by the hypothesis of Theorem 3.8, it is contractive for $T$ small enough, so, we have chosen $z$, we get $u$ and $p$, we obtain a new $z$ and we repeat until the difference between this one and the previous one satisfies the stop test.

The domain $\Omega$ is plotted in the following figure:

We have chosen the following data:

The coefficient for the chemotaxis term

the time step,

and a small final time, $T$,

We solve this problem to get $v$

with

The graphic of the TAF concentration is shown in Figure 1.

The initial density of endothelial cells is

and it is drawn in Figure 2:

The function $V$ which appears in the chemotaxis term is given by

and the function $\beta$ which is the growth rate of the drug $z$, is

we have taken $k = 1$. Finally, the coefficients, $a$ and $b$ are equal to one, that means that we optimize $u(T)$ and $z$ with the same weights. At the final time, the density of ECs is shown in Figure 3:

And the concentration of $z$ at the final time is plotted in Figure 4.

The functional and norms are

When the chemotaxis factor, $\alpha_0$, is smaller, the drug concentration is smaller too, it is necesary less chemical agent to control the angiogenesis process. We take

instead of 50 and the rest of the data the same as before. The density of ECs at the final time is shown in Figure 5.

In this case, the value of $J$ on the optimal and the norms of $u(T)$ and $z$ are

and the concentration of drug at the final times is drawn in Figure 6. If we prioritize to minimize the term of $u(T)$, taking

and we choose a growth function $\beta$ smaller than the previous one, taking

the functions $u(T)$ and $z(T)$ decrease notably, as we can see in Figure 7 and Figure 8.

Now, the functional and the norms are

Conflict of interest

The authors declare no conflict of interest.