Adjoint method for a tumor invasion PDEconstrained optimization problem in 2D using Adaptive Finite Element Method.
Abstract
In this paper we present a method for estimating unknown parameter that appear in a two dimensional nonlinear reactiondiffusion model of cancer invasion. This model considers that tumorinduced alteration of microenvironmental pH provides a mechanism for cancer invasion. A coupled system reactiondiffusion describing this model is given by three partial differential equations for the 2D nondimensional spatial distribution and temporal evolution of the density of normal tissue, the neoplastic tissue growth and 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 reabsorption rate and the destructive influence of H ions in the healthy tissue.
After solving the direct problem, we propose a model for the estimation of parameters by fitting the numerical solution with real data, obtained via in vitro experiments and fluorescence ratio imaging microscopy. We define an appropriate functional to compare both the real data and the numerical solution using the adjoint method for the minimization of this functional.
We apply a splitting strategy joint with Adaptive Finite Element Method (AFEM) to solve the direct problem and the adjoint problem. The minimization problem (the inverse problem) is solved by using a trustregionreflective method including the computation of the derivative of the functional.
keywords:
reactiondiffusion 2D equation, tumor invasion, PDEconstrained optimization, adjoint method, Adaptive Finite Element Method, Splitting Method, Trustregionreflective method1 Introduction
Cancer is one of the diseases causing the most deaths in the world, despite the best efforts of medicine. Human and financial resources are devoted for cancer research, and on several occasions these efforts are successful Adam (1986); Adam and Bellomo (1997); Bellomo et al. (2009, 2008); Byrne (2010); Bellomo and Preziosi (2000).
Some comments on the importance of mathematical modeling in cancer can be found in the literature. In the work Bellomo et al. (2008) the authors say “Cancer modelling has, over the years, grown immensely as one of the challenging topics involving applied mathematicians working with researchers active in the biological sciences. The motivation is not only scientific as in the industrial nations cancer has now moved from seventh to second place in the league table of fatal diseases, being surpassed only by cardiovascular diseases.”
We use the analysis proposed by Gatenby in Gatenby and Gawlinski (1996), which supports the acidmediated invasion hypothesis. Therefore, it can be represented mathematically as a reactiondiffusion system which describes the spatial and temporal evolution of the tumor tissue, normal tissue, and excess concentration of H.
The model simulates a pH gradient extending from the tumorhost interface. The effect of biological parameters that control this transition is supported by experimental and clinical observations Martin and Jain (1994).
Some authors Gatenby and Gawlinski (1996) model tumor invasion in order to find an underlying mechanism by which primary and metastatic cancers invade and destroy normal tissues. They do not attempt to model the genetic changes that lead to the transformation and seek to understand the causes of these changes. Likewise, they do not attempt to model the largescale morphological aspects of tumor necrosis such as central necrosis. Instead, they concentrate on the interactions of microscopic scale populations that occur at the tumorhost interface, arguing that these processes influence the clinically significant manifestations of invasive cancer.
Moreover, in Gatenby and Gawlinski (1996), the authors suppose that transformationinduced 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 the tumor cells survive and proliferate, while normal cells may not remain viable. The following temporal sequence would derive: (a) a high concentration of H ions in tumors will diffuse chemically as a gradient to adjacent normal tissue, exposing these normal cells to an interstitial pH like in the tumor, (b) normal cells immediately adjacent to the edge of the tumor are unable to survive in chronically this acid environment, and (c) progressive loss of normal cell layers in the tumorhost interface facilitates tumor invasion. Key elements of this mechanism of tumor invasion include low pH due to primitive metabolism and reduced viability of normal tissue in an acidic environment.
This model depends only on a small number of cellular and subcellular parameters. The analysis of the equations shows that the model simulates a crossover from a benign tumor to a malignant invasive tumor when some combination of parameters turn over some threshold value.
The structure and dynamics of the tumorhost interface in invasive cancers are controlled by the same parameters which generate a transformation from a benign tumor into malignant tumor. A hypocellular interstitial space, as we can see in Figure 1 (Gatenby and Gawlinski, 1996, Figure 4a), occurs in some cancers.
In this paper we propose a framework via a PDEconstrained optimization problem that can be solved with the splitting method, exploiting the fact that this procedure is easily parallelizable. We follow the PDEbased model by Gatenby Gatenby and Gawlinski (1996) in a twodimensional tissue. We estimate one of the model parameters (the destructive influence of H ions in the healthy tissue) using an inverse problem. It is possible to get data about the concentration of hydrogen ions Martin and Jain (1994) via fluorescence ratio imaging microscopy. In this approach, tumor invasion is modeled via a coupled nonlinear system of partial differential equations, which makes the numerical solution procedure quite challenging. These equations are solved using an Adaptive Finite Element Method (AFEM).
This kind of problem constitutes a particular application of the socalled inverse problems, which are being increasingly used in a broad number of fields in applied sciences. For instance, problems referred to structured population dynamics Perthame and Zubelli (2007), computerized tomography and image reconstruction in medical imaging van den Doel et al. (2011); Zubelli et al. (2003), and more specifically tumor growth Agnelli et al. (2011); Hogea et al. (2008); Knopoff et al. (2013), among many others.
We solve a minimization problem using a gradientbased method considering the adjoint method in order to find the derivative of an objective functional. In this way, we would obtain the best parameter that fits patientspecific data.
This work follows the ideas Quiroga et al. (2013) where the space variable was in a one dimensional space. The extension of the model to two dimensional space allows us to approach the results to more realistic biological hypotheses.
The contents of this paper, which is organized into 5 sections, are as follows: Section 2 consists in some preliminaries about the model, the definition of variational form of the direct and adjoint problems, and the minimization problem. Section 3 deals with suitable numerical algorithms to solve the direct and adjoint problems. In particular, we use the splitting method and the Adaptive Finite Element Method with a computation of a posteriori error. In Section 4 we show numerical simulations of the retrieved parameter and the need of a parallel scheme. Section 5 presents the conclusions and some future work related to the contents of this paper.
Some words about our notation. We use to denote the 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 such that , we denote by the full Fréchetderivative and by and the partial Fréchetderivatives of at . For a linear operator we denote the adjoint operator of . If is invertible, we call the inverse of the adjoint operator .
2 Some preliminaries about the model
A mathematical model of the tumorhost interface based on the acid mediation hypothesis of tumor invasion due to Gatenby and Gawlinski (1996) is given by the following system of partial differential equations (PDEs):
where the variables are in . These equations determine the spatial distribution and temporal evolution of three fields: , the density of normal tissue; , the density of neoplastic tissue; and , the excess concentration of H ions. The units of and are cells/cm and is expressed as a molarity (M). The space and time are given in cm and seconds, respectively.
The biological meaning of each equation can be seen in Gatenby and Gawlinski (1996), and the nondimensional mathematical model is:
(1)  
where the four dimensionless quantities which parameterize the model are given by: , , and .
The interaction parameters between different cells (healthy and tumor) and concentration of H are difficult to measure experimentally. This is the reason for which we propose to estimate them, so we will focus on in this work.
The initial and boundary conditions considered for the nondimensional system are:
From now on, equations (1) with the initial and boundary conditions will be referred to as the direct problem.
2.1 Variational form for the direct problem
Using the variational techniques for obtaining the weak solution of the direct problem Ladyzhenskai︠a︡ and Solonnikov (1988); Kinderlehrer and Stampacchia (1987); Evans (1998), it can be written as where such that
where , , , with
and . We use , with
(2) 
A weak solution is a function that satisfies for all .
2.2 Formulation of the minimization problem
As described above we propose to use an inverse problem technique in order to estimate . Function represents the solution of the direct problem (the components of are the state variables of the problem) for each choice of the parameter .
Let us assume that experimental information is available during the time interval . Then, the inverse mathematical problem can be formulated as:
(3) 
where the objective functional is
(4) 
with , the excess concentration of H ions obtained by solving the direct problem for a certain choice of and , the excess concentration measured experimentally (real data). The constraints are given by , a subset of , the set of admissible values of and is the weak formulation of the direct problem.
We remark that, in general, there is a fundamental difference between the direct and the inverse problems. In fact, the latter is usually illposed in the sense of existence, uniqueness and stability of the solution. This inconvenient is often treated by using some regularization techniques van den Doel et al. (2011); Engl et al. (1996); Kirsch (2011).
2.3 Formulation of the reduced and adjoint problems
In the following, we will consider the socalled reduced problem
(5) 
where is given as the solution of . In order to find a minimum of the continuously differentiable function , it will be important to compute the derivative of this reduced objective function. Hence, we will show a procedure to obtain by using the adjoint approach. According to the theory exposed in Brandenburg et al. (2009); Hinze (2009), the derivative of is given by
(6) 
where solves the socalled adjoint problem
(7) 
Notice that in order to obtain we need first to compute by solving the direct problem, followed by the calculation of by solving the adjoint problem. For computing the second term of (6) it is not necessary to obtain the adjoint of but just its action over .
3 Designing an algorithm to solve the minimization problem
It is worth stressing that obtaining model parameters via minimization of the objective functional is in general an iterative process requiring the value of the derivative. To compute we just solve two weak PDEs problems per iteration: the direct and the adjoint problems. This method is much cheaper than the sensitivity approach Hinze (2009) in which the direct problem is solved many times per iteration. We develop an implementation in MATLAB that solves the direct and adjoint problems. We use the splitting method in order to separate the direct problem in two new problems. The first one consists in a system of ordinary differential equations that contains the reaction terms of the original PDE. The second one is a PDE that contains the diffusion terms of the original PDE. The ODEs are solved by using the RungeKutta method (ode45 MATLAB builtin function). Since we have an ODE system for each spatial point, its resolution can be parallelized accelerating the time execution. The PDE is solved by the Adaptive Finite Element Method. In the next subsections we will explain the splitting method and the adaptive procedure of the FEM. It is wellknown Nocedal and Wright (2006) that gradientbased optimization algorithms require the evaluation of the gradient of the functional. The optimization problem is solved by using a Sequential Quadratic Programming (SQP) method , using the builtin function fmincon.
For the direct problem, Figure 3 shows the density of health cells, tumor cells and excess concentration of H at fixed time () in terms of variable.
3.1 Solving the direct problem
3.1.1 Splitting method
A multiscale operator splitting
We proceed like in (Estep et al., 2008, Section 2). For the time discretization, we introduce a theoretical framework in which each component (the reaction component and the diffusion component ) is solved exactly. We define a piecewise continuous approximate solution:
for , with the nodal values obtained from the following procedure. We first discretize into with diffusion time step , for . For each diffusion step, we choose a (small) time step where , with , and the nodes (see figure 4). We associate the time intervals and with these discretizations.
3.1.2 Adaptive FEM
The adaptive procedure for FEM consists in a four step loop: (a) solve the PDE using the FEM discretization, (b) estimate a posteriori error of the discrete solution, (c) mark the elements to be refined according to the relative error size of the a posteriori error, and (d) refine the marked elements keeping the mesh conformity (for more details see Nochetto et al. (2009)).
Given a mesh at time , the element residual and the jump residual are defined as:
(10)  
(11) 
were are the edges of and is the strong form of the operator as defined in (2).
We define the local error indicator by
were is the diameter of and is the length of the edge . If is an edge of an element, then
The residualtype error estimator of with respect to the mesh is
3.1.3 Algorithm

