An Hybrid Model for Rectal Tumour Response Prediction during Radiotherapy

Sena Apeke^1,2, Laurent Gaubert^1,3, Nicolas Boussion^1,4, Dimitris Visvikis^1,4, Olivier Saut⁵, Thierry Colin⁵, Philippe Lambin⁶, Vincent Rodin², Pascal Redou^1,3

¹LaTIM, Unité Mixte de Recherche 1101, Institut National de la Santé et de la Recherche Médicale, Brest, France.
²Immersive Natural User Interaction Team, Lab-STICC, Unité Mixte de Recherche CNRS 6285, Computer Science Department, Université de Bretagne Occidentale, Brest, France.
³Centre Européen de Réalité Virtuelle, Brest, France.
⁴Centre Hospitalier Régional Universitaire de Brest, Brest, France.
⁵Institut de Mathématiques de Bordeaux, Talence, France.
⁶The D-Lab & M-Lab, Department of Precision Medicine, GROW—School of Oncology and Developmental Biology, Maastricht University, Maastricht, The Netherlands.

DOI: 10.4236/ojbiphy.2022.124012                

Abstract

A hybrid model is proposed in this study to predict rectal tumour response during radiotherapy treatment. As the oxygen partial pressure distribution (pO₂) is a data which is naturally represented at the microscopic scale, we firstly estimate the optimal pO₂ distribution using both a diffusion equation and a discrete multi-scale model (that we proposed in a previous study). The aim is to use the effectiveness in algorithmic complexity of the discrete model and its multi-scale aspect in this work to estimate biological information at cellular scale and then construct them at macroscopic scale. Secondly, the obtained pO₂ distribution results are used as an input of a biomechanical model in order to simulate tumour volume evolution during radiotherapy. FDG PET images of 21 rectal cancer patients undergoing radiotherapy are used to simulate the tumour evolution during the treatment. The simulated results using the proposed hybride model, allow the interpretation of tumour aggressiveness.

Keywords

Tumour, Treatment, Response, Discrete, Density, FDG PET, SUV, PDE, Simulation

1. Introduction

The simulation of tumour growth and tumour response to radiation therapy remains a major research topic due to the overall impact of cancer [1] [2] . According to the World Health Organization (WHO), cancer burden rises to 18.1 million new cases and 9.6 million deaths worldwide in 2018. In this context, research on the treatment of malignant tumours is obviously a major concern and requires the mobilization of multidisciplinary research communities. In the literature, a wide series of mathematical or computational approaches to tumour growth modelling is proposed [3] - [8] . Such studies are generally theoretical and rarely tackle the individual specificities of real patients. Therefore, building relevant relationships between mathematical models and actual data like morphological or anatomical images corresponding to standard clinical use remains a real challenge. In addition to this macroscopic aspect, the microscopic scale of the tumour microenvironment should also be taken into account. For example, partial oxygen pressure (pO₂) is a very important local factor to consider when simulating tumour growth [9] [10] . Unfortunately, reliable and precise data concerning local pO₂ remain difficult to obtain for a given patient. The absence of such information makes it difficult to realistically evaluate simulation models. A potential solution is to build an autonomous model which calculates the oxygen partial pressure from available information and then uses it as an input for the main model, i.e. the one which predicts tumour growth. In the literature, several approaches to pO₂ modelling exist both at the microscopic scale [11] [12] and macroscopic scale [13] . Partial oxygen pressure is also of crucial importance for radiotherapy outcomes since hypoxic tumours, with a low level of oxygen, are known to be more resistant to radiation than non-hypoxic tumours. Ideally the dynamics and the heterogeneity of pO₂ distribution inside the tumour must be manageable. Furthermore, the impact of pO₂ and the response to radiation may differ according to the tumour cell types. The behaviour of the different types of tumour cells, either proliferative, hypoxic or quiescent, can be distinguished by the way they manage pO₂. In this context, taking oxygen partial pressure heterogeneity into consideration in tumour growth modelling during treatment is important, because the microenvironment has a significant impact on the efficacy of radiotherapy. Generally, pO₂ modelling does not take into account the dynamic nature of his evolution during the entire treatment [14] , it is considered constant across the whole tumour, or not taken into account at all [15] . We previously described a multi-scale approach allowing to take images and cellular cycle phases into account but pO₂ was considered constant inside the volume of interest and the tumour surface evolution was not addressed [16] . Typical cancer treatment generally consists of a combination of surgery, radiotherapy and chemotherapy. When radiotherapy is associated with the treatment, it is sometimes performed to reduce the size of the tumour before surgery [17] . After surgery, radiotherapy can still be used to reduce local recurrence risks. Surgery obviously requires careful planning, and indications leading to the least complicated therapeutic strategy are welcome. For this purpose, morphologic and metabolic images are used in order to get some characteristics of the tumour without necessarily going through biopsies. Furthermore, a CT scan (x-ray computed tomography) and FDG PET images (fluorodeoxyglucose positron emission tomography) are acquired before radiotherapy in order to delineate target and organs at risk, and also for optimizing treatment planning through relevant ballistics. In cases of potential tumour evolution during radiotherapy, additional images are acquired in order to adapt treatment according to the new size or shape of the tumour. In this context, the main objective of this study is to build a numerical tool able to predict the evolution of the tumour during radiotherapy by using images of the patient as an input. The proposed methodology is built on our previously described multi-scale and discrete framework [16] . In this new approach, a dynamic and heterogeneous map of pO₂ is generated and combined with a system of partial differential equations for modelling the tumour volume evolution. This hybrid process is evaluated using real FDG PET images acquired at different moments of radiotherapy in 21 rectal cancer cases.

