Adjoint method for a tumor invasion PDE-constrained optimization problem using FEM
In this paper we present a method for estimating unknown parameter that appear on a non-linear reaction-diffusion model of cancer invasion. This model considers that tumor-induced alteration of micro-enviromental pH provides a mechanism for cancer invasion. A coupled system reaction-diffusion describing this model is given by three partial differential equations for the 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 reabsorption rate and the destructive influence of H ions in the healthy tissue.
After solving the forward problem properly, we use the 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 Finite Element Method (FEM) to solve both the direct and inverse problem, computing the a posteriori error.
keywords:reaction-diffusion equation, tumor invasion, PDE-constrained optimization, adjoint method, Finite Element Method, a posteriori error
Cancer is one of the greatest killers in the world although medical activity has been successful, despite great difficulties, at least for some pathologies. A great effort of human and economical resources is devoted, with successful outputs, to cancer research, 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 in this work the mathematical analyses first proposed by Gatenby and Gawlinski (1996) 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.
The model predicts a pH gradient extending from the tumor-host interface. The effect of biological parameters critical to controlling this transition is supported by experimental and clinical observations Martin and Jain (1994).
In Gatenby and Gawlinski (1996) 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 the genetic changes which result in transformation nor do they seek to understand the causes of these 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 cells survive and proliferate, whereas normal cells are unable to remain viable. The following temporal 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 crossover from 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 controlled by the same biological parameters which generate the transformation from benign to malignant growth. A hypocellular interstitial gap, as we can see in Figure 1 (Gatenby and Gawlinski, 1996, Figure 4a), at the interface is predicted to occur in some cancers.
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 Martin and Jain (1994). We propose a framework via a PDE-constrained optimization problem, following the PDE-based model by Gatenby Gatenby and Gawlinski (1996). In this approach, tumor invasion is modeled via a coupled nonlinear system of partial differential equations, which makes 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 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.
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. Section 3 deals with the variational formulation of the direct problem. Section 4 considers the formulation of the minimization problem. Section 5 introduces the reduced and adjoint problem, deriving the optimality conditions for the problem. Section 6 finds the derivative of the solution of a functional with respect to a parameter that does not appear explicitly in the equation. Section 7 deals with the numerical solution of the adjoint problem, designing a suitable algorithm to solve it. In particular, we use the Finite Element Method with a computation of a posteriori error. In Section 8 we show some numerical simulations to give information on the behavior of the functional and its dependence on the 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 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.
We present the mathematical model of the tumor-host interface based on the acid mediation hypothesis of tumor invasion due to Gatenby and Gawlinski (1996). For convenience we reproduce the equations here, which 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 excess H ion concentration is expressed as a molarity (M), and are the position (in cm) and time (in seconds), respectively.
where the variables are in .
In equation (1) the behavior of the normal tissue is determined by the logistic growth of with growth rate and carrying capacity , and the interaction of with excess H ions leading to a death rate proportional to . The number is the excess acid concentration, dependent death rate in accord with the well-described decline in the growth rate of normal cells, due to the reduction of pH from its optimal value of . The constants , and have units of s, l/(M s) and cells/cm, 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 with growth rate and carrying capacity . The diffusion term depends on the absence of healthy tissue with a diffusion constant . Constants , and have units of s, cells/cm and cm/s, respectively.
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-scale vascular evacuation Gatenby and Gawlinski (1996)). Constant is the production rate (M cm/(cell s)), is the reabsorption rate (1/s), and is the H ion diffusion constant (cm/s).
All the parameter values can be found in Table 1.
Following the ideas exposed in Gatenby and Gawlinski (1996), and considering that , the mathematical model is rescaled and the domain is transformed onto the interval . Hence, let us define the following functions:
for , where the four dimensionless quantities which parameterize the model are given by:
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 other parameters can be estimated by different techniques (see Table 1).
2.2 Initial and boundary conditions.
At we will consider the tumor at a certain stage of its evolution. Hence the initial conditions are:
for all . 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 , we have
3 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), we define the following weak formulation:
A weak solution is a function that satisfies (15) for all and .
4 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 problem can be formulated as:
Find a parameter able to generate data that best match the available (experimental) information over time .
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 .
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 Martin and Jain (1994); Gatenby et al. (2006) at certain times , . For example, Figure 2 (Gatenby et al., 2006, Figure 4) shows 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 could be defined as:
where is the excess concentration of H ions obtained by solving the direct problem for a certain choice of and is the excess concentration measured experimentally (real data).
Let us define such that
where , and .
In this way we can rewrite the weak formulation (15) as .
The parameter that best matches the experimental information with the generated data provided by the direct problem can be computed by solving a PDE-constrained optimization problem, namely:
where denotes the set of admissible values of . In our case, should be a subset of . Notice that a solution must satisfy the constraint , which constitutes 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).
5 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 . The existence of the function is obtained by the implicit function theorem. According to the ideas exposed in Brandenburg et al. (2009); Hinze (2009), this can be done since is a nonempty, closed and convex set, and are continuously Fréchet-differentiable functions, and assuming that for each there exists a unique corresponding solution such that and the derivative is a continuous linear operator continuously invertible for all .
In order to find a minimum of the continuosly 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. Since
Let us consider as the solution of the so-called adjoint problem:
where is the adjoint operator of . Note that each term in (21) is an element of the space .
An equation for the derivative is obtained by differentiating the equation with respect to :
where is the zero vector in .
By using (20) we have that:
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 (23) it is not necessary to obtain the adjoint of but just its action over .
6 Getting the derivative of the functional.
In order to obtain the adjoint operator of , we have to find such that:
where is the direction of descent for the state variables , and , respectively, then
After some algebraics, it can be shown that is given by:
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 satisfying
for all and . As we show in the Appendix we can define .
Equation (26) shall be solved in order to get . Notice that the adjoint equations are posed backwards in time, with a final condition at , while the state equations are posed forward in time, with an initial condition at .
7 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 by using a Finite Element Method and 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 concetration of H at fixed time () in terms of variable.
It is well-known Nocedal and Wright (2006) that gradient-based optimization algorithms require the evaluation of the gradient of the functional. One important advantage of evaluating the gradient through adjoints is that it requires to solve the adjoint problem only once per iteration, regardless the number of inversion 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 can be summarized as follows:
Adjoint-based minimization method.
Give an initial guess for the parameter.
Given in step , solve the direct and adjoint problems at this step.
Obtain the derivative of the functional, i.e. , using (27).
Obtain by performing one iteration of the SQP method.
Stop using the criteria of fmincon.
Do an implicit Euler step to find the state variables : , where , is a nonlinear functional and the intial condition is .
Use FEM to make a discretization of , , are the linear shape function and we note , , where is the number of uniform distributed nodes for the spatial meshgrid for .
Use the Newton method to solve: find such as , where is the discretization of .
Do an implicit Euler step to find the adjoint variable : , and the final condition is .
Use FEM to make a discretization of and solve the linear problem , where is the discretization of .
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 is generated via the forward model, for a choice of the model parameters , , and . We chose different values of because each one of these shows a different behavior of tumor invasion, according to Gatenby and Gawlinski (1996).
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 problem has a solution or, in that case, if it is unique or if the method converges to another local minima.
We have run the Algorithm 7.1 for different values of taking the initial condition randomly, as we can see from the Table 2 the retrieved parameter is obtained very accurately since the standard deviation is small. For Algorithms 7.2 and 7.3 we use the following algorithmic parameters and , and .
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 method to obtain values of is by using fluorescence ratio imaging microscopy Martin and Jain (1994). As it is well-known that measurements are often affected by perturbations, usually random ones.
Then we perform numerical experiments where is perturbed by using Gaussian random noise whith mean zero and standard deviation . In the next tables 3-6, for each value of , we show the average of 30 values of , the standard deviation and the relative error .
9 Final conclusions and future work.
A miscellany of new strategies, experimental techniques and theoretical approaches are emerging in the ongoing battle against cancer. Nevertheless, as new, ground-breaking discoveries relating to many and diverse areas of cancer research are made, scientists often have recourse to mathematical modelling in order to elucidate and interpret these experimental findings, Adam and Bellomo (1997); Bellomo et al. (2008); Byrne (2010); Araujo and McElwain (2004), and it became clear that these models are expected to success if the parameters involved in the modeling process are known. Or eventually, taking into account that some biological parameters may be unknown (especially in vivo), the model can be used to obtain them Agnelli et al. (2011); van den Doel et al. (2011).
This paper, as already mentioned in Section 1, aims at offering a mathematical tool for the obtention of phenomenological parameters which can be identified by inverse estimation, by making suitable comparisons with experimental data. The inverse problem was stated as a PDE-constrained optimization problem, which was solved by using the adjoint method. In addition, the gradient of the proposed functional is obtained and can be extended, in principle, to any number of unknown parameters.
We remark that the parameter estimation via PDE-constrained optimization is a general approach that can be used, for instance, to consider the effects of nonlinear interaction between the health and tumor cells McGillen et al. (2013).
As a future work we are interested in the dependence of the on time and in the dependence of the diffusivity coefficient of excess of the H concentration with respect to the space variable , as in Martin et al. (2010). Also we propose to solve the problem in two dimensional space, where the importance of using adaptive FEM will be crucial.
We appreciate the courtesy of Claudio Padra, from the Grupo de Mecánica Computacional - CNEA Bariloche - Argentina, who strongly contributed with information above FEM and a posteriori error.
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).
- 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, American Mathematical Society, 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.
- 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.
- Verfürth (????) R. Verfürth, A posteriori error estimates for nonlinear problems. finite element discretizations of elliptic equations, in: 475 (1994) MR 94j:65136, p. 445.
- Babuška and Rheinboldt (1978) I. Babuška, W. C. Rheinboldt, A-posteriori error estimates for the finite element method, International Journal for Numerical Methods in Engineering 12 (1978) 1597–1615.
- Araujo and McElwain (2004) R. Araujo, D. McElwain, A history of the study of solid tumour growth: the contribution of mathematical modelling, Bulletin of mathematical biology 66 (2004) 1039–1091.
- McGillen et al. (2013) J. B. McGillen, E. A. Gaffney, N. K. Martin, P. K. Maini, A general reaction–diffusion model of acidity in cancer invasion, Journal of mathematical biology (2013) 1–26.
- Martin et al. (2010) N. K. Martin, E. A. Gaffney, R. A. Gatenby, P. K. Maini, Tumour–stromal interactions in acid-mediated invasion: a mathematical model, Journal of theoretical biology 267 (2010) 461–470.
Appendix A Appendix: obtaining the adjoint problem.
In this section we show the calculations involved in order to obtain the adjoint equations (26). As stated in Section 6, the adjoint equations constitute a weak formulation of the adjoint problem, with unknown , given by (21). Here, is obtained by using (24). In what follows, we shall obtain equivalent expressions for each of the six terms of the summation , which are associated with the six constraints given by in (17).