A radiobiological model of radiotherapy response and its correlation with prognostic imaging variables

Abstract

Radiobiological models of tumour control probability (TCP) can be personalized using imaging data. We propose an extension to a voxel-level radiobiological TCP model in order to describe patient-specific differences and intra-tumour heterogeneity. In the proposed model, tumour shrinkage is described by means of a novel kinetic Monte Carlo method for inter-voxel cell migration and tumour deformation. The model captures the spatiotemporal evolution of the tumour at the voxel level, and is designed to take imaging data as input. To test the performance of the model, three image-derived variables found to be predictive of outcome in the literature have been identified and calculated using the model's own parameters. Simulating multiple tumours with different initial conditions makes it possible to perform an in silico study of the correlation of these variables with the dose for 50% tumour control ([Formula: see text]) calculated by the model. We find that the three simulated variables correlate with the calculated [Formula: see text]. In addition, we find that different variables have different levels of sensitivity to the spatial distribution of hypoxia within the tumour, as well as to the dynamics of the migration mechanism. Finally, based on our results, we observe that an adequate combination of the variables may potentially result in higher predictive power.

Figure 1

Schematic of the key biological processes of the model. Two neighbouring tumourlets with arbitrary compartment sizes are depicted at three different time points. Cells in the Proliferative comparment (P) are normoxic and proliferating. Cells in the Intermediately hypoxic compartment (I) are hypoxic and non-proliferating, but still metabolizing. Cells in the extremely Hypoxic compartment (H) are starving and non-proliferating. Steps 1 and 2 include processes that happen only within each tumourlet, following the method described in Jeong et al. 2013. In step 1, a fraction of the cells in each compartment gets ‘doomed’ due to radiation. After some of the doomed cells die after attempting mitosis, the remaining cells move to more nutrient-rich compartments in step 2. Step 3 represents the inter-tumourlet migration introduced in the present study. The figures are meant to be suggestive; the actual arrangement of various cells is highly variable.

Figure 2. Details of the simulation process

(a) Simplified schematic of the main steps of the rfKMC algorithm. The goal is to produce a succession of transitions with the right frequencies. In order to do this, the algorithm builds the cumulative sum vector of all transition rates (that is, [r₁₂, r₁₂ + r₁₃,…, R], where R = Σ_i≠jr_ij. Then, a random number uniformly distributed between 0 and R is generated. The number will fall between two elements of the vector: this determines which transition occurs next. Because likely transitions have large transition rates, it is more likely that a random number will fall within their range. This is how the right dynamics is obtained. (b) Flow chart of the full simulation process. ‘Cell cycle’ indicates mitosis and cell death, and ‘reclassify’ refers to the process in which cells move to a compartment with more nutrients if space is available. The ellipsis indicates that the flow lines corresponding to the remaining tumourlets have been omitted.

Figure 3

Sample maps of the effective transition rates used as input to the rfKMC algorithm, overlayed on contour plots of the number of tumour cells in the I compartment (left) and the total number of cells (right). Red lines indicate regions of higher number of intermediately hypoxic tumour cells in the left panel, and higher total number of cells in the right panel. Each effective transition rate (i.e. each arrow) is obtained by adding vectorially all the transition rates originating from the same tumourlet and ending in its nearest and next-to-nearest neighbours.

Figure 4

Two examples of the profile distribution of I compartment sizes for two simulations with different values of σ_I and I_max. Compartment sizes are normalized to 1.

Figure 5

Time evolution of the transverse profile of the number of tumour cells (left) and of the maximum FDG uptake (right) of a simulated tumour with I_max=45% and σ_I = 0.3, both with and without cell migration. The numbers in the left panel indicate the number of days passed.

Figure 6

Evolution of a simulated tumour with I_max=45% and σ_I = 0.3, with (left) and without (right) cell migration. The hypoxic core is off-centre in this particular simulation to better illustrate the properties of the simulation. Each panel contains the evolution in time of the total number of hypoxic tumour cells, the number of stromal cells, and the simulated FDG uptake. Dark blue equals zero density.

Figure 7

Correlation of several imaging-based variables with TCD₅₀: δFDG (left), $f_P^{(10\%)}$ (centre) and M (right). Each marker corresponds to a simulated tumour with different initial conditions of hypoxia, defined in terms of I_max and σ_I.

Figure 8

Correlation of (a) δFDG and (b) M with TCD₅₀, plotted as a function of σ_I and colour-coded according to I_max. The simulations do not include migration.

Figure 9

Three-dimensional scatter plot of δFDG, $f_P^{(10\%)}$ and M for the simulated tumours. Points falling in the region where δFDG is approximately linearly correlated with TCD₅₀ are marked with a blue circle. Similarly, points falling in the region where $f_P^{(10\%)}$ and M are linearly correlated with TCD₅₀ are marked with green squares and red crosses, respectively. The ranges are δFDG < 1, $f_P^{(10\%)}<0.5$, and M< 6, and they are determined from visual inspection of the correlation plots between each of the variables and TCD₅₀.