STEP 0: Set an initial condition on the coarse uniform mesh . Set .

STEP 1: Given do the following steps to compute if .

STEP 2: Compute satisfying the reaction equation:
for and for all .

STEP 3: Compute satisfying the diffusion equation:
for and for all . Set .

STEP 4: Compute the a posteriori error . If , set and go to STEP 1.

STEP 5: Mark and refine, and go to STEP 2.
In STEP 2, we compute , the reaction component, by using the ode45 MATLAB builtin function for each node of the current mesh, allowing a parallelization strategy. In STEP 3, the diffusion component is solved by using FEM. In STEP 5, according to Byfut et al. (2007), we use the bulk algorithm to mark and the RedGreenBlue algorithm to refine. The bulk algorithm defines the set of marked edges such that
or it contains all the edges of marked elements that satisfy
where is the set of edges of and .
3.2 Solving the adjoint problem
In order to solve the adjoint problem we shall use FEM. The spatial discretization is the coarse mesh used for the initial mesh in the direct problem. We denote for .
3.2.1 Algorithm

STEP 0: Set the final condition on the initial mesh .

STEP 1: Given do the following steps to compute if .

STEP 2: Do an implicit Euler step in time, and FEM in space to approximate the adjoint variable by solving the linear system , where is the discretization of as defined in (8).

