Adjoint method for a tumor invasion PDE-constrained optimization problem using FEM

0.27 MB·

by  Andrés Quiroga

Adjoint method for a tumor invasion PDE-constrained optimization
problem using FEM
A. A. I. Quirogaa,b,∗, D. R. Ferna´ndeza,b, G. A. Torresa,b, C. V. Turnera,b
aFacultad de Matem´atica, Astronom´ıa y F´ısica, Medina Allende s/n, 5000 Co´rdoba, Argentina
bCentro de Investigaciones y Estudios enMatem´atica - CONICET, Medina Allende s/n, 5000 Co´rdoba, Argentina
4
1 Abstract
0
In this paper we present a method for estimating unknown parameter that appear on a non-linear
2
reaction-diffusion model of cancer invasion. This model considers that tumor-induced alteration of
n
micro-enviromentalpHprovidesamechanismforcancerinvasion. Acoupledsystemreaction-diffusion
a
J describing this model is given by three partial differential equations for the non dimensional spatial
2 distribution and temporal evolution of the density of normal tissue, the neoplastic tissue growth and
1 the excess concentration of H+ ions. Each of the model parameters has a corresponding biological
interpretation, for instance, the growth rate of neoplastic tissue, the diffusion coefficient, the reab-

] sorption rate and the destructive influence of H+ ions in the healthy tissue.
C
After solving the forward problem properly, we use the model for the estimation of parameters
O
by fitting the numerical solution with real data, obtained via in vitro experiments and fluorescence
.
h ratioimagingmicroscopy. We define anappropriatefunctionaltocompareboth the realdata andthe
t numerical solution using the adjoint method for the minimization of this functional.
a
m We applyFinite ElementMethod(FEM)to solveboththe directandinverseproblem,computing
the a posteriori error.
[

1 Keywords: reaction-diffusion equation, tumor invasion, PDE-constrained optimization, adjoint
v method, Finite Element Method, a posteriori error
5
2
6
1. Introduction.
2
.
1
Cancer is one of the greatest killers in the world although medical activity has been successful,
0
despite great difficulties, at least for some pathologies. A great effort of human and economical
4
1 resources is devoted, with successful outputs, to cancer research, [1, 2, 3, 4, 5, 6].
: Some comments on the importance of mathematical modeling in cancer can be found in the
v
literature. In the work [4] the authors say “Cancer modelling has, over the years, grown immensely
i
X as one of the challenging topics involving applied mathematicians working with researchers active in
r the biological sciences. The motivation is not only scientific as in the industrial nations cancer has
a
now movedfromseventhto secondplace in the league table offatal diseases,being surpassedonly by
cardiovascular diseases.”
We use in this work the mathematical analyses first proposed by [7] which supports the acid-
mediated invasion hypothesis, hence it is acquiescent to mathematical representation as a reaction-
diffusion system at the tissue scale, describing the spatial distribution and temporal development of
tumor tissue, normal tissue, and excess H+ ion concentration.
ThemodelpredictsapHgradientextendingfromthetumor-hostinterface. Theeffectofbiological
parameterscriticaltocontrollingthistransitionissupportedbyexperimentalandclinicalobservations
[8].
In [7] the authors model tumor invasion in an attempt to find a common, underlying mechanism
by which primary and metastatic cancers invade and destroy normal tissues. They are not modeling
∗Correspondingauthor: FacultaddeMatema´tica,Astronom´ıayF´ısica,MedinaAllendes/n,5000Co´rdoba,Argentina
Email address: [email protected] (A.A.I.Quiroga)
Preprint submitted to Elsevier January 14, 2014
the genetic changeswhichresultintransformationnordothey seekto understandthe causesofthese
changes. Similarly, they do not attempt to model the large-scale morphological features of tumors
such as central necrosis. Rather, they concentrate on the microscopic scale population interactions
occurring at the tumor-host interface, reasoning that these processes strongly influence the clinically
significant manifestations of invasive cancer.
Specifically, the authors hypothesize that transformation-induced reversion of neoplastic tissue to
primitive glycolytic metabolic pathways, with resultant increased acid production and the diffusion
of that acid into surrounding healthy tissue, creates a peritumoral microenvironment in which tumor
cellssurviveandproliferate,whereasnormalcellsareunabletoremainviable. Thefollowingtemporal
sequence would derive: (a) high H+ ion concentrations in tumors will extend, by chemical diffusion,
as a gradient into adjacent normal tissue, exposing these normal cells to tumor-like interstitial pH;
(b) normal cells immediately adjacent to the tumor edge are unable to survive in this chronically
acidic environment; and (c) the progressive loss of layers of normal cells at the tumor-host interface
facilitates tumor invasion. Key elements of this tumor invasion mechanism are low interstitial pH
of tumors due to primitive metabolism and reduced viability of normal tissue in a pH environment
favorable to tumor tissue.
These model equations depend only on a small number of cellular and subcellular parameters.
Analysis of the equations shows that the model predicts a crossoverfrom a benign tumor to one that
is aggressively invasive as a dimensionless combination of the parameters increases through a critical
value.
The dynamics and structure of the tumor-host interface in invasive cancers are shown to be con-
trolledbythesamebiologicalparameterswhichgeneratethetransformationfrombenigntomalignant
growth. A hypocellular interstitial gap, as we can see in Figure 1 [7, Figure 4a], at the interface is
predicted to occur in some cancers.
Figure1: Amicrographsofthetumor-hostinterfacefromhumansquamous cellcarcinomasoftheheadandneck[7].
In this paper we estimate one of these parameters (the destructive influence of H+ ions in the
healthy tissue) using an inverse problem. Moreover, via fluorescence ratio imaging microscopy, it is
possible get data about the concentration of hydrogen ions [8]. We propose a framework via a PDE-
constrained optimization problem, following the PDE-based model by Gatenby [7]. In this approach,
tumorinvasionismodeledviaacouplednonlinearsystemofpartialdifferentialequations,whichmakes
the numerical solution procedure quite challenging.
This kind of problem constitutes a particular application of the so-called inverse problems, which
are being increasingly used in a broad number of fields in applied sciences. For instance, problems
referred to structured population dynamics [9], computerized tomography and image reconstruction
in medical imaging [10, 11], and more specifically tumor growth [12, 13, 14], among many others.
Wesolveaminimizationproblemusingagradient-basedmethodconsideringtheadjointmethodin
ordertofindthederivativeofanobjectivefunctional. Inthisway,wewouldobtainthebestparameter
that fits patient-specific data.
The contents of this paper, which is organized into 9 sections and an Appendix, are as follows:
Section 2 consists in some preliminaries about the model and the definition of the direct problem.
2
Section 3 deals with the variationalformulationofthe direct problem. Section4 considersthe formu-
lation of the minimization problem. Section 5 introduces the reduced and adjoint problem, deriving
theoptimalityconditionsfortheproblem. Section6findsthederivativeofthesolutionofafunctional
with respect to a parameter that does not appear explicitly in the equation. Section 7 deals with the
numericalsolutionoftheadjointproblem,designingasuitablealgorithmtosolveit. Inparticular,we
use the Finite Element Method with a computation of a posteriori error. In Section 8 we show some
numericalsimulationstogiveinformationonthebehaviorofthefunctionalanditsdependence onthe
parameters including the corresponding tables. Section 9 presents the conclusions and some future
work related to the contents of this paper. In the Appendix we include all algebraics of Section 6.
Some words about our notation. We use h·,·i to denote the L2 inner product (the space is always
clear from the context) and we consider the sum of inner products for a cartesian product of spaces.
For a function F :V ×U →Z such that (u,δ )7→F(u,δ ), we denote by F′(u,δ ) the full Fr´echet-
ad 1 1 1
derivative and by ∂F(u,δ ) and ∂F(u,δ ) the partial Fr´echet-derivativesof F at (u,δ ). For a linear
operator T : V →∂Zu we d1enote T∂δ∗1 : Z∗1→ V∗ the adjoint operator of T. If T is inv1ertible, we call
T−∗ the inverse of the adjoint operator T∗.
2. Some preliminaries about the model.
We present the mathematical model of the tumor-host interface based on the acid mediation
hypothesis of tumor invasion due to [7]. For convenience we reproduce the equations here, which
determine the spatial distribution and temporal evolution of three fields: N (x,t), the density of
1
normal tissue; N (x,t), the density of neoplastic tissue; and L(x,t), the excess concentration of H+
2
ions. TheunitsofN andN arecells/cm3 andexcessH+ ionconcentrationisexpressedasamolarity
1 2
(M), x and t are the position (in cm) and time (in seconds), respectively.
∂N N
1 1
= r N 1− −d LN , (1)
1 1 1 1
∂t (cid:18) K (cid:19)
1
∂N N N
2 2 1
= r N 1− +∇· D 1− ∇N , (2)
∂t 2 2(cid:18) K (cid:19) (cid:18) N2(cid:18) K (cid:19) 2(cid:19)
2 1
∂L
= r N −d L+D ∆L, (3)
∂t 3 2 3 N3
where the variables are in Ω×[0,T].
In equation (1) the behavior of the normal tissue is determined by the logistic growth of N with
1
growth rate r and carrying capacity K , and the interaction of N with excess H+ ions leading to a
1 1 1
deathrateproportionaltoL. The numberd Lis the excessacidconcentration,dependentdeathrate
1
inaccordwiththewell-describeddeclineinthegrowthrateofnormalcells,duetothereductionofpH
from its optimal value of 7.4. The constants r , d and K have units of 1/s, l/(M s) and cells/cm3,
1 1 1
respectively.
For equation (2), the neoplastic tissue growth is described by a reaction-diffusion equation. The
reaction term is governed by a logistic growth of N with growth rate r and carrying capacity K .
2 2 2
ThediffusiontermdependsontheabsenceofhealthytissuewithadiffusionconstantD . Constants
N2
r , K and D have units of 1/s, cells/cm3 and cm2/s, respectively.
2 2 N2
In equation (3), it is assumed that excess H+ ions are produced at a rate proportional to the
neoplastic cell density, and diffuse chemically. An uptake term is included to take account of the
mechanisms for increasing local pH (e.g., buffering and large-scalevascular evacuation [7]). Constant
r is the production rate (M cm3/(cell s)), d is the reabsorption rate (1/s), and D is the H+ ion
3 3 N3
diffusion constant (cm2/s).
All the parameter values can be found in Table 1.
2.1. Nondimensionalization.
Followingtheideasexposedin[7],andconsideringthatΩ⊂ ,themathematicalmodelisrescaled
R
and the domain is transformed onto the interval [0,1]×[0,T]. Hence, let us define the following
functions:
3
Parameter Estimate
K 5×107/cm3
1
K 5×107/cm3
2
r 1×10−6/s
1
r 1×10−6/s
2
D 2×10−10cm2/s
N2
D 5×10−6cm2/s
N3
r 2.2×10−17M cm3/s
3
d 1.1×10−4/s
3
Table1: Parametervaluesusedin[7].
N N
1 2
u = u =
1 2
K K
1 2
L
u3 = τ = r1t (4)
L
0
r
1
ξ = x
rDN3
where L = r K /d . We will continue denoting x and t instead of ξ and τ, respectively. Using the
0 3 2 3
transformation (4) the dimensionless form of the equations (1)-(3) become
∂u
1
= u (1−u )−δ u u , (5)
1 1 1 1 3
∂t
∂u ∂ ∂u
2 2
= ρ u (1−u )+ D (1−u ) , (6)
2 2 2 2 1
∂t ∂x(cid:18) ∂x (cid:19)
∂u ∂2u
3 3
= δ (u −u )+ , (7)
∂t 3 2 3 ∂x2
for (x,t) ∈ (0,1)×(0,T], where the four dimensionless quantities which parameterize the model are
given by:
d r K r D d
δ = 1 3 2, ρ = 2, D = N2, δ = 3.
1 2 2 3
d r r D r
3 1 1 N3 1
The interaction parameters between different cells (healthy and tumor) and concentration of H+
are difficult to measure experimentally. This is the reasonfor which we propose to estimate them, so
we will focus on δ in this work. The other parameters can be estimated by different techniques (see
1
Table 1).
2.2. Initial and boundary conditions.
Att=0wewillconsiderthe tumorata certainstageofitsevolution. Hence the initialconditions
are:
u (x,0) = u0(x), (8)
1 1
u (x,0) = u0(x), (9)
2 2
u (x,0) = u0(x), (10)
3 3
for all x∈[0,1]. We assume that the tumor is on the left of the domain, in the sense that the tumor
cells are not moving. Then, for all t∈[0,T], we have
∂u
1
(0,t)=0, u (1,t)=1, (11)
1
∂x
4
∂u
2
(0,t)=0, u (1,t)=0, (12)
2
∂x
∂u
3
(0,t)=0, u (1,t)=0. (13)
3
∂x
From now on, equations (5)-(13) will be referred to as the direct problem.
3. Variational form for the direct problem.
Using the variationaltechniques for obtainingthe weaksolutionofthe directproblem[15, 16,17],
we define the following weak formulation:
T 1 ∂u
1
0 = λ −u (1−u )+δ u u dxdt+
1 1 1 1 1 3
Z Z (cid:20) ∂t (cid:21)
0 0
T 1 ∂u ∂ ∂u
2 2
λ −ρ u (1−u )− D (1−u ) dxdt+
2 2 2 2 2 1
Z Z (cid:20) ∂t ∂x(cid:18) ∂x (cid:19)(cid:21)
0 0
T 1 ∂u ∂2u
3 3
λ −δ (u −u )− dxdt, (14)
Z Z 3(cid:20) ∂t 3 2 3 ∂x2 (cid:21)
0 0
where λ=(λ ,λ ,λ ),
1 2 3
∂v
λ ,λ ,λ ∈W = v ∈L2(0,T;H1((0,1))) and ∈L2(0,T;(H1((0,1)))∗) ,
1 2 3 (cid:26) D ∂t D (cid:27)
L2(0,T;H1((0,1)))= v(x,·)∈L2((0,T)) and v(·,t)∈H1((0,1))
D D
(cid:8) (cid:9)
and
H1 = v ∈H1((0,1)):v =0 on Γ ={1} .
D D
(cid:8) (cid:9)
Using integration by parts and boundary condition for λ and u in (14) we get the following weak
formulation of (5)-(13):
T 1 ∂u
1
0 = λ −u (1−u )λ +δ u u λ dxdt+
1 1 1 1 1 1 3 1
Z Z (cid:18) ∂t (cid:19)
0 0
T 1 ∂u ∂u ∂λ
2 2 2
λ −ρ u (1−u )λ +D (1−u ) dxdt+
2 2 2 2 2 2 1
Z Z (cid:18) ∂t ∂x ∂x (cid:19)
0 0
T 1 ∂u ∂u ∂λ
3 3 3
λ +δ u λ −δ u λ + dxdt. (15)
3 3 3 3 3 2 3
Z Z (cid:18) ∂t ∂x ∂x (cid:19)
0 0
A weak solution u = [u ,u ,u ]T ∈ V = W3 is a function that satisfies (15) for all λ ∈ V and
1 2 3
u(x,0)=u0(x)=[u0(x),u0(x),u0(x)].
1 2 3
4. Formulation of the minimization problem.
As described above we propose to use an inverse problem technique in order to estimate δ .
1
Function u represents the solution of the direct problem (the components of u are the state variables
of the problem) for each choice of the parameter δ .
1
Letusassumethatexperimentalinformationisavailableduringthetimeinterval0≤t≤T. Then,
the inverse problem can be formulated as:
Find a parameter δ able to generate data u=[u ,u ,u ]T that best match the available
1 1 2 3
(experimental) information over time 0≤t≤T.
5
Figure2: AmapofperitumoralH+flowusingvectorsgeneratedfromthepHdistributionaroundthetumor,[18,Figure
4].
For this purpose, we should construct an objective functional which gives us a notion of distance
between the experimental (real) data and the solution of the system of PDEs for each choice of the
parameter δ .
1
First of all, it is important to decide which variables are capable to be measured experimentally.
For instance, the excess concentration of H+ ions can be measured using fluorescence ratio imaging
microscopy [8, 18] at certain times t , k = 1,...,M. For example, Figure 2 [18, Figure 4] shows
k
a map of peritumoral H+ flow using vectors generated from the pH distribution around the tumor.
Such experiments could help to determine optimal variables and the parameter in order to control
real tumor invasion.
So, the functional J :V ×U →R could be defined as:
ad
1 T 1
J(u,δ )= [u (x,t)−uˆ (x,t)]2dxdt, (16)
1 3 3
2Z Z
0 0
where u (x,t) is the excess concentration of H+ ions obtained by solving the direct problem for a
3
certain choice of δ and uˆ (x,t) is the excess concentration measured experimentally (real data).
1 3
Let us define E :V ×U →V∗×Z∗ such that
ad
T 1 ∂u
1
hE(u,δ ),ζi = λ −u (1−u )λ +δ u u λ dxdt+
1 1 1 1 1 1 1 3 1
Z Z (cid:18) ∂t (cid:19)
0 0
T 1 ∂u ∂u ∂λ
2 2 2
λ −ρ u (1−u )λ +D (1−u ) dxdt+
2 2 2 2 2 2 1
Z Z (cid:18) ∂t ∂x ∂x (cid:19)
0 0
T 1 ∂u ∂u ∂λ
3 3 3
λ +δ u λ −δ u λ + dxdt+
3 3 3 3 3 2 3
Z Z (cid:18) ∂t ∂x ∂x (cid:19)
0 0
1 1
(u (x,0)−u0(x))γ dx+ (u (x,0)−u0(x))γ dx+
1 1 1 2 2 2
Z Z
0 0
1
(u (x,0)−u0(x))γ dx
3 3 3
Z
0
∂u
= ,λ +hF(u),λi + u(x,0)−u0(x),γ , (17)
(cid:28)∂t (cid:29) V∗,V Z∗,Z
V∗,V (cid:10) (cid:11)
where ζ =[λ,γ], γ =[γ ,γ ,γ ]∈Z and Z = H1((0,1)) 3.
1 2 3 D
In this way we can rewrite the weak formu(cid:0)lation (15) a(cid:1)s E(u,δ1)=0.
6
The parameter that best matches the experimental information with the generateddata provided
by the directproblemcanbe computed by solvinga PDE-constrainedoptimizationproblem, namely:
minimize J(u,δ )
1
δ1
subjectto E(u,δ )=0, (18)
1
δ ∈U ,
1 ad
where U denotes the set of admissible values of δ . In our case, U should be a subset of (0,∞).
ad 1 ad
Notice that a solution (u,δ ) must satisfy the constraint E(u,δ ) = 0, which constitutes the direct
1 1
problem.
We remark that, in general, there is a fundamental difference between the direct and the inverse
problems. In fact, the latter is usually ill-posed in the sense of existence, uniqueness and stability of
the solution. This inconvenient is often treated by using some regularizationtechniques [10, 19, 20].
5. Formulation of the reduced and adjoint problems.
In the following, we will consider the so-called reduced problem
minimize J˜(δ )=J(u(δ ),δ )
1 1 1
δ1 (19)
subjectto δ ∈U ,
1 ad
where u(δ ) is given as the solution of E(u(δ ),δ ) = 0. The existence of the function u is obtained
1 1 1
by the implicit function theorem. According to the ideas exposed in [21, 22], this can be done since
U = [0,L] is a nonempty, closed and convex set, J and E are continuously Fr´echet-differentiable
ad
functions,andassumingthatforeachδ ∈U thereexistsauniquecorrespondingsolutionu(δ )such
1 ad 1
that E(u(δ ),δ ) = 0 and the derivative ∂E(u(δ ),δ ) is a continuous linear operator continuously
1 1 ∂u 1 1
invertible for all δ ∈U .
1 ad
In order to find a minimum of the continuosly differentiable function J˜, it will be important to
compute the derivative of this reduced objective function. Hence, we will show a procedure to obtain
J˜′ by using the adjoint approach. Since
J˜′(δ )= u′(δ ) ∗∂J(u(δ ),δ )+ ∂J (u(δ ),δ ). (20)
1 1 1 1 1 1
∂u ∂δ
(cid:0) (cid:1) 1
Let us consider ζ ∈V ×Z as the solution of the so-called adjoint problem:
∂J ∂E ∗
(u(δ ),δ )+ (u(δ ),δ ) ζ =0. (21)
1 1 1 1
∂u (cid:18)∂u (cid:19)
where ∂E(u,δ ) ∗ is the adjoint operator of ∂E(u,δ ). Note that each term in (21) is an element of
∂u 1 ∂u 1
the spa(cid:0)ce V∗. (cid:1)
An equation for the derivative u′(δ ) is obtained by differentiating the equation E(u(δ ),δ ) = 0
1 1 1
with respect to δ :
1
∂E ∂E
(u(δ ),δ )u′(δ )+ (u(δ ),δ )=0, (22)
1 1 1 1 1
∂u ∂δ
1
where 0 is the zero vector in V∗×Z∗.
By using (20) we have that:
J˜′(δ ) = u′(δ ) ∗∂J(u(δ ),δ )+ ∂J (u(δ ),δ )
1 1 1 1 1 1
∂u ∂δ
(cid:0) (cid:1) 1
∗ −∗
∂E ∂E ∂J ∂J
= − (u(δ ),δ ) (u(δ ),δ ) (u(δ ),δ )+ (u(δ ),δ )
1 1 1 1 1 1 1 1
(cid:18)∂δ (cid:19) (cid:18)∂u (cid:19) ∂u ∂δ
1 1
∗
∂E ∂J
= (u(δ ),δ ) ζ+ (u(δ ),δ ),
1 1 1 1
(cid:18)∂δ (cid:19) ∂δ
1 1
7
where in the second equation we used (22) and for the last equation we used (21). Then:
∗
∂J ∂E
J˜′(δ )= (u(δ ),δ )+ (u(δ ),δ ) ζ. (23)
1 1 1 1 1
∂δ (cid:18)∂δ (cid:19)
1 1
NoticethatinordertoobtainJ˜′(δ )weneedfirsttocomputeu(δ )bysolvingthedirectproblem,
1 1
followed by the calculation of ζ by solving the adjoint problem. For computing the second term of
(23) it is not necessary to obtain the adjoint of ∂E(u(δ ),δ ) but just its action over ζ.
∂δ1 1 1
6. Getting the derivative of the functional.
In order to obtain the adjoint operator of ∂E, we have to find ∂E ∗ such that:
∂u ∂u
(cid:0) (cid:1)
∗
∂E ∂E
η,ζ = η, ζ , (24)
(cid:28)∂u (cid:29) (cid:28) (cid:18)∂u(cid:19) (cid:29)
where η = [η ,η ,η ]T is the direction of descent for the state variables u , u and u , respectively,
1 2 3 1 2 3
then
∂E hE(u+µη,δ ),ζi−hE(u,δ ),ζi
1 1
(u,δ )η,ζ = lim .
1
(cid:28)∂u (cid:29) µ→0+ µ
After some algebraics, it can be shown that ∂E (u,δ )η is given by:
∂u 1
∂E T 1 ∂η
1
(u,δ )η,ζ = −η (1−2u )+δ η u +δ u η λ dxdt+
1 1 1 1 1 3 1 1 3 1
(cid:28)∂u (cid:29) Z Z (cid:18) ∂t (cid:19)
0 0
T 1 ∂η
2
−ρ η (1−2u ) λ dxdt+
2 2 2 2
Z Z (cid:18) ∂t (cid:19)
0 0
T 1 ∂u ∂η ∂λ
2 2 2
−D η +D (1−u ) dxdt+
2 1 2 1
Z Z (cid:18) ∂x ∂x (cid:19) ∂x
0 0
T 1 ∂η T 1 ∂η ∂λ
3 3 3
−δ (η −η ) λ dxdt+ dxdt+
3 2 3 3
Z Z (cid:18) ∂t (cid:19) Z Z ∂x ∂x
0 0 0 0
1 1 1
η (x,0)γ (x)dx+ η (x,0)γ (x)dx+ η (x,0)γ (x)dx. (25)
1 1 2 2 3 3
Z Z Z
0 0 0
An inspection over equations (24) and (25) shows that, roughly speaking, we should remove the
spatial and temporal derivatives from η and pass them to λ.
The calculations make use of successive integration by parts to express each derivative of η in
terms of a derivative of λ. Omitting here the details, that are shown in the Appendix, we obtain the
following expression of the adjoint problem (21), which consists in finding λ∈V satisfying
T 1 ∂λ ∂u ∂λ
1 2 2
0 = − η −η (1−2u )λ +δ η u λ −D η dxdt+
1 1 1 1 1 1 3 1 2 1
Z Z (cid:18) ∂t ∂x ∂x (cid:19)
0 0
T 1 ∂λ ∂λ ∂η
2 2 2
− η −ρ η (1−2u )λ +D (1−u ) −δ η λ dxdt+
2 2 2 2 2 2 1 3 2 3
Z Z (cid:18) ∂t ∂x ∂x (cid:19)
0 0
T 1 ∂λ ∂λ ∂η T 1
3 3 3
− η +δ η λ + +δ u η λ dxdt+ η (u −uˆ )dxdt
3 3 3 3 1 1 3 1 3 3 3
Z Z (cid:18) ∂t ∂x ∂x (cid:19) Z Z
0 0 0 0
∂λ
= − ,η +hH(λ),ηi , (26)
(cid:28) ∂t (cid:29) V∗,V
V∗,V
for all η ∈V and λ(x,T)=0. As we show in the Appendix we can define γ(x)=λ(x,0).
Equation (26) shall be solved in order to get λ. Notice that the adjoint equations are posed
backwards in time, with a final condition at t = T, while the state equations are posed forward in
time, with an initial condition at t=0.
8
Inordertoobtainthe derivativeofthe functional,accordingto (23), wemustcomputethe deriva-
tive of E with respect to δ . Since
1
∂E hE(u,δ +µq),ζi−hE(u,δ ),ζi
1 1
(u,δ )q,ζ = lim ,
1
(cid:28)∂δ1 (cid:29) µ→0+ µ
for q ∈U , then
ad
∂E T 1
(u,δ )q,ζ =q u u λ dxdt.
1 1 3 1
(cid:28)∂δ1 (cid:29) Z0 Z0
Thus, since ∂J =0, we obtain an expression for (23), that is
∂δ1
∂E ∗ T 1
J˜′(δ )= (u(δ ),δ ) ζ = u u λ dxdt. (27)
1 1 1 1 3 1
(cid:18)∂δ1 (cid:19) Z0 Z0
7. Designing an algorithm to solve the minimization problem.
Itis worthstressingthat obtainingmodel parametersvia minimization ofthe objectivefunctional
J˜ is in general an iterative process requiring the value of the derivative. To compute J˜′ we just
solve two weak PDEs problems per iteration: the direct and the adjoint problems. This method is
muchcheaperthanthe sensitivity approach[22]inwhichthe directproblemis solvedmany times per
iteration. We developanimplementation in MATLAB that solvesthe direct andadjointproblems by
usingaFiniteElementMethodandtheoptimizationproblemissolvedbyusingaSequentialQuadratic
Programming (SQP) method , using the built-in function fmincon. For the direct problem, Figure 3
shows the density of health cells, tumor cells and excess concetration of H+ at fixed time (t=20) in
terms of x variable.
1 1
0.9 0.9
0.8 0.8
0.7 0.7
sity0.6 uu1 sity0.6 uu1
n0.5 2 n0.5 2
e u e u
D 3 D 3
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x x
Figure3: Densityofhealth,tumorcellsandexcessconcetrationofH+ atfixedtime(t=20)intermsofxvariable,for
δ1=0.5(left)andδ1=12.5(right).
It is well-known [23] that gradient-based optimization algorithms require the evaluation of the
gradient of the functional. One important advantage of evaluating the gradient through adjoints is
thatitrequirestosolvetheadjointproblemonlyonceperiteration,regardlessthenumberofinversion
variables. Note that the derivative of the functional can be approximated by using Finite Element
Method.
The method we will use for minimizing the functional J˜can be summarized as follows:
Algorithm 7.1. Adjoint-based minimization method.
1. Give an initial guess δ0 for the parameter.
1
2. Given δk in step k, solve the direct and adjoint problems at this step.
1
9
3. Obtain the derivative of the functional, i.e. J˜′(δk), using (27).
1
4. Obtain δk+1 by performing one iteration of the SQP method.
1
5. Stop using the criteria of fmincon.
Algorithm 7.2. Direct problem.
∂u u(·,t )−u(·,t )
n n−1
1. Do an implicit Euler step to find the state variables u: (·,t ) ≈ =
n
∂t τ
F(u(·,t )), where t =t +τ, F(u(·,t )) is a nonlinear functional and the intial condition is
n n n−1 n
u0(x)=u(x,0).
nod
2. Use FEM to make a discretization of u (x,t ) ≈ un φ (x), i = 1,2,3, φ (x) are the linear
i n i,j j j
jP=1
shape function and we note Un = [un ,···,un ,···,un ] ∈ Rnod, Un = [Un,Un,Un] ∈ Rq,
i i,1 i,j i,nod 1 2 3
where nod is the number of uniform distributed nodes for the spatial meshgrid for [0,1].
3. Use the Newton method to solve: find Un ∈ Rq such as Un−Un−1−τG(Un) = 0, where G is
the discretization of F.
Algorithm 7.3. Adjoint problem.
∂λ λ(·,t )−λ(·,t )
n n−1
1. Do an implicit Euler step to find the adjoint variable λ: − (·,t ) ≈ − =
n
∂t τ
H(λ(·,t )), and the final condition is λ(·,T)=0.
n−1
2. UseFEMtomakeadiscretizationofλ(·,t )andsolvethelinearproblemλn−1−λn−τK(λn−1)=
n
0, where K is the discretization of H.
8. Numerical experiments.
The goal of this section is to test and evaluate the performance of an adjoint-based optimization
method, by executing some numerical simulations of Algorithm 7.1 for some test-cases.
First consider an optimization problem that consist in minimizing the functional defined in (19),
where uˆ (x,t) is generated via the forward model, for a choice of the model parameters ρ = 1,
3 2
D = 4×10−5, δ = 1 and δˆ = 0.5, 4, 12.5, 16. We chose different values of δˆ because each one of
2 3 1 1
these shows a different behavior of tumor invasion, according to [7].
Figure4showsthevaluethatthefunctionaldefinedin(19)takesfordifferentvaluesofδ ,remaining
1
the other parameters constant. It is worth mentioning that J˜looks convex with respect to δ .
1
x 10−3
16
14
12
10
˜J
8
6
4
2
8 9 10 11 12 δ 13 14 15 16 17
1
Figure4: ThefunctionalJ˜foruˆ3 generated withδˆ1=12.5.
The idea of this test case is to investigate how close the original value of the parameter can be
retrieved. However, it is not a trivial one, because we do not know, for instance, if the optimization
10