ADJOINT SENSITIVITY ANALYSIS OF A TUMOR GROWTH MODEL AND ITS APPLICATION TO SPATIOTEMPORAL RADIOTHERAPY OPTIMIZATION

. We investigate a spatial model of growth of a tumor and its sensitivity to radiotherapy. It is assumed that the radiation dose may vary in time and space, like in intensity modulated radiotherapy (IMRT). The change of the ﬁnal state of the tumor depends on local diﬀerences in the radiation dose and varies with the time and the place of these local changes. This leads to the concept of a tumor’s spatiotemporal sensitivity to radiation, which is a function of time and space. We show how adjoint sensitivity analysis may be applied to calculate the spatiotemporal sensitivity of the ﬁnite diﬀerence scheme resulting from the partial diﬀerential equation describing the tumor growth. We demonstrate results of this approach to the tumor proliferation, invasion and response to radiotherapy (PIRT) model and we compare the ac-curacy and the computational eﬀort of the method to the simple forward ﬁnite diﬀerence sensitivity analysis. Furthermore, we use the spatiotemporal sensitivity during the gradient-based optimization of the spatiotemporal radiation protocol and present results for diﬀerent parameters of the model.


1.
Introduction. Intensity modulated radiotherapy (IMRT) is an advanced type of radiation therapy used to treat tumors. It is a type of conformal radiation, which shapes radiation beams to closely approximate the shape of the tumor. In general, it allows to deliver a different radiation dose to different points in 3D space. The goal of IMRT is to adjust the radiation dose to the target and to reduce exposure of healthy tissue to minimize the side effects of treatment. Conformal radiation may be used in one or more sessions. Radiotherapy can also be a component of more advanced combined anticancer therapies [14].
In current clinical practice IMRT is planned based on the current state (shape and volume) of the tumor and there are very few studies which try to use mathematical modeling to predict the response to the therapy. One of the first such attempts is an approach proposed by Corwin and co-authors in [1]. They used a spatial mathematical model of tumor growth to simulate and predict the response to a dose of particular shape, and based on this information they optimized the therapy using an evolutionary optimization algorithm. It is important to stress that in [1] the mathematical model was simulated only forward in time, which is, of course, a typical way of using any mathematical model.
The problem of optimizing therapy is related to the problem of sensitivity of an objective function, measuring the quality of the response, with respect to the therapy protocol.
In this work we use adjoint sensitivity analysis, which is a sensitivity method that utilizes the adjoint system constructed based on the analyzed model and simulated backward in time. Adjoint sensitivity analysis, in a special structural form, is especially useful for models described by block diagrams and was originally developed for neural networks [2] and afterwards was used for different models described by ordinary differential equations [3,4,8,9], systems with delays [5] and agestructured models [7]. Adjoint sensitivity analysis is especially useful in analyzing MISO (Multiple-Input Single-Output) systems due to the much lower computational complexity when compared to forward sensitivity analysis. For example, it is useful in a parameter estimation process where the system analyzed has many inputs (parameters) and usually one scalar output describing the system (objective function). In the present work adjoint sensitivity analysis is used to compute the spatiotemporal gradient of the objective function with respect to the input signal of the finite difference scheme being the result of the discretization of a given partial differential equation.