STEP 3: Set and go to STEP 1.
3.3 Solving the minimization problem
The fmincon MATLAB builtin function was used to solve the minimization problem. The chosen algorithm in the fmincon function was the trustregionreflective method, where the derivative of the objective function was computed according to 9.
3.3.1 Algorithm
The method we will use for minimizing the functional can be summarized as follows:

STEP 0: Give an initial guess for the parameter.
In order to compute and is necessary to solve the direct and adjoint problems.
4 Numerical experiments
The goal of this section is to test and evaluate the performance of an adjointbased optimization method, by executing some numerical simulations of Algorithm 3.3.1 for some test cases. The experiments were run in MATLAB, in a PC running Linux, Intel(R) Core(TM) i73770K CPU, 3.50GHz.
First consider an optimization problem that consists in minimizing the functional defined in (5), where is generated via the direct problem for some with the choice of model parameters , and . We choose several values of , for instance , because they show a different behavior of tumor invasion, according to Gatenby and Gawlinski (1996).
Figure 5 shows the value that the functional defined in (5) takes for different values of , for generated with . It is worth mentioning that, even when we do not know in advance if the optimization problem has a unique solution, looks convex with respect to .
The idea of these test cases is to investigate how close the original value of the parameter can be retrieved (even in the presence of noise), and how efficiently these computations can be done. Regarding computational efficiency, one of the most expensive parts of Algorithm 3.3.1 is the resolution of a system of ODEs (STEP 2 of Algorithm 3.1.3). Since we have a system of ODEs for each node of a current mesh, a parallel strategy (each processor solves a system of ODES for one node) is the best and natural option to reduce time execution. For this particular inverse problem, parallelization is not an option, but it is a need, because the direct problem could be called many times by the optimization solver. For example, figure 6(a) shows how many seconds takes to solve the direct problem. Figure 6(b) shows the speedup. In addition, since we have to go through all the nodes in order to compute the a posteriori error, we have also parallelized this computation.
(a)  (b) 
On the other hand, the method used for the minimization algorithm (3.3.1) is the trustregionreflective method Coleman and Li (1996); matlab , where the option GradObj is on by default (the gradient of the objective function must be supplied). If we would use another algorithm (like the activeset or the interiorpoint algorithms) where the gradient can be estimated by finite differences, then the computational cost of this algorithm makes it not practical. That is the reason for which we compute the exact derivative of the functional using the adjoint method.
We have run the Algorithm 3.3.1 for different values of taking the initial condition randomly. Averaging the different solutions and taking the standard deviation of all of these experiments, the retrieved parameter is obtained very accurately.
The algorithmic parameters for Algorithm 3.1.3 and 3.2.1 are: , , the initial coarse mesh has triangular elements, , and .
The algorithmic parameters for Algorithm 3.3.1 are: the feasible set for the optimization problem (5) is , the method used is trustregionreflective, the option GradObj is on, and the maximum of function evaluations is 100.
4  4.2666  7.0640 