2. Materials and Methods

2.1. Model Description

The new dynamic model for tumour response to radiotherapy that is proposed in this study comprises two main steps. As a first step, a diffusion equation is used to model the pO₂ evolution in the tumour at each time step [13] [18] . This equation is incorporated into a previously described multi-scale framework in order to take the personalized images of the patient into account [16] . As a second step, an advection reaction equation is proposed in order to predict the tumour volume evolution during radiotherapy treatment. These two steps compose the proposed hybrid approach and are described separately in the following sections.

2.1.1. pO₂ Evolution Modelling

Oxygen transport in tissues via blood vessels obviously depends on vessels structure but also on other biological constraints linked to the tumour behaviour. In the present study, Equation (1) was chosen to locally model the pO₂ evolution [18] :

$\frac{\partial p{O}_2}{\partial t}=\frac{2P_m}{R}\left(p_{cap}-p{O}_2\right)+\nabla \cdot \left(D\nabla p{O}_2\right)-c_{\max }\cdot \frac{p{O}_2}{p{O}_2+P_h}\cdot \rho$ (1)

where, $\rho$ represents the cell density in a voxel, with:

· $\frac{2P_m}{R}\left(p_{cap}-p{O}_2\right)$ : the source term, $P_m$ is the blood vessels permittivity and

R the radius;

· $\nabla \cdot \left(D\nabla p{O}_2\right)$ : the diffusion term, D is the diffusion coefficient;

· $c_{\max }\cdot \frac{p{O}_2}{p{O}_2+P_h}$ : the oxygen consumption per unit cell density, $c_{\max }$ is the maximum consumption of pO₂, and $P_h$ is the pO₂ at $\frac{c_{\max }}{2}$ .

In order to provide personalized simulation according to patient specific data coming from images, a parameter controlling the source term was introduced. By denoting $\mu$ this parameter, pO₂ distributions for a given patient are obtained by solving Equation (2). The values of the fixed parameters are given in Table 1.

$\frac{\partial p{O}_2}{\partial t}=\mu \left(p_{cap}-p{O}_2\right)+\nabla \cdot \left(D\nabla p{O}_2\right)-c_{\max }\cdot \frac{p{O}_2}{p{O}_2+P_h}\cdot \rho$ (2)

Since pO₂ consumption is a process that takes place at the cellular scale, a multi-scale approach is used in order to estimate the optimal distributions from macroscopic PET image data. The diagram shown in Figure 1 describes the operating steps of this multi-scale stochastic methodology, but full details can be

Table 1. List of parameters used in this work.

Figure 1. Diagram showing all the processes involved in the stochastic model, $\left[t_0\mathrm{,}{t}_f\right]$ is the time interval between the beginning and the end of a cell cycle. In this case $t_f=28\text{h}$ .

found in [16] . The model is based on transitions between the successive phases of the cellular cycle. The total number of active tumour cells ( $N_{cell}$ ) inside a voxel of the image is directly calculated from the voxel intensity ( $i_v$ ). Then, the number of tumour cells in each cellular cycle phase ( $N_{G_1},N_S,N_{G_2},N_M,N_{G_0}$ ) are deduced using a precomputed distribution ( $\lambda$ ).

2.1.2. Estimation of the Optimal pO₂ Distributions

Optimal pO₂ distributions are estimated by minimizing a cost function which depends on the number of cells given by simulated and real data. These numbers are calculated from clinical and simulated FDG PET [19] images obtained at the 8^th day following the beginning of radiotherapy. The simulated number of cells was derived from the multi-scale stochastic model. Let us now describe the optimization algorithm which consists of five steps:

1) an acceptability criterion is defined by means of a cost function F (Equation (3)); for every voxel index $\left(l\mathrm{,}n\mathrm{,}m\right)$ :

$F\left(P_1,P_2\right)\left(l,m,n\right)={\left(1-\frac{N_8^s\left(l,m,n\right)}{N_8^c\left(l,m,n\right)}\right)}^2+{\left(1-\frac{N_{15}^s\left(l,m,n\right)}{N_{15}^c\left(l,m,n\right)}\right)}^2$ (3)

where, $N_8^s$ and $N_{15}^c$ are the total numbers of tumour cells calculated respectively from simulated and clinical images. $P_1$ and $P_2$ are pO₂ distributions at day 0 and at day 8 respectively;

2) The capillary pressure is initialized at $t=t_0$ ( $t_0$ corresponds to day 0) as $P_{cap}=40\text{mmHg}$ (pO₂ in the arteries [20] );

3) Equation (2) is solved at each time step using finite difference method;

4) pO₂ distributions obtained in 3) are used as an input for the discrete model [16] ; Then the results are used to calculate the cost function 1);

5) If the result of the cost function is less than a set threshold, the value of the parameter $\mu$ and the pO₂ distributions are saved, and the algorithm is stoped. Otherwise, the parameter $\mu$ is modified using simulated annealing method [21] and the algorithm is resumed from step 2.

The final optimal pO₂ maps are used as input of a biomechanical model that we describe now.

2.1.3. Description of the Biomechanical Model

We denote by $A\subset R^3$ an image containing a patient tumour at time $t\in \left[\mathrm{0,}\Upsilon \right]$ , where $\Upsilon$ is the time between the beginning of the first image acquisition (before treatment) and the last acquisition (15 days after the beginning of the irradiation for 17 patients, or after radiotherapy and just before surgery for 4 patients). Then, we denote by $\Omega \left(t\right)$ , the tumour zone in the image $A$ (See Figure 2). $\Omega \left(t\right)$ is given by the set of standardized uptake values (SUV) calculated from FDG PET images [22] :

