Optimal control of a chemotaxis equation arising in angiogenesis

http://www.aimspress.com/journal/mine Mathematics in Engineering, 4(6): 1–25. DOI:10.3934/mine.2022047 Received: 03 February 2021 Accepted: 02 October 2021 Published: 04 November 2021 Research article Optimal control of a chemotaxis equation arising in angiogenesis† M. Delgado, I. Gayte and C. Morales-Rodrigo∗ Dpto. Ecuaciones Diferenciales y An´ alisis Num´ erico, Fac. Matem´ aticas, Univ. de Sevilla, Calle Tarfia s/n, 41012, Sevilla, Spain †This contribution is part of the Special Issue: Advances in the analysis of chemotaxis systems Guest Editor: Michael Winkler Link: www.aimspress.com/mine/article/6067/special-articles *Correspondence: Email: [email protected]. Abstract: In this paper we consider an optimal control for an equation that models a crucial step in the tumor development, the angiogenesis. We show the existence of an optimal control, we characterize the optimal control as a solution of the optimality system and we show the uniqueness of the optimal control for short times. Keywords: chemotaxis; anti-angiogenic therapy; optimal control 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 zthe concentration of 2 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,zand vare defined in a bounded domain Ω⊂IR3. Here we are interested in the minimum of the functional defined by J(u,z)=a 2ZΩ u(T)2+b 2ZΩ×(0,T) z2,(1.1) where aand bare positive parameters and (u,z) is a solution of the nonlinear differential equation            ut−div (∇u−V(u)α(z)∇v)=β(z,v)u−u2in Ω×(0,T), ∂nu−V(u)α(z)∂nv=0 on ∂Ω×(0,T), u(0,x)=u0(x) in Ω. (1.2) and z∈ Uad where Uad :={z∈L2(Ω×(0,T)) : z(x,t)≥0 for almost (x,t)∈Ω×(0,T)}. 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∈ Uad, is a problem of the optimal control theory. It is classical in this theory to call ustate variable and to call zcontrol 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∈ Uad which has empty interior in L2(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 Tsmall enough and this turns out the uniqueness of solution of the optimal problem. The chemotaxis is described in (1.2) by the term div(V(u)α(z)∇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 zbut with an additional equation for chemoattractant is given in [4]. The case of a therapy zsatisfying 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 Tis small enough are studied in Section 3. In Section 4 we show some numerical simulations to illustrate the theoretical results. Mathematics in Engineering Volume 4, Issue 6, 1–25. 3 2. Global existence and uniqueness for the equation Let Ω⊂IRNbe a regular bounded domain and T>0 a fixed number. We will denote Q= Ω ×(0,T) and Γ = ∂Ω×(0,T). In this paper, we consider the problem            ut−div (∇u−V(u)α(z)∇v)=β(z,v)u−u2in Q, ∂nu−V(u)α(z)∂nv=0 on Γ, u(0,x)=u0(x) in Ω. (2.1) Here, v∈W2,∞(Ω) is a known function, z∈L2(Q), and u0∈L2(Ω). The function V: [0,+∞)→ [0,+∞) verify V(0) =0 and V∈C2([0,∞)) ∩L∞([0,∞)). If it is needed we will extend this function by zero for negative values. The functions α: [0,+∞)→[0,+∞) is in W1,∞([0,+∞)) and β: [0,+∞)×[0,+∞)→[0,+∞) is in C2([0,∞)) ∩L∞([0,∞)). We would like to point out here that the main difficulty of the problem is that zis not defined on ∂Ω. Let T>0 and Xis a Banach space, we define the space Lp(0,T;X) of equivalence classes of measurable functions u:I→Xsuch that t∈(0,T)→ kukXbelongs to Lp(I) which is a Banach space for the norm kukLp(X)=           ZT 0kukp Xdt!1/p if 1 ≤p<∞ Ess supt∈(0,T)ku(t)kXif p=∞ For instance, we will use kukL2(L2)= ZT 0ZΩ|u(x)|2dx dt!1/2 , kukL2(H1)= ZT 0ZΩ|∇u(x)|2+|u(x)|2dx dt!1/2 . The definition of a weak solution for the problem is Definition 2.1. We say that uis a weak solution of (1.2) if the following conditions are verified 1) u∈W(0,T) where W(0,T) :={u∈L2(0,T;H1(Ω)) such that ut∈L2(0,T; (H1(Ω))0)}. 2) ∀w∈C1([0,T]ׯ Ω) : w(T,x)=0,∀x∈¯ Ω, it holds ZT 0hut,wi(H1)0,H1+ZT 0ZΩ∇u·∇w−ZT 0ZΩ V(u)α(z)∇v·∇w= =ZT 0ZΩ β(z,v)u w −ZT 0ZΩ u2w+ZΩ u0(x)w(0,x).(2.2) Theorem 2.2. If u0∈L∞(Ω), u0≥0, then there is a unique positive global weak solution to the problem (1.2) that also satisfies u∈L∞(0,T;L∞(Ω)) . Mathematics in Engineering Volume 4, Issue 6, 1–25. 4 Proof. We denote by BR={u∈L2(0,T;L2(Ω)) : kukL2(L2)≤R}, and we define the mapping S:BR→L2(0,T;L2(Ω)) such that for each ¯u∈BR,S(¯u)=uis the weak solution of the linear problem            ut=div (∇u−V(¯u)α(z)∇v)+β(z,v)u−Tk(¯u)uin Q, ∂nu−V(¯u)α(z)∂nv=0 on Γ, u(0,x)=u0(x) in Ω. (2.3) where Tk(ϕ)=(kif ϕ > k, ϕ+if ϕ≤k. We will prove the existence of a fixed point for Swhich 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). Step 1 First of all we show that the operator Sis 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])            ut=div (∇u)+PN i=1∂ifi(x,t)+c0(x,t)uin Q ∂nu+PN i=1fiνi=0 on Γ u(0,x)=u0(x) in Ω. (2.4) (where fi=−V(¯u)α(z)∂ivand c0(x,t)=β(z,v)−Tk(¯u) in (2.3) and ν=(ν1, .., νN) is the outward normal to Ω), has a unique weak solution i.e. u∈L2(0,T;H1(Ω)) ∩L∞(0,T;L2(Ω)) such that ZT 0ZΩ      −uwt+|∇u|·|∇w|− N X i=1 fi∂iw−c0(x,t)uw      dx dt =ZΩ u0(x)v(0,x)dx, for each w∈C1([0,T]ׯ Ω) such that w(T,x)=0,∀x∈Ω. Note that, if ut∈L2(0,T; (H1(Ω))0) then the previous definition it is exactly the one given in Definition (2.1). Therefore we should show that ut∈L2(0,T; (H1(Ω)0). hut, ϕi=ZΩ∇u·∇ϕdx dt +ZΩ V(¯u)α(z)∇v·∇ϕdx +ZΩ c0(x,t)uϕdx ≤ k∇ukL2k∇ϕkL2+kV(¯u)α(z)∇vkL2k∇ϕkL2+kc0(x,t)ukL2kϕkL2 ≤k∇ukL2+kV(¯u)α(z)∇vkL2+kc0(x,t)ukL2kϕkH1, where h·,·i stands for the duality pairing between (H1(Ω))0and H1(Ω). Therefore kutk(H1(Ω))0=sup 0,ϕ∈H1(Ω) hut, ϕi kϕkH1≤ k∇ukL2+kV(¯u)α(z)∇vkL2+kc0(x,t))ukL2, and kutkL2(0,T;(H1(Ω))0)= ZT 0kutk2 (H1(Ω))0dt!1/2 ≤C(2.5) Mathematics in Engineering Volume 4, Issue 6, 1–25. 5 i.e., ut∈L2(0,T; (H1(Ω))0). Hence, Sis well defined. In what follows we will apply the Schauder fixed point theorem to get the existence of a fixed point for Sin BR. Step 2 We claim that there exists R>0 and T>0 such that S(BR)⊂BR. In fact, multiplying the equation of (2.3) by uand integrating in Ω, it results 1 2 d dt ZΩ u2=−ZΩ|∇u|2+ZΩ α(z)V(¯u)∇u·∇v+ZΩ β(z)u2−ZΩ Tk(¯u)u2 ≤ −1 2ZΩ|∇u|2+1 2ZΩ α(z)2V(¯u)2|∇v|2+kβk∞ZΩ u2 ≤ kβk∞ZΩ u2+C. Thus for y(t)=ZΩ u2we obtain the problem (y0(t)≤ay(t)+b, y(0) =RΩu2 0, for some positive constants a,b. Consequently, 0≤y(t)≤ y0+b a!eat −b a. We can choose any R>y(0) =ku0k2and determine T>0 such that y(t)=ZΩ u2(x,t)dx ≤R,∀t∈(0,T),(2.6) and kukL2(L2)= ZT 0ZΩ|u(x,t)|2dx dt!1/2 ≤R1/2T1/2≤R, if we take R>1 and T<1. So, for R≥ ku0k2+1 there exists T<1 such that S(BR)⊂BR. On the other hand, taking into account that y0(t)+1 2ZΩ|∇u|2≤ay(t)+b, we can choose R, determine Tand integrate the inequality on the interval (0,T) to obtain y(T)−y(0) +ZT 0ZΩ|∇u|2dx dt ≤aZT 0 y dt +bZT 0 dt , ZT 0ZΩ|∇u|2dx dt ≤y(0) +aZT 0ZΩ u2dx dt +bT ≤ ku0k2+aR2+bT.(2.7) Consequently, kukL2(H1)is bounded. Step 3 We claim that Sis a continuous mapping. We will prove that kS(¯u1)−S(¯u2)kL2(L2)≤Φ(k,k¯u1−¯u2kL2(L2))∀¯u1,¯u2∈BR,(2.8) Mathematics in Engineering Volume 4, Issue 6, 1–25. 6 with Φ(k,s)→0 when s→0 for each fixed k, and Tas in Step 2. If we denote u1=S(¯u1) and u2=S(¯u2), it holds ((u1)t=div (∇u1−V(¯u1)α(z)∇v)+β(z,v)u1−Tk(¯u1)u1, (u2)t=div (∇u2−V(¯u2)α(z)∇v)+β(z,v)u2−Tk(¯u2)u2, and taking w=u1−u2=S(¯u1)−S(¯u2), we have wt= ∆w−div ((V(¯u1)−V(¯u2)α(z)∇v)+β(z,v)w+(Tk(¯u2)u2−Tk(¯u1)u1) = ∆w−div ((V(¯u1)−V(¯u2)α(z)∇v)+β(z,v)w+Tk(¯u1)w+ +(Tk(¯u2)−Tk(¯u1))u2. On multiplying the previous inequality by wand integrating on Ωwe obtain 1 2 d dt ZΩ w2=−ZΩ|∇w|2−ZΩ (V(¯u1)−V(¯u2))α(z)∇v·∇w+ +ZΩ β(z,v)w2+ZΩ Tk(¯u1)w2+ZΩ (Tk(¯u2)−Tk(¯u1))u2w. (2.9) Next, we provide a bound to the terms in the right hand side of (2.9) ZΩ (V(¯u1)−V(¯u2))α(z)∇v·∇w≤1 2ZΩ|∇w|2+1 2ZΩ (V(¯u1)−V(¯u2))2α(z)2|∇v|2 ≤1 2ZΩ|∇w|2+Ckαk2 ∞k∇vk2 ∞ZΩ|¯u1−¯u2|2, ZΩ β(z,v)w2≤ kβk∞ZΩ w2, ZΩ Tk(¯u1)w2≤kZΩ w2, ZΩ (Tk(¯u2)−Tk(¯u1))u2w≤1 2ZΩ w2+1 2ZΩ (Tk(¯u2)−Tk(¯u1))2u2 2. The Sobolev embedding H1(Ω),→L6(Ω) together with the estimate u2∈L2(0,T;H1(Ω)) implies that u2∈L2(0,T;L6(Ω)). On the other hand if we apply the H¨ older inequality with exponents 3 and 3/2 to the last term of the above inequality we obtain ZΩ (Tk(¯u2)−Tk(¯u1))2u2 2≤ ku2 2kL3/2k(Tk(¯u2)−Tk(¯u1))2kL3 =ku2k2 L3k(Tk(¯u2)−Tk(¯u1))k2 L6. Using the interpolation inequality, ku2k2 L3≤ ku2kL2ku2kL6≤C(R)ku2kL6, and, thus ZΩ (Tk(¯u2)−Tk(¯u1))2u2 2≤C(R)ku2kL6k(Tk(¯u2)−Tk(¯u1))k2 L6. Mathematics in Engineering Volume 4, Issue 6, 1–25. 7 Then (2.9) becomes 1 2 d dt ZΩ w2≤C1k¯u1−¯u2k2 L2+C2kwk2 L2+C(R)ku2kL6k(Tk(¯u2)−Tk(¯u1))k2 L6.(2.10) We denote by y(t)=ZΩ w(x,t)2dx. Since y(0) =0 it follows from (2.10) that y(t)≤ect ZT 0 e−cs C1k¯u1−¯u2k2 L2+C(R)ku2kL6k(Tk(¯u2)−Tk(¯u1))k2 L6ds ≤ecT C3k¯u1−¯u2kL2(L2)+ZT 0 C(R)e−csku2kL6k(Tk(¯u2)−Tk(¯u1))k2 L6ds!. (2.11) By the H¨ older inequality it holds ZT 0 C(R)e−csku2kL6k(Tk(¯u2)−Tk(¯u1))k2 L6ds ≤ ≤ ZT 0ku2k2 L6!1/2 ZT 0kTk(¯u2)−Tk(¯u1)k2 L63!1/3 ZT 0 (C(R)e−cs)6ds!1/6 . Since ZT 0kTk(¯u2)−Tk(¯u1)k2 L63!1/3 = ZT 0ZΩ|Tk(¯u2)−Tk(¯u1)|6dx ds!1/3 = = ZT 0ZΩ|Tk(¯u2)−Tk(¯u1)|4|Tk(¯u2)−Tk(¯u1)|2dx ds!1/3 ≤k4/3 ZT 0kTk(¯u2)−Tk(¯u1)k2 L2ds!1/3 ≤k4/3 ZT 0k¯u2−¯u1k2 L2ds!1/3 =k4/3k¯u2−¯u1k2/3 L2(L2), it results from (2.11), y(t)≤C4k¯u2−¯u1k2 L2(L2)+C5k4/3k¯u2−¯u1k2/3 L2(L2).(2.12) Integrating on (0,T), we obtain kS(¯u1)−S(¯u2)kL2(L2)≤T(C4k¯u2−¯u1k2 L2(L2)+C5k4/3k¯u2−¯u1k2/3 L2(L2)), which proves the desired continuity. Step 4 We claim that Sis a compact mapping in L2(Q). We know that for each u∈BRS(u)=uis bounded in L2(0,T;H1(Ω)) and S(u)t=utis bounded in L2(0,T; (H1(Ω))0) (see (2.7) and (2.5)) then by the Lions-Aubin Lemma (see for instance [18]) S(BR) is embedded compactly in L2(0,T;L2(Ω)). Step 5 By the Schauder Theorem we have a solution to the problem            ut=div (∇u−V(u)α(z)∇v)+β(z,v)u−Tk(u)uin Q, ∂nu−V(u)α(z)∂nv=0 on Γ, u(0,x)=u0(x) in Ω. (2.13) Mathematics in Engineering Volume 4, Issue 6, 1–25. 8 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 d 2dt ZΩ (u−)2≤0. Integrating on the space variable infer 0≤ZΩ u−(t)2≤ZΩ ((u0)−)2=0 Therefore, u−(t)=0 a.e. in Ω. Step 6 We prove that a solution to (2.3) satisfies u∈L∞(0,T;L∞(Ω)) independently of the k, as a consequence a solution of (2.13) is a solution to (1.2) for ksufficiently large. We multiply (2.3) by pTk(u)p−1for p≥2, a cut offfunction that approximate to pup−1(see for instance [3, p. 1196]) to get d dt ZΩ Tk(u)p+4(p−1) pZΩ|∇Tk(u)p/2|2=pZΩ β(z,v)up−pZΩ Tk(u)p+1+ +2ZΩ Tk(u)p/2−1V(u)∇(up/2). Since V(0) =0 and Vis a Lipschitz function then, d dt ZΩ Tk(u)p+2ZΩ|∇Tk(u)p/2|2≤ ≤pkβk∞ZΩ Tk(u)p+2ZΩ Tk(u)p/2−1|V(u)−V(0)||∇(Tk(u)p/2)| ≤pkβk∞ZΩ Tk(u)p+2CZΩ Tk(u)p/2|∇(Tk(u)p/2)| ≤pkβk∞ZΩ Tk(u)p+2C ZΩ Tk(u)p!1/2 ZΩ|∇(Tk(u)p/2)|2!1/2 Adding on both sides of the above inequality RΩTk(u)pwe deduce that d dt ZΩ Tk(u)p+ZΩ Tk(u)p+3 2ZΩ|∇(Tk(u)p/2)|2≤ kβk∞ p+1+2C2 kβk∞!ZΩ Tk(u)p.(2.14) At this point we will provide a bound for RΩup. We claim that for any  > 0 ZΩ Tk(u)p≤ZΩ|∇Tk(u)p/2|2+C() ZΩ Tk(u)p/2!2 .(2.15) Next lines are devoted to the proof of the above inequality. Let z=1 |Ω|RΩz. We know that ZΩ z2=ZΩ (z−z)2+1 |Ω| ZΩ z!2 . Mathematics in Engineering Volume 4, Issue 6, 1–25. 9 On the other hand, by the H¨ older inequality and the Poincare-Wintinger inequality we have ZΩ (z−z)2≤ kz−zk3/2 3kz−zk1/2 1 =kz−zk3/2 3√2kzk1/2 1≤0kz−zk2 3+C(0) ZΩ|z|!2 ≤ZΩ|∇z|2+C(0) ZΩ|z|!2 . Therefore if we take z=Tk(u)p/2the claim easily follows. By (2.15) we get from (2.14) that d dt ZΩ Tk(u)p+ZΩ Tk(u)p≤C(p+C) ZΩ Tk(u)p/2!2 . For any 0 ≤t≤T<Tmax, where Tmax stand for the maximal existence time, the above inequality asserts ZΩ Tk(u)p(t)≤ZΩ up 0e−t+C(p+C) sup 0≤t≤T ZΩ Tk(u)p/2!2 ≤C(p+C) max       ku0kp ∞,sup 0≤t≤T ZΩ Tk(u)p/2!2       . (2.16) Let us define θ(p/2) :=max       ku0k∞,sup 0≤t≤T ZΩ Tk(u)p/2!2/p       . From (2.16), by a recursive procedure, we have ZΩ Tk(u)p(t)!1/p ≤C(p+k)1/pθ(p/2) ≤C(p+C)1/pC(p/2+C)2/pθ(p/4) In particular, for p=2jwe deduce ZΩ Tk(u)2j(t)!2−j ≤Csj j Y i=0 (2i+C)2−iθ(1), where sj=Pj i=02−i. Taking into account that θ(1) ≤max (ku0k∞,sup 0≤t≤TZΩ u):=C1, Therefore we can take k→+∞to get ZΩ u2j(t)!2−j ≤Csj j Y i=0 (2i+C)2−iC1 Mathematics in Engineering Volume 4, Issue 6, 1–25. 16 u(x,0) =u0(x) in Ω,(3.13) −pt−∆p=α0V0(u)∇v·∇p+β(z,v)p−2up in Ω×(0,T),(3.14) ∂np=0 on ∂Ω×(0,T),(3.15) p(T)=u(T) in Ω,(3.16) z=[−a b(∂1β(z,v)up)]+.(3.17) Lemma 3.6. We have that ku(T)kW1,3(Ω)≤C. Proof. Using the variations of constants formula in (3.11)–(3.13) we get u(t)=e−tAu0+Zt 0 e−(t−s)Ah(u,v)ds ,(3.18) where A=−∆ + Iwith domain {ψ∈W2,3(Ω) : ∂nu=0}and h(u,v) :=−∇·(α0V(u)∇v)+(β(z,v)+1)u−u2. Let Xβ,β∈(0,1) the fractional powers of X. By [8, Theorem 1.6.1] we have that Xβ,→W1,3, for any β∈(1/2,1). Taking the norm W1,3in (3.18) and applying the above embedding in the right hand side we infer ku(T)kW1,3(Ω)≤ ke−tAu0kXβ+ZT 0ke−(t−s)AhkXβ. By [8, Theorem 1.4.3] we deduce ku(T)kW1,3(Ω)≤Cβe−δtt−βku0kL3+ZT 0 e−δ(t−s)(t−s)−βkh(u,v)kL3 Let us notice that by the regularity of β,Vand vand the embedding W1,3(Ω),→L6(Ω) we have kh(u,v)kL3≤CkukW1,3(Ω). 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 p(t)= p(T−t). Therefore, the problem is            pt−∆p=α0V0(u)∇v·∇p+β(z,v)p−2up in Ω×(0,T), ∂np=0 on ∂Ω×(0,T), p(0) =u(T) in Ω. (3.19) Lemma 3.7. We have that max t∈[0,T]kp(t)kW1,3(Ω)≤C. Mathematics in Engineering Volume 4, Issue 6, 1–25. 17 Proof. As in the previous Lemma we apply the variation of constants formula for pto get p(t)=e−tAu(T)+Zt 0 e−(t−s)Af(u,v,p)ds , where f(u,v,p) :=α0V0(u)∇v·∇p+(β(z,v)+1)p−2up . We take W1,3(Ω) in the variations of constants formula for pand we apply [8, Theorem 1.4.3, Theorem 1.6.1] and ( [10, p. 59]) to infer kp(t)kW1,3(Ω)≤Cku(T)kW1,3(Ω)+ZT 0 (t−s)−βe−δ(t−s)kf(u,v,p)kL3ds with β∈(1/2,1). Since kf(u,v,p)kL3≤CkpkW1,3(Ω)then by the Gronwall Lemma (see [20]) we can conclude the result.  Theorem 3.8. Let be α=α0a constant function, the TAF concentration v satisfies that ∂nv=0on Γ, the function βis twice differentiable with respect to the first variable and it verifies that ∂2 11βis bounded and its norm in L∞is small. Then, the optimality system (3.7) has a unique solution when the time T is small enough. Proof. Using the hypothesis α=α0, the optimality system is given by these equations: z=[−a b(∂1β(z,v)up)]+.(3.20) We call Tto the following operator: T:L2(L2)→L2(L2) z7→ T(z)=[−a b∂1β(z,v)up]+, where L2(L2) is the space L2(0,T;L2(Ω)). The function uis obtained by the data z: T1:L2(L2)→W(0,T) z7→ u,the solution of (3.11), (3.12), (3.13) Having zand u, the function pcan be solved: T2:L2(L2)×W(0,T)→W(0,T) (z,T1(z)) 7→ p,the solution of (3.14), (3.15), (3.16) With this notation, the operator Tis T(z)=[−a b∂1β(z,v)T1(z)T2(z,T1(z))]+. The existence of an optimal control is the existence of a fixed point of T. This is guaranteed by the Theorem 3.1. We want to prove that Tis contractive. Then, there will be a unique fixed point, ˆz, and Mathematics in Engineering Volume 4, Issue 6, 1–25. 18 we can say that the optimality system has a unique solution and so, the optimal control problem has a unique optimal control. Let be z1,z2∈L2(L2). kT(z1)−T(z2)kL2(L2)≤a bk∂1β(z1,v)u1p1−∂1β(z2,v)u2p2kL2(L2), where we have called ui=T1(zi) and pi=T2(zi,T1(zi)), i=1,2. Adding and subtracting the appropriate terms and using the triangular inequality we have kT(z1)−T(z2)kL2(L2)≤a b[k∂2 11βkL∞(L∞)kz1−z2kL2(L2)ku1kL∞(L∞)kp1kL∞(L∞)+ +k∂1βkL∞(L∞)ku1−u2kL2(L2)kp1kL∞(L∞)+k∂1βkL∞(L∞)kp1−p2kL2(L2)ku2kL∞(L∞)]. By Theorem 2.2, it is known that kukL∞(L∞), where uis a solution of (3.11), (3.12) and (3.13), and kpkL∞(L∞),pis a solution of (3.14), (3.15) and (3.16), are bounded. Let w=u1−u2. Writing the equation for w, multiplying by wand integrating in Ωwe obtain the following estimate: 1 2 d d tkwk2 L2+1 2k∇wk2 L2≤1 2ZΩ α2 0(V(u1)−V(u2))2|∇v|2+ZΩ β(z1,v)w2+ +ZΩ (β(z1,v)−β(z2,v))u2w−ZΩ (u1+u2)w2. Applying the medium value theorem, V0,β,∂1βand ku2kL∞(Q)are bounded and the Young’s inequality, we have d d tkwk2 L2≤Akwk2 L2+Bkz1−z2k2 L2,A,B>0. Then, kw(t)k2 L2≤BeAtkz1−z2k2 L2(L2).(3.21) Therefore, ku1−u2k2 L2(L2)≤B A(eAT −1)kz1−z2k2 L2(L2). We do a similar reasoning with p1−p2. Let η=p1−p2. Writing the equation for η, multiplying it by ηand integrating in Ω, we get −d dt ZΩ η2+ZΩ|∇η|2=α0ZΩ V0(u1)∇v·∇η η +α0ZΩ (V0(u1)−V0(u2))∇v·∇p2η+ +ZΩ β(z1,v)η2+ZΩ (β(z1,v)−β(z2,v))p2η−2ZΩ u1η2−2ZΩ wp2η. (3.22) We multiply the previous equality by −1 and we apply the Young inequality and the H¨ older inequality to obtain d dt ZΩ η2+ZΩ|∇η|2≤ZΩ|∇η|2+C()ZΩ η2+ +k(u1−u2)∇vkL2k∇p2kL3kηkL6+Ckz1−z2k2 L2+CZΩ w2. Mathematics in Engineering Volume 4, Issue 6, 1–25. 19 Hence, the Sobolev embedding and (3.21) entails d dt ZΩ η2≤C()ZΩ η2+C()BeAT kz1−z2kL2(L2)+Ckz1−z2kL2. Next, we solve the differential equation to get ZΩ η2≤eC()tku1(T)−u2(T)kL2(Ω)+ec()tC()BeAT kz1−z2kL2(L2)C(T)+ +eC()tZt 0 e−C()sCkz1−z2k2 L2. We apply (3.21) for t=Tto obtain ZΩ η2≤eC()tC(T)kz1−z2kL2(L2). Therefore, after integration on the interval [0,T] we have kηkL2(L2)≤C(T)eC()T−1 C()kz1−z2kL2(L2). After all of these inequalities we get the following equation kT(z1)−T(z2)kL2(L2)≤C(C(T)+k∂2 11βkL∞ku1kL∞kp1kL∞)kz1−z2kL2(L2), where C(T) goes to zero as Tgoes to zero. If we assume that k∂2 11βkL∞is small enough we can say that the operator Tis contractive if Tis small enough. Therefore, there exists a unique optimal control if Tis 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 uand 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 pwe have to solve a backward equation. The optimality system is a recursive equation, by the hypothesis of Theorem 3.8, it is contractive for Tsmall enough, so, we have chosen z, we get uand p, we obtain a new zand we repeat until the difference between this one and the previous one satisfies the stop test. The domain Ωis plotted in the following figure: Mathematics in Engineering Volume 4, Issue 6, 1–25. 20 We have chosen the following data: The coefficient for the chemotaxis term α0=50, the time step, dt =0.1 and a small final time, T, T=0.5. We solve this problem to get v −∆v+v=g,Ω ∂v ∂n=0, ∂Ω with g=1 0.01 +x2+y2. The graphic of the TAF concentration is shown in Figure 1. Figure 1. TAF concentration. Mathematics in Engineering Volume 4, Issue 6, 1–25. 21 The initial density of endothelial cells is u0=x2+y2+100, and it is drawn in Figure 2: Figure 2. Density of ECs at the initial time. The function Vwhich appears in the chemotaxis term is given by V(u)=u 1+u2, and the function βwhich is the growth rate of the drug z, is β(z,v)=k exp(−z)v 1+v2, we have taken k=1. Finally, the coefficients, aand bare equal to one, that means that we optimize u(T) and zwith the same weights. At the final time, the density of ECs is shown in Figure 3: Figure 3. Density of ECs at the final time. Mathematics in Engineering Volume 4, Issue 6, 1–25. 22 Figure 4. Concentration of the drug at the final time. And the concentration of zat the final time is plotted in Figure 4. The functional and norms are J(ˆu,ˆz)=65.1,ku(T)kL2(Ω)=6.27,kzkL2(Q)=9.53 When the chemotaxis factor, α0, is smaller, the drug concentration is smaller too, it is necesary less chemical agent to control the angiogenesis process. We take α0=5, 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. Figure 5. Density of ECs at the final time with α0=5. Mathematics in Engineering Volume 4, Issue 6, 1–25. 23 Figure 6. Concentration of the drug at the final time with α0=5. In this case, the value of Jon the optimal and the norms of u(T) and zare J(ˆu,ˆz)=66.46,ku(T)kL2(Ω)=5.7,kzkL2(Q)=10, 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 a=10,b=1 and we choose a growth function βsmaller than the previous one, taking k=0.1 the functions u(T) and z(T) decrease notably, as we can see in Figure 7 and Figure 8. Figure 7. Density of ECs at the final time with α0=5 and k=0.1. Mathematics in Engineering Volume 4, Issue 6, 1–25. 24 Figure 8. Concentration of the drug at the final time with α0=5 and k=0.1 Now, the functional and the norms are J(ˆu,ˆz)=70,ku(T)kL2(Ω)=3,76,kzkL2(Q)=0,46. Conflict of interest The authors declare no conflict of interest. References 1. A. R. A. Anderson, M. A. J. Chaplain, Continuous and discrete mathematical models of tumorinduced angiogenesis, Bull. Math. Biol.,60 (1998), 857–899. 2. E. P. Avakov, Necessary extremum conditions for smooth abnormal problems with equality and inequality-type constraints, Mathematical Notes of the Academy of Sciences of the USSR,45 (1989), 431–437. 3. P. Biler, W. Hebisch, T. Nadzieja, The Debye system: existence and large time behavior of solutions, Nonlinear Anal.,23 (1994), 1189–1209. 4. M. Delgado, I. Gayte, C. Morales-Rodrigo, A. Su´ arez, An angiogenesis model with nonlinear chemotactic response and flux at the tumor boundary, Nonlinear Anal.,72 (2010), 330–347. 5. M. Delgado, C. Morales-Rodrigo, A. Su´ arez, Anti-angiogenic therapy based on the binding to receptors, DCDS,32 (2012), 3871–3894. 6. I. Gayte, F. Guillen-Gonz´ alez, M. Rojas-Medar, Dubovitskii-Milyutin formalism applied to optimal control problems with constraints given by the heat equation with final data, IMA J. Math. Control Inf.,27 (2010), 57–76. 7. I. V. Girsanov, Lectures on mathematical theory of extremum problem, Springer-Verlag, 1970. 8. D. Henry, Geometric theory of semi linear parabolic equations, Berlin-New York: Springer-Verlag, 1981. Mathematics in Engineering Volume 4, Issue 6, 1–25. 25 9. T. Hillen, K. Painter, Global existence for a parabolic chemotaxis model with prevention of overcrowding, Adv. Appl. Math.,26 (2001), 280–301. 10. D. Horstmann, M. Winkler, Boundedness vs. blow-up in a chemotaxis system, J. Differ. Equations, 215 (2005), 52–107. 11. W. J¨ ager, S. Luckhaus, On explosions of solutions to a system of partial differential equations modelling chemotaxis, Trans. Amer. Math. Soc.,329 (1992), 819–824. 12. E. F. Keller, L. A. Segel, Initiation of slime mold aggregation viewed as an instability, J. Theor. Biol.,26 (1970), 399–415. 13. U. Ledzewicz, On distributed parameter control systems in the abnormal case and in the case of nonoperator equality constraints, International Journal of Stochastic Analysis,6(1993), 704189. 14. G. M. Lieberman, Second order parabolic differential equations, River Edge, NJ: World Scientific Publishing Co., Inc., 1996. 15. J. L. Lions, Quelques m´ethodes de r´esolution des probl`emes aux limites non lin´eaires, (French), Paris: Dunod Gauthier-Villars, 1969. 16. N. V. Mantzaris, S. Webb, H. G. Othmer, Mathematical modeling of tumor induced angiogenesis, J. Math. Biol.,49 (2004), 111–187. 17. C. Morales-Rodrigo, A therapy inactivating the tumor angiogenic factors, Math. Biosci. Eng.,10 (2013), 185–198. 18. J. Simon, Compact sets in the space Lp(0,T;B), Ann. Mat. Pura Appl.,146 (1987), 65–96. 19. S. Walczak, Some properties of cones in normed spaces and their application to investigating extremal problems, J. Optim. Theory Appl.,42 (1984), 561–582. 20. C.-L. Wang, A short proof of a Greene theorem, Proc. Amer. Math. Soc.,69 (1978), 357–358. c 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0) Mathematics in Engineering Volume 4, Issue 6, 1–25.