2.
Mathematical model of growth of a tumor. In this paper we derive the spatiotemporal sensitivity to radiation for a mathematical model for tumor proliferation, invasion and response to radiotherapy (PIRT model) published previously in [1,10,11]. The model is formalized as the following reaction-diffusion equation: where: D is a diffusion coefficient (tumor invasiveness), ρ means proliferation coefficient, k t is a environment capacity for tumor cells and R is a term describing the effect of radiation on the tumor cells. In this work R is formalized as a linearquadratic term: R = 1 − e −αd−βd 2 (2) where d is the radiation dose or, more precisely, spatiotemporal radiation dose distribution. Later in this article d is called shorter as the radiation dose distribution or simply as the dose distribution. In this work we use and optimize radiation as a continuous function of time, hence the term R has interpretation of rate constant with units 1/time. The interpretation of the term R has been a topic of a recent discussion -see articles [12,13] and the discussion therein. α and β are parameters of the linear-quadratic model.
We used the finite-difference method to discretize the model (1) which gives the following difference equation: The discretized model (3) is an example of a one dimensional finite difference scheme usually used in numerical simulations of spatiotemporal models. The simplest form of such a difference scheme is as follows: where: y j i is a discretized scalar variable of the original PDE at the time j and at the spatial point i, u j i is the value of an input signal affecting the system (if such a signal exists) at the time j and at the spatial point i, f (·) is a general non-linear scalar function, i max is the number of spatial points, and j max is the final discrete time.
The difference scheme (4) is very simple: it is only one dimensional, y and f (·) are scalar, and y j+1 i depends on the state of the difference scheme only at the previous j-th time and only in the closest neighborhood (i − 1, i, i + 1). In general, the difference scheme resulting from the PDE may be more complicated; it may be higher dimensional, y and f (·) may be multidimensional, and y j+1 i may depend on a longer history horizon of the difference scheme and a larger neighborhood.
The difference scheme (4) can also be viewed as a coupled map lattice or a special case of a cellular automaton with continuous state.
3. Problem formulation. Let us consider the difference equation (4) and let us assume that the minimized objective function is a functional and depends on the solution of this equation as follows: where y j is the state of the difference scheme at time j and g j (·) are given non-linear scalar functions. The whole spatiotemporal input signal U can be aggregated and presented as a matrix: The problem formulated and solved in this paper consists of finding a gradient G of an objective function J in a space of the spatiotemporal input signal: which is a matrix of partial derivatives: Once the gradient G is obtained it may be utilized during a gradient-descent iterative minimization procedure of the objective function appearing, for example, in parameter estimation or control optimization problems.