$\Omega \left(t\right)=\left\{\left(x\mathrm{,}y\mathrm{,}z\right)\in A\mathrm{<mo>\; /\; </mo>}\rho \left(x\mathrm{,}y\mathrm{,}z\mathrm{,}t\right)>0\right\}$ (4)

In this study, it is assumed that the tumour tissue is composed either of proliferative cells, quiescent cells, or necrotic cells [23] [24] . These cells densities

Figure 2. An example of data obtained from clinical images: (a) before the beginning of treatment $A_0$ , (b) after one week of treatment $A_8$ and (c) after two weeks of treatment $A_{15}$ . The corresponding tumour zones are respectively $\Omega \left(t_0\right)$ , $\Omega \left(t_8\right)$ , $\Omega \left(t_{15}\right)$ .

are denoted by ${\rho }_p\left(x\mathrm{,}t\right)$ , ${\rho }_q\left(x\mathrm{,}t\right)$ and ${\rho }_N\left(x\mathrm{,}t\right)$ respectively, ( $x=\left(x\mathrm{,}y\mathrm{,}z\right)\in \Omega \left(t\right)$ ). They satisfy an advection-reaction-diffusion equation [15] [25] :

$\frac{\partial {\rho }_l}{\partial t}+\nabla \cdot K=Sr\left({\rho }_l\right)+\nabla \cdot J-Tr\left({\rho }_l\right)$ (5)

where, $l=p,q,N$ and:

· Sr is the modelling tumour cells natural death and birth phenomena;

· $\nabla \cdot J$ models migratory movements; with $J=D\left(x\right)\nabla {\rho }_l$ , D is the diffusion coefficient;

· T models the cell death caused by radiotherapy;

· $\nabla \cdot$ is the divergence operator. $K=v\cdot {\rho }_l$ is the physical flow given by the system advection (all cells have the same advective velocity $v$ ).

Since rectal tumours are solid tumours, we assume that there is no migratory phenomenon caused by cells movement. Therefore, $J$ in Equation (5) vanishes. The only cell movements considered are the movements caused by the tumour volume variation. The source term Sr is modelled by a logistic function:

$S\left({\rho }_l\right)=\zeta \left(p{O}_2\right){\rho }_l\left(1-\frac{{\rho }_l}{\Phi }\right)$ (6)

with,

$\zeta \left(p{O}_2\right)=\sigma \cdot \frac{1+\tanh \left(p{O}_2-p{O}_2^h\right)}{2}$ (7)

where, $\sigma >0$ is the intrinsic growth rate, $\Phi$ is the maximum capacity of a voxel. $p{O}_2^h$ is the hypoxic threshold and tanh represents the classical hyperbolic tangent function. The term $\zeta$ gives a distinction between proliferating and quiescent cells, and was inspired by [26] . Cells survival probability after irradiation is given by the time-dependent linear quadratic model [27] .

$SF\left(p{O}_2\right)=\exp \left(-\alpha d\cdot p{O}_2\cdot \left(1+\iota \cdot \frac{\beta }{\alpha }\cdot d\cdot p{O}_2\right)\right)$ (8)

where, $\iota$ is an adjustment parameter of radiotherapy d accumulation while $\alpha$ and $\beta$ are classical radio-sensitivities parameters. Thus, tumour cells densities killed by irradiation are modelled as:

$T\left({\rho }_l\right)=\left(1-SF\left(p{O}_2\right)\right)\cdot {\rho }_l$ (9)

In this study, a macroscopic representation of the oxygen partial pressure in the tumour is considered. Indeed, according to the value of the pO₂ in the voxel, for methodological reasons, it exists either a density of proliferating and necrotic cells, or only a density of quiescent and necrotic cells. After normalization, we obtain the following Relation (10):

${\rho }_p\left(x\mathrm{,}t\right)+{\rho }_N\left(x\mathrm{,}t\right)=1\text{or}{\rho }_q\left(x\mathrm{,}t\right)+{\rho }_N\left(x\mathrm{,}t\right)=1$ (10)

As one can notice, Equation (5) is not closed because it has two unknown variables: cell density and velocity $v$ ( $K=v\cdot {\rho }_l$ ). To close this equation, it is assumed that tumour environment is isotropic and porous, allowing to use Darcy’s law:

$v=-\nabla \Pi$ (11)

where, $\Pi$ represents the local pressure. Based on the previous descriptions, tumour cells densities evolution equation is given by ( $\rho \in \left\{{\rho }_p\mathrm{,}{\rho }_q\right\}$ ):