12.5  12.4937  7.1875 
16  16.6246  1.8826 
We emphasize that we have retrieved accurately the value of independently of the value of . Thus, in the next experiment we will consider a fixed value .
It is wellknown that the presence of noise in the data may imply the appearance of strong numerical instabilities in the solution of an inverse problem Bertero and Piana (2006).
One of the experimental methods to obtain values of is by using fluorescence ratio imaging microscopy Martin and Jain (1994).
As it is wellknown, measurements are often affected by perturbations, usually random ones. Then we perform numerical experiments where is perturbed by using Gaussian random noise with zero mean and standard deviation . In Tables 24, for each value of , we show the average for 10 values of , the standard deviation and the relative error .
0.1000  3.1125  0.8624  0.2219 
0.1500  3.5409  1.8611  0.1148 
0.2000  3.4471  2.3701  0.1382 
0.1000  11.6235  2.6314  0.0701 
0.1500  12.3825  4.6561  0.0094 
0.2000  11.9537  5.5749  0.0437 
0.1000  16.6996  2.1280  0.0437 
0.1500  17.0926  5.1026  0.0683 
0.2000  17.6308  2.4753  0.1185 
5 Final conclusions and future work
In this paper we have solved a parameter estimation problem following the model proposed by Gatenby and Gawlinski (1996) in a twodimensional space. The inverse problem is formulated as an optimization problem in order to find the parameter (the destructive influence of H ions in the healthy tissue).
The direct problem was solved by the splitting technique together with Adaptive Finite Element Method, for the purpose of controlling the numerical error and defining a parallel strategy. A gradientbased method was used to solve the optimization problem. The derivative of the objective functional was computed using the solution of the adjoint problem.
The experiments were run in MATLAB recovering several values of the parameter representing different scenarios. Also, a stability analysis was performed using random noise to simulate perturbations in the experimental data.
We consider that the results are accurately enough. In most cases the parameters are retrieved with a relative error less than .
Acknowledgments
We appreciate the courtesy of Sebastian Pauletti, from IMAL, Santa Fe, Argentina, who strongly contributed with information above splitting method.
The work of the authors was partially supported by grants from CONICET, SECYTUNC and PICTFONCYT.
References
 Adam (1986) J. Adam, A simplified mathematical model of tumor growth, Mathematical Biosciences 81 (1986) 229–244.
 Adam and Bellomo (1997) J. Adam, N. Bellomo, A survey of models for tumor immune systems dynamics, Modeling and simulation in science, engineering & technology, Birkhäuser, 1997.
 Bellomo et al. (2009) N. Bellomo, M. Chaplain, E. De Angelis, Selected Topics on Cancer Modeling  Genesis  Evolution  Immune Competition  Therapy, Birkhäuser, Boston, 2009.
 Bellomo et al. (2008) N. Bellomo, N. Li, P. Maini, On the foundations of cancer modelling: selected topics, speculations, and perspectives, Mathematical Models and Methods in Applied Sciences 18 (2008) 593–646.
 Byrne (2010) H. M. Byrne, Dissecting cancer through mathematics: from the cell to the animal model, Nature Reviews Cancer 10 (2010) 221–230.
 Bellomo and Preziosi (2000) N. Bellomo, L. Preziosi, Modelling and mathematical problems related to tumor evolution and its interaction with the immune system, Mathematical and Computer Modelling 32 (2000) 413–452.
 Gatenby and Gawlinski (1996) R. A. Gatenby, E. T. Gawlinski, A reactiondiffusion model of cancer invasion, Cancer Research 56 (1996) 5745–5753.
 Martin and Jain (1994) G. R. Martin, R. K. Jain, Noninvasive measurement of interstitial pH profiles in normal and neoplastic tissue using fluorescence ratio imaging microscopy, Cancer Research 54 (1994) 5670–5674.
 Perthame and Zubelli (2007) B. Perthame, J. Zubelli, On the inverse problem for a sizestructured population model, Inverse Problems 23 (2007) 1037–1052.
 van den Doel et al. (2011) K. van den Doel, U. M. Ascher, D. K. Pai, Source localization in electromyography using the inverse potential problem, Inverse Problems 27 (2011) 025008.
 Zubelli et al. (2003) J. Zubelli, R. Marabini, C. Sorzano, G. Herman, Threedimensional reconstruction by Chahine’s method from electron microscopic projections corrupted by instrumental aberrations, Inverse Problems 19 (2003) 933–949.
 Agnelli et al. (2011) J. Agnelli, A. Barrea, C. Turner, Tumor location and parameter estimation by thermography, Mathematical and Computer Modelling 53 (2011) 1527–1534.
 Hogea et al. (2008) C. Hogea, C. Davatzikos, G. Biros, An imagedriven parameter estimation problem for a reactiondiffusion glioma growth model with mass effects, J. Math. Biol. 56 (2008) 793–825.
 Knopoff et al. (2013) D. Knopoff, D. Fernández, G. Torres, C. Turner, Adjoint method for a tumour growth PDEconstrained optimization problem, Computers and Mathematics with Applications (accepted 2013).
 Quiroga et al. (2013) A. A. I. Quiroga, D. R. Fernández, G. A. Torres, C. V. Turner, Adjoint method for a tumor invasion PDEconstrained optimization problem using FEM, Serie A (2013).
 Ladyzhenskai︠a︡ and Solonnikov (1988) O. A. Ladyzhenskai︠a︡, V. A. Solonnikov, Linear and quasilinear equations of parabolic type, volume 23, American Mathematical Soc., 1988.
 Kinderlehrer and Stampacchia (1987) D. Kinderlehrer, G. Stampacchia, An introduction to variational inequalities and their applications, Society for Industrial and Applied Mathematics, 1987.
 Evans (1998) L. C. Evans, Partial differential equations (1998).
 Gatenby et al. (2006) R. A. Gatenby, E. T. Gawlinski, A. F. Gmitro, B. Kaylor, R. J. Gillies, Acidmediated tumor invasion: a multidisciplinary study, Cancer Research 66 (2006) 5216–5223.
 Engl et al. (1996) H. W. Engl, M. Hanke, A. Neubauer, Regularization of inverse problems, volume 375, Kluwer Academic Pub, 1996.
 Kirsch (2011) A. Kirsch, An introduction to the mathematical theory of inverse problems, volume 120, Springer Science+ Business Media, 2011.
 Brandenburg et al. (2009) C. Brandenburg, F. Lindemann, M. Ulbrich, S. Ulbrich, A continuous adjoint approach to shape optimization for Navier Stokes flow, in: Optimal control of coupled systems of partial differential equations, Springer, 2009, pp. 35–56.
 Hinze (2009) M. Hinze, Optimization with PDE constraints, volume 23, Springer, 2009.
 Nocedal and Wright (2006) J. Nocedal, S. J. Wright, Numerical optimization, Springer Science+ Business Media, 2006.
 Estep et al. (2008) D. Estep, V. Ginting, D. Ropp, J. N. Shadid, S. Tavener, An a posterioria priori analysis of multiscale operator splitting, SIAM Journal on Numerical Analysis 46 (2008) 1116–1146.
 Nochetto et al. (2009) R. H. Nochetto, K. G. Siebert, A. Veeser, Theory of adaptive finite element methods: an introduction, in: Multiscale, nonlinear and adaptive approximation, Springer, 2009, pp. 409–542.
 Byfut et al. (2007) A. Byfut, J. Gedicke, D. Günther, J. Reininghaus, S. Wiedemann, et al., FFW documentation, Humboldt University of Berlin, Germany (2007).
 Coleman and Li (1996) T. F. Coleman, Y. Li, A reflective newton method for minimizing a quadratic function subject to bounds on some of the variables, SIAM Journal on Optimization 6 (1996) 1040–1058.
 (29) MATLAB, URL: http://www.mathworks.com/products/matlab.
 Bertero and Piana (2006) M. Bertero, M. Piana, Inverse problems in biomedical imaging: modeling and methods of solution, in: Complex systems in biomedicine, Springer, 2006, pp. 1–33.
 McGILLEN et al. (2012) J. B. McGILLEN, N. K. Martin, I. F. Robey, E. A. Gaffney, P. K. Maini, Applications of mathematical analysis to tumour acidity modelling, RIMS Kokyuroku Bessatsu 31 (2012) 31–59.
 Martin et al. (2012) N. Martin, I. Robey, E. Gaffney, R. Gillies, R. Gatenby, P. Maini, Predicting the safety and efficacy of buffer therapy to raise tumour phe: an integrative modelling study, British journal of cancer 106 (2012) 1280–1287.