Problem solution.
To calculate the spatiotemporal gradient of the objective function J with respect to the whole input signal U we will utilize the concept of the adjoint system and adjoint variables known from the discrete formulation of the Pontryagin maximum principle.
Let us introduce a scalar function called Hamiltonian: where p j i , i = 1, 2, . . . , i max , j = 1, 2, . . . , j max are adjoint variables (or Lagrange multipliers), each corresponding to one scalar constraint (4) and satisfying the following equations of the adjoint system: Now, let us calculate p j i based on (10) and taking into account that y j i appears under the sum symbol in (10) three times: for i as y j i , for i − 1 as y j i+1 and for i + 1 as y j i−1 . This observation leads to the following equation of the adjoint system ADJOINT SENSITIVITY ANALYSIS OF A TUMOR GROWTH MODEL 1135 The adjoint system (12) should be solved backward in time with the initial condition Then the elements of the searched gradient (9) can be calculated as partial derivatives of the Hamiltonian with respect to u j i : ∂J 5. Results of gradient calculation -spatiotemporal sensitivity. The objective function J is defined as the total density of tumor cells at the final time j max : where j max = 200, i max = 200. Acquired from adjoint sensitivity analysis, the spatiotemporal gradient of the objective function (15) with respect to the radiation dose distribution d for model (3) is shown in Fig. 2. We have calculated this gradient twice, for two nominal spatiotemporal irradiation protocols: for zero irradiation (left panels) and for uniform irradiation d ≡ 1 (right panels). The results of simulation of the model compared to the antigradient (the gradient multiplied by −1) for particular time moments j = 1, 100, and 200 are shown in Fig. 3. In this Figure we can observe that the shape of the spatiotemporal sensitivity of J is different from the shape of the growing tumor. This is in contradiction, to some extent, with our intuition that J should be more sensitive for a higher tumor density. This phenomena and its possible reasons will be discussed further in next sections where the optimal spatiotemporal radiation will be presented.
In order to confirm that the gradient is properly calculated we compared it with the gradient computed by the finite difference method. In this method we change the radiation dose distribution only at one time j and one space point i by ∆d j i and check the resulting variation of the objective function ∆J. Then the gradient at time j and space point i can be approximated by the following formula: This has to be repeated i max × j max times to obtain the approximation of the whole gradient. Comparisons of the gradient calculated by the finite difference method and the adjoint sensitivity analysis for time j = 150 are shown in Figs. 4 and 5. We can also compare the time needed to execute both methods. On a laptop with Intel Core i7-4700MQ CPU and with implementation of the algorithms in MAT-LAB, the execution time that is needed to calculate the spatiotemporal gradient is about 2 minutes when using the finite difference method. Using the adjoint sensitivity analysis described in this paper we can calculate the whole spatiotemporal gradient in about 0.3 second. One can see that the adjoint sensitivity approach is especially useful for this type of task, especially when we take into account that during gradient-based therapy optimization the gradient has to be calculated many times.
The MATLAB code of the adjoint sensitivity analysis and the finite difference method used in this paper can be downloaded from the MathWorks File Exchange site. 6. Spatiotemporal radiation optimization. To optimize the radiation dose distribution d we used the gradient descent method with additional constraints: • Maximum cumulative dose S max : • Maximum point dose d max : To incorporate these constraints during optimization we modified the defined objective function J and the way in which the input signal of radiation dose distribution is generated. In order to implement the maximum cumulative dose constraint (17) we applied an external penalty function approach by adding an extra term to the objective function J in cases when the minimization algorithm exceeds S max : where k is a scalar positive coefficient increasing during the optimization. The maximum point dose constraint (18) was implemented by adding a new signal v which is an argument of a hyperbolic tangent function (having an upper limit equal to one). The value of this function multiplied by d max gives us the radiation dose distribution d at which all values do not exceed d max : After adding these modifications we calculated the spatiotemporal gradient of the modified objective function J (instead of J) with respect to the new signal v (instead of d) during the optimization procedure.   For the first set of parameters D = 27.7, ρ = 3.59 the values of the objective function J before (for radiation equals zero) and after optimization were 1.12e3 and 163.90, respectively. For the second set of parameters D = 8.9, ρ = 50.29 (the most aggressive tumor) the value of the objective function J before optimization was 3.30e7 and after optimization decreased only to 3.27e7. These results are summarized in Table 1. One can see that in the second case the radiation therapy has little effect on tumor progression. Further, especially for the second (the worse) case, the shape of the optimal spatiotemporal radiation dose distribution is visibly different than the tumor density profile. This raises an interesting question: what is the reason for such an effect? In general, the radiation is applied at points and times where the tumor density reaches average values (non-zero and not saturated -less than k t ). When the density of the tumor reaches the maximal environment capacity k t then the optimal irradiation level decreases. This is because of the form of the mathematical model of the tumor growth; the last term in (1), responsible for the killing action, disappears at the maximal density of the tumor and the optimal spatiotemporal radiation protocol is not applied to these places and times. Of course, this raises the question of the need to further improve the model used here, and in our future work we plan to apply our approach to different models and to assess the results.
Parameters D=27.7, ρ=3.59 D=8.9, ρ=50.29 J (d ≡ 0) 1.12e3 3.30e7 J (d opt ) 163.90 3.27e7 Table 1. Values of the objective function J without irradiation and with optimized radiation dose distribution for two sets of model parameters 7. Conclusions. In this work we propose a novel approach to optimization of spatiotemporal radiotherapy using a mathematical model of tumor growth described by partial differential equations. This allows to control the radiation dose which can be different for different points of the tumor and surrounding tissue and for different times. The optimization is based upon adjoint sensitivity analysis which gives information about the spatiotemporal sensitivity (being a function of time and space) of a given scalar objective function characterizing the model's solution. The original partial differential model is transformed into a discrete difference scheme and then the adjoint sensitivity analysis and the radiation protocol is optimized.
In future work, we plan to apply our approach to optimization of fractioned radiotherapy which is more realistic and applicable in clinical reality. We will also try to minimize other forms of the objective function by using the same adjoint sensitivity analysis.