$\{\begin{array}{l}\nabla \cdot \left(v\left(x\mathrm{,}t\right)\rho \left(x\mathrm{,}t\right)\right)=S\left(\rho \left(x\mathrm{,}t\right)\right)-T\left(\rho \left(x\mathrm{,}t\right)\right)\\ {\partial }_t{\rho }_N+\nabla \cdot \left(v{\rho }_N\right)=0\\ v\left(x\mathrm{,}t\right)=-\nabla \Pi \left(x\mathrm{,}t\right)\end{array}$ (12)

The determination of pressure $\Pi$ will lead to the knowledge of $v$ allowing to close the system (12). By summing the first two equations of this system and by using (10), we obtain:

$\nabla \cdot v=S\left(\rho \right)-T\left(\rho \right)$ (13)

Then by replacing $v$ by its expression given by (11) in (13), we obtain:

$-\Delta \Pi =S\left(\rho \right)-T\left(\rho \right)$ (14)

The system (15) summarizes all the equations needed for the simulation:

$\{\begin{array}{l}{\partial }_t\rho \left(x\mathrm{,}t\right)+\nabla \cdot \left(v\left(x\mathrm{,}t\right)\rho \left(x\mathrm{,}t\right)\right)=S\left(\rho \left(x\mathrm{,}t\right)\right)-T\left(\rho \left(x\mathrm{,}t\right)\right)\\ {\partial }_t{\rho }_N+\nabla \cdot \left(v{\rho }_N\right)=0\\ v\left(x\mathrm{,}t\right)=-\nabla \Pi \left(x\mathrm{,}t\right)\\ -\Delta \Pi =S\left(\rho \right)-T\left(\rho \right)\\ \Pi \left(x\mathrm{,}t\right)=0\\ \rho \left(x\mathrm{,0}\right)={\rho }_0\left(x\right)\end{array}$ (15)

2.2. Simulation of the Biomechanical Model

2.2.1. Meshing and Simulation Algorithm

For simulation purposes we used a finite volumes based method and a 3D cartesian grid $\left[\mathrm{0,}{I}_x\right]\times \left[\mathrm{0,}{I}_y\right]\times \left[\mathrm{0,}{I}_z\right]$ , where $I_x$ , $I_y$ and $I_z$ are voxels numbers in the x, y and z directions, respectively (see Figure 3). This grid corresponds to the distribution of cells densities in the SUV medical images. The simulation of the first equation in system (15) was performed by using the Strang splitting method [28] [29] and then by determining the fields of proliferating and quiescent cells densities. For this latter purpose, Equations (16) and (17) are simulated:

$\frac{\partial \rho }{\partial t}\left(x\mathrm{,}t\right)+\nabla \left(\rho \left(x\mathrm{,}t\right)v\left(x\mathrm{,}t\right)\right)=0$ (16)

$\frac{\text{d}\rho }{\text{d}t}\left(x\mathrm{,}t\right)=S\left(\rho \left(x\mathrm{,}t\right)\right)-T\left(\rho \left(x\mathrm{,}t\right)\right)$ (17)

In practice, the two equations above are simulated by running the following algorithm:

1) Equation (14) is used to determine the pressure field $\Pi$ ;

2) Equation (11) is used to compute the velocity field $v$ , knowing the pressure field $\Pi$ ;

3) Proliferating ( ${\rho }_p$ ) and quiescent ( ${\rho }_q$ ) cells densities are computed by simulating Equations (16) and (17);

4) If necessary, necrotic cells densities were computed as follows: ${\rho }_N=1-{\rho }_p-{\rho }_q$ .

The discretization and numerical schemes used are presented in the Subsection 6.1.

2.2.2. Evaluation of the Model

The proposed model (Equation (15)) contains three parameters: $\sigma$ , $\Phi$ and $\iota$ , which are analyzed using Sobol sensitivity indices [30] .

The tumour volume at time t is calculated by:

$V\left(t\right)={\displaystyle {\int }_{\Omega \left(t\right)}}{\text{1}}_{\left\{x\mathrm{<mo>\; /\; </mo>}\rho \left(t\mathrm{,}x\right)>0\right\}}\text{d}x$ (18)

where, ${\text{1}}_{\left\{x\mathrm{<mo>\; /\; </mo>}\rho \left(t\mathrm{,}x\right)>0\right\}}$ is the indicator function defined on $\left\{x\mathrm{<mo>\; /\; </mo>}\rho \left(t\mathrm{,}x\right)>0\right\}$ .

The correlation formula for clinical and simulated volumes comparison is:

$Corr\left(\mathrm{<mi>\; \%\; </mi>}\right)=\left(1-\frac{\left|V_s-V_c\right|}{V_c}\right)\cdot 100$ (19)

where, $Corr$ is the correlation result, $V_s$ and $V_c$ are the simulated and clinical volumes, respectively.

Figure 3. Meshing of the study domain.

3. Results

A clinical database containing 17 patients is used to evaluate the proposed model. First, optimal pO₂ distributions are estimated using the diffusion equation and the stochastic multi-scale model, as explained above. Second, the obtained pO₂ results are used as an input to the biomechanical model in order to simulate tumour volume evolution during radiotherapy. FDG PET images are used for both the pO₂ estimation and the tumour volume evolution for each patient.

3.1. Estimation of the pO₂ Distributions

The specific data for patient 5 and 12 are given as examples of obtained results. For these two patients, values for the parameter $\mu$ are 0.13 and 0.8 respectively. The distribution of oxygen partial pressure over time is given in Figure 4 and Figure 5. In these figures, (a) and (c) show a cross-section of the pO₂ images at day 8 and day 15 after the beginning of irradiation; (b) and (d) are histograms over the whole volume of each of these pO₂ images. The estimation of pO₂ distributions for the whole set of patients is shown in Table 2. From this table one can see that oxygen partial pressure is increased for almost all patients after one week of treatment. This is an expected result but a validation cannot be put forward in the absence of a real pO₂ measurement. Nevertheless, the results suggest that the evolution of pO₂ should be taken into account for the tumour growth simulation, supporting published data [31] .

Figure 4. (a) and (c) sections of pO₂ distribution images, obtained at the 8^th and 15^th days, (b) and (d) histograms of these images for the patient 5 example.

Figure 5. (a) and (c) sections of pO₂ distribution images, obtained at the 8^th and 15^th days, (b) and (d) histograms of these images of the patient 12 example.

Table 2. pO₂ optimal distributions results, at 8^th (day 8) and 15^th days (day 15) after the start of treatments. $\overline{p{O}_2}8$ and $\overline{p{O}_2}15$ are the pO₂ average values, respectively at the 8^th and 15^th days.

