Adjoint method for a tumor invasion PDE-constrained optimization problem in 2D using Adaptive Finite Element Method.
In this paper we present a method for estimating unknown parameter that appear in a two dimensional non-linear reaction-diffusion model of cancer invasion. This model considers that tumor-induced alteration of micro-environmental pH provides a mechanism for cancer invasion. A coupled system reaction-diffusion describing this model is given by three partial differential equations for the 2D non-dimensional 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 re-absorption 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 trust-region-reflective method including the computation of the derivative of the functional.
keywords:reaction-diffusion 2D equation, tumor invasion, PDE-constrained optimization, adjoint method, Adaptive Finite Element Method, Splitting Method, Trust-region-reflective method
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 acid-mediated invasion hypothesis. Therefore, it can be represented mathematically as a reaction-diffusion 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 tumor-host 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 large-scale morphological aspects of tumor necrosis such as central necrosis. Instead, they concentrate on the interactions of microscopic scale populations that occur at the tumor-host interface, arguing that these processes influence the clinically significant manifestations of invasive cancer.
Moreover, in Gatenby and Gawlinski (1996), the authors suppose 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 micro-environment 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 tumor-host 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 sub-cellular 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 tumor-host interface in invasive cancers are controlled by the same parameters which generate a transformation from a benign tumor into malignant tumor. A hypo-cellular 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 PDE-constrained optimization problem that can be solved with the splitting method, exploiting the fact that this procedure is easily parallelizable. We follow the PDE-based model by Gatenby Gatenby and Gawlinski (1996) in a two-dimensional 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 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 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 gradient-based 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 patient-specific 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échet-derivative and by and the partial Fréchet-derivatives 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 tumor-host 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 non-dimensional mathematical model is:
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 non-dimensional 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
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:
where the objective functional is
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 ill-posed 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 so-called reduced problem
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
where solves the so-called adjoint problem
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 Runge-Kutta method (ode45 MATLAB built-in 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 well-known Nocedal and Wright (2006) that gradient-based 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 built-in 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:
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 residual-type error estimator of with respect to the mesh is
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 built-in 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 .
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 built-in function was used to solve the minimization problem. The chosen algorithm in the fmincon function was the trust-region-reflective method, where the derivative of the objective function was computed according to 9.
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 adjoint-based 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) i7-3770K 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 speed-up. 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.
On the other hand, the method used for the minimization algorithm (3.3.1) is the trust-region-reflective 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 active-set or the interior-point 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.3.1 are: the feasible set for the optimization problem (5) is , the method used is trust-region-reflective, the option GradObj is on, and the maximum of function evaluations is 100.
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 well-known 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 well-known, 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 2-4, for each value of , we show the average for 10 values of , the standard deviation and the relative error .
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 two-dimensional 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 gradient-based 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 .
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, SECYT-UNC and PICT-FONCYT.
- 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 reaction-diffusion 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 size-structured 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, Three-dimensional 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 image-driven parameter estimation problem for a reaction-diffusion 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 PDE-constrained 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 PDE-constrained optimization problem using FEM, Serie A (2013).
- Ladyzhenskai︠a︡ and Solonnikov (1988) O. A. Ladyzhenskai︠a︡, V. A. Solonnikov, Linear and quasi-linear 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, Acid-mediated 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 posteriori-a 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.