3.2. Simution of the Tumour Volume Evolution

According to the value of pO₂ in a voxel, the system of Equation (15) is simulated by considering only the proliferating and necrotic cells or only the necrotic and quiescent cells in the voxel. Sensitivities of biomechanical model parameters according to Sobol analysis are presented in Table 3. The results show that, contrary to parameters $\sigma$ and $\iota$ , the parameter $\Phi$ can be fixed because a small disturbance of the latter does not influence the output of the model. Throughout the simulation and for all patients, we fixed $\Phi =2$ . The two other parameters are estimated using the annealing method [21] for each patient. As an illustration, tumour volume evolution during the treatment is showed, Figure 6 illustrates the case of patient number 5 (and Figure 7 shows a 2D section). Also, the

Figure 6. Tumour volume evolution, case of patient 5.

Figure 7. (a) Clinical SUV image and (b) Simulated SUV image, case of patient 5.

optimal adjustment parameters obtained for patients 5 and 12 are given in Table 4. The aim of Table 5 is to give an overview of the correlations (using Equation (19)) between simulated and clinical images at day 8 and at day 15 after the beginning of radiotherapy. For day 15 one can observe that 10/17 of patients have a correlation superior to 90%, 4/17 of patients have an average correlation of 70% - 80%, 2/17 of patient have an average correlation of 60% - 70%, and 1/17 has an overall correlation <60.

Table 3. Sobol total sensitivity indices for the biomechanical model parameters.

Table 4. Optimal parameters obtained for patients 5 and 12.

Table 5. Results of correlations between clinical and simulated images at 8^th and 15^th days after the beginning of the irradiation.

In order to compare the aggressiveness of tumours, the example of patient 12 is also given (see Figure 9 and Figure 10). Additionally, Figures 8-11 give a 3D representation of the simulation results for patients 5 and 12 respectively.

One of the objectives of this study is to build a model able to propose the potential extent of the tumour a few days before surgery in order to help resection planning. To this end, the model is also applied to four other patients who have one additional PET image a few days prior to surgery. Results can be seen on Table 6 where correlations between simulated and clinical images are presented at day 8, day 15 and day 90 after the beginning of radiotherapy. Except for patient 19 who has a low correlation <50%, the other three have a correlation >90%.

Figure 8. 3D illustration of clinical and simulated volumes, case of patient 5.

Figure 9. Tumour volume evolution, case of patient 12.

Figure 10. (a) Clinical SUV image and (b) Simulated SUV image, case of patient 12.

Figure 11. 3D illustration of clinical and simulated volumes, case of patient 12.

Table 6. Correlations results, between clinical and simulated images at the 8^th (day 8), 15^th (day 15) and 90^th (day 90) days after the beginning of irradiation.

4. Discussion

The impact of pO₂ remains one of the most studied biological phenomena in simulations of tumour growth and tumour response to radiotherapy [32] [33] . Most of the latest studies focus either on space-time aspects via reaction-diffusion equations [34] or on purely biological aspects at the cellular level [35] . These algorithms are generally compared with each other [36] or challenged against theoretical models [37] or empirical information [38] , but they are rarely confronted with real data. In the present work, we proposed a method based on fully-described reaction-diffusion equations driven by a cellular-based stochastic approach [16] . The obtained platform can be adapted to available clinical data and is evaluated by using a series of FDG PET images from 21 patients. The distribution of pO₂ inside the tumour is known to have an impact on radiotherapy effects, but it remains diffcult to obtain personalized and reliable pO₂ information from images. Often, functional imaging during cancer monitoring only consists of FDG PET images depicting glucose metabolism. For this reason, we separated our simulation approach into two parts working concomitantly. The first part is a multi-scale model which simulates temporal and spatial evolution of oxygen concentration from available patient images. The second part simulates the tumour growth using a biomechanical approach based on the reaction-diffusion equations and the pO₂ knowledge as provided by the first step.

Concerning this second part, the source term is a logistic function whose growth rate depends on the oxygen partial pressure. This made the entire reaction term dependent on pO₂, leading to an increased complexity. The challenge is to find the oxygen distribution that would provide the best prediction in terms of tumour volume evolution. This hybrid methodology is applied to a set of real data and the obtained results appeared satisfactory since a reasonably good agreement is observed between real and computed data (see Table 5). In addition to obtain information on the tumour volume evolution during radiotherapy, it is worth mentioning that the model also suggests a qualitative overview of the tumour aggressiveness. As an illustration, one can observe that the tumour of patient 5 is found more radioresistant than that of patient 12 (see Figures 7-10) while in the same time the simulated levels of oxygen partial pressure for patient 5 are found lower than those of patient 12 (Table 2). This difference in aggressiveness is also markedly observed before the beginning of treatments since tumour 5 appears to grow very rapidly while tumour 12 appears constant in volume. This general observation is in accordance with the obtained intrinsic growth rates, with $\sigma$ values of 0.028 for patient 5 and 0.00026 for patient 12 (Table 4).

The above encouraging results suggest that this study will benefit from a validation on a larger database, including a comprehensive clinical follow-up of the patients. However, despite the variation of tumour cell densities between voxels, as driven by the computed flux, this model does not directly take into account tumour deformations. This is illustrated by the fact that according to Dice

Table 7. Dice index of simulated images on the 90^th day after the beginning of irradiation, cases of patients 18, 19, 20 and 21.

metric [39] the proposed hybrid model does not predict the tumour shape evolution accurately. Also known as an overlap index, the Dice metric is a similarity index that is widely used for volume comparison purposes. The values of Dice index for patients 18, 19, 20 and 21 are given in Table 7 and show that simulated and real images at day 90 do not overlap (Dice < 0.6).

The mechanical constraints that are at the origin of the shape alteration are indeed difficult to control. This point clearly represents a potential improvement of our method, and in this context, it will be valuable to use anatomical information provided by CT or MRI images for example.

5. Conclusion

In this study, we proposed a methodology for simulating tumour growth and tumour response to radiation therapy. The adopted approach is based on the synergy between a discrete multi-scale stochastic model and a continuous model based on advection-reaction equations. This image-based process can be personalised according to available clinical data. The evaluation of the method on actual FDG images of patients suffering from rectal cancer is encouraging and opens several opportunities for improvement. The use of multi-modal images providing additional functional information instead of the single modality as presented here will certainly reinforce the robustness and the reliability of the simulations. Also, the introduction of morphological images like X-ray computed tomography is expected to help manage the mechanical constraints that can modify the shape of the tumour and influence its deformation.

Acknowledgements

The authors would like to thank the Bretagne Region and the LaTIM for their financial support.

Appendix

Discretization and Numerical Schemes

Denote by ${\rho }_{i\mathrm{,}j\mathrm{,}k}^n=\rho \left(\left(x_i\mathrm{,}{y}_j\mathrm{,}{z}_k\right)\mathrm{,}{t}_n\right)$ , the local cells densities at time $t_n$ and at the position $\left(x_i\mathrm{,}{y}_j\mathrm{,}{z}_k\right)$ , by $v_{i\mathrm{,}j\mathrm{,}k}^n=v\left(\left(x_i\mathrm{,}{y}_j\mathrm{,}{z}_k\right)\mathrm{,}{t}_n\right)$ and ${\Pi }_{i\mathrm{,}j\mathrm{,}k}^n=\Pi \left(\left(x_i\mathrm{,}{y}_j\mathrm{,}{z}_k\right)\mathrm{,}{t}_n\right)$ the corresponding advection velocity and local pressure respectively. The numerical approaches used for the simulation are as follows:

· A finite difference method based on an implicit Euler scheme is used for simulating Equation (14);

· A finite volume method with the 5th order WENO scheme (Weighted Essentially Non-Oscillatory [40] ) is used for simulation of Equation (16). By rewriting this equation in the form:

${\partial }_t\rho +{\partial }_x\left(v_x\rho \right)+{\partial }_y\left(v_y\rho \right)+{\partial }_z\left(v_z\rho \right)=0$ (20)

and by integrating it, we obtain:

${\bar{\rho }}_{ijk}^{n+1/2}={\bar{\rho }}_{ijk}^n-\frac{\Delta t}{\Delta x}\left(F_{i+1/2}-F_{i-1/2}\right)-\frac{\Delta t}{\Delta y}\left(G_{j+1/2}-G_{j-1/2}\right)-\frac{\Delta t}{\Delta z}\left(H_{k+1/2}-H_{k-1/2}\right)$ (21)

where, ${\bar{\rho }}_{ijk}^n$ is the mean local tumour cells densities, $F_i$ , $G_j$ and $H_k$ , are their numerical flow, respectively in the x, y and z directions. For all i, j and k, we wrote $\Delta x_i=\Delta x$ , $\Delta y_j=\Delta y$ , and $\Delta z_k=\Delta z$ , with $\Delta x$ , $\Delta y$ , and $\Delta z$ representing voxel dimensions in the x, y and z directions. $\Delta t$ is the time step;

· A finite differences method based on an explicit Euler scheme is used for simulating Equation (17), with a discretization given by:

${\bar{\rho }}_{ijk}^{n+1}={\bar{\rho }}_{ijk}^{n+1/2}+\Delta t\cdot \left(S_{i\mathrm{,}j\mathrm{,}k}^{n+1/2}-T_{i\mathrm{,}j\mathrm{,}k}^{n+1/2}\right)$ (22)

where, $S_{i\mathrm{,}j\mathrm{,}k}^{n+1/2}=S\left(\bar{\rho }\left(x_i\mathrm{,}{y}_j\mathrm{,}{z}_k\mathrm{,}t+1/2\right)\right)$ , and $T_{i\mathrm{,}j\mathrm{,}k}^{n+1/2}=T\left(\bar{\rho }\left(x_i\mathrm{,}{y}_j\mathrm{,}{z}_k\mathrm{,}t+1/2\right)\right)$ .

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Bekker, R.A., Kim, S., Pilon-Thomas, S. and Enderling, H. (2022) Mathematical Modeling of Radiotherapy and Its Impact on Tumor Interactions with the Immune System. Neoplasia, 28, Article ID: 100796. https://doi.org/10.1016/j.neo.2022.100796
[2]	Borkenstein, K., Levegrn, S. and Peschke, P. (2004) Modeling and Computer Simulations of Tumor Growth and Tumor Response to Radiotherapy. Radiation Research, 162, 71-83. https://doi.org/10.1667/RR3193
[3]	Ribba, B., Colin, T. and Schnell, S. (2006) A Multiscale Mathematical Model of Cancer, and Its Use in Analyzing Irradiation Therapies. Theoretical Biology and Medical Modelling, 3, 1-19. https://doi.org/10.1186/1742-4682-3-7
[4]	Lo, C.F. (2007) Stochastic Gompertz Model of Tumour Cell Growth. Journal of Theoretical Biology, 248, 317-321. https://doi.org/10.1016/j.jtbi.2007.04.024
[5]	Iwata, K., Kawasaki, K. and Shigesada, N. (2000) A Dynamical Model for the Growth and Size Distribution of Multiple Metastatic Tumors. Journal of Theoretical Biology, 203, 177-186. https://doi.org/10.1006/jtbi.2000.1075
[6]	Titz, B. and Jeraj, R. (2008) An Imaging-Based Tumour Growth and Treatment Response Model: Investigating the Effect of Tumour Oxygenation on Radiation Therapy Response. Physics in Medicine & Biology, 53, 4471-4488. https://doi.org/10.1088/0031-9155/53/17/001
[7]	Titz, B., Kozak, K.R. and Jeraj, R. (2012) Computational Modelling of Anti-Angiogenic Therapies Based on Multiparametric Molecular Imaging Data. Physics in Medicine & Biology, 57, 6079-6101. https://doi.org/10.1088/0031-9155/57/19/6079
[8]	Gerisch, A. and Chaplain, M.A.J. (2008) Mathematical Modelling of Cancer Cell Invasion of Tissue: Local and Non-Local Models and the Effect of Adhesion. Journal of Theoretical Biology, 250, 684-704. https://doi.org/10.1016/j.jtbi.2007.10.026
[9]	Vaupel, P. (2004) Tumor Microenvironmental Physiology and Its Implications for Radiation Oncology. Seminars in Radiation Oncology, 14, 198-206. https://doi.org/10.1016/j.semradonc.2004.04.008
[10]	Chen, Y.J., Cairns, R., Papandreou, I., Koong, A. and Denko, N.C. (2009) Oxygen Consumption Can Regulate the Growth of Tumors, a New Perspective on the Warburg Effect. PLOS ONE, 4, e7033. https://doi.org/10.1371/journal.pone.0007033
[11]	Espinoza, I., Peschke, P. and Karger, C.P. (2013) A Model to Simulate the Oxygen Distribution in Hypoxic Tumors for Different Vascular Architectures. Medical Physics, 40, Article ID: 081703. https://doi.org/10.1118/1.4812431
[12]	Grimes, D.R., Kannan, P., Warren, D.R., Markelc, B., Bates, R., Muschel, R. and Partridge, M. (2016) Estimating Oxygen Distribution from Vasculature in Three-Dimensional Tumour Tissue. Journal of the Royal Society Interface, 13, Article ID: 20160070. https://doi.org/10.1098/rsif.2016.0070
[13]	Daşu, A., Toma-Daşu, I. and Karlsson, M. (2003) Theoretical Simulation of Tumour Oxygenation and Results from Acute and Chronic Hypoxia. Physics in Medicine & Biology, 48, 2829-2842. https://doi.org/10.1088/0031-9155/48/17/307
[14]	Powathil, G.G., Gordon, K.E., Hill, L.A. and Chaplain, M.A.J. (2012) Modelling the Effects of Cell-Cycle Heterogeneity on the Response of a Solid Tumour to Chemotherapy: Biological Insights from a Hybrid Multiscale Cellular Automaton Model. Journal of Theoretical Biology, 308, 1-19. https://doi.org/10.1016/j.jtbi.2012.05.015
[15]	Mi, H.M., Petitjean, C., Dubray, B., Vera, P. and Ruan, S. (2014) Prediction of Lung Tumor Evolution during Radiotherapy in Individual Patients with PET. IEEE Transactions on Medical Imaging, 33, 995-1003. https://doi.org/10.1109/TMI.2014.2301892
[16]	Apeke, S., Gaubert, L., Boussion, N., Lambin, P., Visvikis, D., Rodin, V. and Redou, P. (2018) Multi-Scale Modeling and Oxygen Impact on Tumor Temporal Evolution: Application on Rectal Cancer during Radiotherapy. IEEE Transactions on Medical Imaging, 37, 871-880. https://doi.org/10.1109/TMI.2017.2771379
[17]	Gérard, A., Buyse, M., Nordlinger, B., Loygue, J., Pène, F., Kempf, P., Bosset, J.F., Gignoux, M., Arnaud, J.P. and Desaive, C. (1988) Preoperative Radiotherapy as Adjuvant Treatment in Rectal Cancer. Final Results of a Randomized Study of the European Organization for Research and Treatment of Cancer (EORTC). Annals of Surgery, 208, 606-614. https://doi.org/10.1097/00000658-198811000-00011
[18]	Toma-Daşu, I. and Daşu, A. (2013) Modelling Tumour Oxygenation, Reoxygenation and Implications on Treatment Outcome. Computational and Mathematical Methods in Medicine, 2013, Article ID: 141087. https://doi.org/10.1155/2013/141087
[19]	Gambhir, S.S., Czernin, J., Schwimmer, J., Silverman, D.H., Coleman, R.E. and Phelps, M.E. (2001) A Tabulated Summary of the FDG PET Literature. Journal of Nuclear Medicine, 42, 1S-93S.
[20]	Bertout, J.A., Patel, S.A. and Simon, M.C. (2008) The Impact of O2 Availability on Human Cancer. Nature Reviews Cancer, 8, 967-975. https://doi.org/10.1038/nrc2540
[21]	Kirkpatrick, S., Gelatt, C.D. and Vecchi, M.P. (1983) Optimization by Simulated Annealing. Science, 220, 671-680. https://doi.org/10.1126/science.220.4598.671
[22]	Kim, C.K., Gupta, N.C., Chandramouli, B. and Alavi, A. (1994) Standardized Uptake Values of FDG: Body Surface Area Correction Is Preferable to Body Weight Correction. Journal of Nuclear Medicine, 35, 164-167.
[23]	Ferreira, S.C., Martins, M.L. and Vilela, M.J. (2002) Reaction-Diffusion Model for the Growth of Avascular Tumor. Physical Review E, 65, Article ID: 021907. https://doi.org/10.1103/PhysRevE.65.021907
[24]	Wallace, D.I. and Guo, X.Y. (2013) Properties of Tumor Spheroid Growth Exhibited by Simple Mathematical Models. Frontiers in Oncology, 3, Article No. 51. https://doi.org/10.3389/fonc.2013.00051
[25]	Konukoglu, E., Sermesant, M., Clatz, O., Peyrat, J.-M., Delingette, H. and Ayache, N. (2007) A Recursive Anisotropic Fast Marching Approach to Reaction Diffusion Equation: Application to Tumor Growth Modeling. 20th International Conference, IPMI 2007, Kerkrade, 2-6 July 2007, 687-699.
[26]	Cornelis, F., Saut, O., Cumsill, P., Lombardi, D., Iollo, A., Palussiere, J. and Colin, T. (2013) La modélisation mathématique in Vivo de la croissance tumorale sur les données de l’imagerie: Un avenir proche? Journal de Radiologie Diagnostique et Interventionnelle, 94, 610-617. https://doi.org/10.1016/j.jradio.2013.02.008
[27]	Brenner, D.J., Hlatky, L.R., Hahnfeldt, P.J., Huang, Y. and Sachs, R.K. (1998) The Linear-Quadratic Model and Most Other Common Radiobiological Models Result in Similar Predictions of Time-Dose Relationships. Radiation Research, 150, 83-91. https://doi.org/10.2307/3579648
[28]	Castella, F., Chartier, P., Descombes, S. and Vilmart, G. (2009) Splitting Methods with Complex Times for Parabolic Equations. BIT Numerical Mathematics, 49, 487-508. https://doi.org/10.1007/s10543-009-0235-y
[29]	Glowinski, R. and Le Tallec, P. (1989) Augmented Lagrangian and Operator-Splitting Methods in Nonlinear Mechanics. Society for Industrial and Applied Mathematics, University City. https://doi.org/10.1137/1.9781611970838
[30]	Nossent, J., Elsen, P. and Bauwens, W. (2011) Sobol’ Sensitivity Analysis of a Complex Environmental Model. Environmental Modelling & Software, 26, 1515-1525. https://doi.org/10.1016/j.envsoft.2011.08.010
[31]	Lyng, H., Sundfør, K. and Rofstad, E.K. (2000) Changes in Tumor Oxygen Tension during Radiotherapy of Uterine Cervical Cancer: Relationships to Changes in Vascular Density, Cell Density, and Frequency of Mitosis and Apoptosis. International Journal of Radiation Oncology, Biology, Physics, 46, 935-946. https://doi.org/10.1016/S0360-3016(99)00497-6
[32]	Petit, S.F., Lambin, P., Seigneuric, R., Murrer, L., Riel, N.A.W., van Dekker, A.L.J. and Wouters, B.G. (2007) Realistic Cellular Oxygen Model Reveals High Potential Gains with Dose Painting and Offers a General Solution for Incorporating Cellular Distributions That Underlie Imaging Data. International Journal of Radiation Oncology, Biology, Physics, 69, S69. https://doi.org/10.1016/j.ijrobp.2007.07.126
[33]	Harriss-Phillips, W.M., Bezak, E. and Potter, A. (2016) Stochastic Predictions of Cell Kill during Stereotactic Ablative Radiation Therapy: Do Hypoxia and Reoxygenation Really Matter? International Journal of Radiation Oncology, Biology, Physics, 95, 1290-1297. https://doi.org/10.1016/j.ijrobp.2016.03.014
[34]	Roque, T., Risser, L., Kersemans, V., Smart, S., Allen, D., Kinchesh, P., Gilchrist, S., Gomes, A.L., Schnabel, J.A. and Chappell, M.A. (2018) A DCE-MRI Driven 3-D Reaction-Diffusion Model of Solid Tumor Growth. IEEE Transactions on Medical Imaging, 37, 724-732. https://doi.org/10.1109/TMI.2017.2779811
[35]	Zangooei, M.H. and Habibi, J. (2017) Hybrid Multiscale Modeling and Prediction of Cancer Cell Behavior. PLOS ONE, 12, e0183810. https://doi.org/10.1371/journal.pone.0183810
[36]	Elaff, I. (2018) Comparative Study between Spatio-Temporal Models for Brain Tumor Growth. Biochemical and Biophysical Research Communications, 496, 1263-1268. https://doi.org/10.1016/j.bbrc.2018.01.183
[37]	Sundstrom, A., Bar-Sagi, D. and Mishra, B. (2016) Simulating Heterogeneous Tumor Cell Populations. PLOS ONE, 11, e0168984. https://doi.org/10.1371/journal.pone.0168984
[38]	Forster, J.C., Douglass, M.J.J., Harriss-Phillips, W.M. and Bezak, E. (2017) Development of an in Silico Stochastic 4D Model of Tumor Growth with Angiogenesis. Medical Physics, 44, 1563-1576. https://doi.org/10.1002/mp.12130
[39]	Taha, A.A. and Hanbury, A. (2015) Metrics for Evaluating 3D Medical Image Segmentation: Analysis, Selection, and Tool. BMC Medical Imaging, 15, Article No. 29. https://doi.org/10.1186/s12880-015-0068-x
[40]	Jiang, G.-S. and Wu, C.-C. (1999) A High-Order WENO Finite Difference Scheme for the Equations of Ideal Magnetohydrodynamics. Journal of Computational Physics, 150, 561-594. https://doi.org/10.1006/jcph.1999.6207