Abstract
Background: High-grade gliomas with a widespread infiltration beyond the lesion detectable on diagnostic images are increasingly treated with Gamma Knife™ Radiosurgery (GKRS). The aim of this study was to assess the cell infiltration impact on the GKRS outcome for invasive gliomas. Materials and Methods: Tumor cell distribution was predicted using a novel algorithm whose computations are iterated until they reach an agreement with histopathology results. Treatment plans with different combinations of dose prescription (20 Gy at 50%-20% isodose) and targets [Gross Tumour Volume (GTV), zone 1 with 100%-60% of the GTV cell density and zone 2 with 60%-0% of the GTV cell density] were evaluated using standard conformity indexes (CI) and radiobiological parameters. Results: Considerable differences in terms of tumor control probability were found between plans having similar CI but different targets. Conclusion: To account for tumor cell infiltration outside the target is of key importance in GKRS and a radiobiological evaluation should accompany well-established CI.
High-grade gliomas (grades III and IV) are brain tumours with widespread infiltration in the surrounding brain tissue (1). The treatment of this aggressive type of tumour usually involves surgery, radiotherapy, and chemotherapy (2, 3). Gamma Knife™ Radiosurgery (GKRS) alone or in conjunction with other techniques, has also been used for salvage treatment of recurrent high-grade gliomas (4, 5).
The steep dose fall-off of the radiosurgery treatment and the definition of the target as a lesion detectable in the diagnostic images [Gross Tumour Volume (GTV) without margins] make the presence of the microscopic tumour spread a matter of concern when dealing with high-grade gliomas (6).
With the imaging techniques currently available, regions with a density of infiltrated cells lower than a certain threshold, such as 8,000 cells/mm3 for enhanced Magnetic Resonance Imaging (MRI) (7), are not visible.
The maximal extent of the microscopic spread of cancerous cells varies in density while its maximal extent in the peripheral region has been estimated to be 3-4 cm away from the tumour boundary (8, 9). The region of infiltration has a variable cell density, which decreases on average with increasing distance from the tumour core. Yamahara et al. (10) have identified three areas of infiltration, with cell densities equal to 100-60% (zone A), 60-20% (zone B), and 20-0% (zone C) of the tumour mass. The extension of the infiltration has been estimated to be at least 1 cm from the tumour borders for zone A and, in agreement with two other groups (8, 9), 3-4 cm considering the union of zones B and C.
The type of tissue nearby the tumour region also strongly affects the extent of tumour infiltration in addition to its cell density, due to the presence of natural barriers in the brain (e.g. ventricles and cerebral falx) (11) as well as to the differential motility of the tumour cells in white and gray matter (higher motility, and thus preferential tumour growth, is observed in white matter) (12, 13).
A study by Sandström et al. (6) has suggested that it might be useful to accompany a physical assessment of the quality of the GKRS treatment, which is usually performed by calculating a series of indexes, such as the Conformity Index (CI), the Selectivity Index (SI), the Gradient Index (GI), or the Paddick Conformity Index (PCI) (14), with a radiobiological assessment of the treatment outcome. The radiobiological assessment could be performed by calculating the Tumour Control Probability (TCP) and considering the probability of tumour cells' presence outside the target region. A similar approach has also been proposed by Le et al. (15) in the context of fractionated radiotherapy, who developed a computational model of glioblastoma growth that was combined with an exponential cell survival model to describe the effect of radiotherapy.
The aim of this work was to develop a computational and mathematical model for glioma tumour cell invasion outside the target border, based on diagnostic imaging modalities currently available, and to investigate the impact of tumour cells' presence beyond the target boundaries on the treatment outcome of radiosurgery. Thus, a novel mathematical model able to predict the spatial density distribution of cancerous cells outside the tumour boundary is presented and implemented here on a virtual brain tumour for different planning scenarios and the treatment outcome is simulated.
Materials and Methods
Simulation of virtual tumour. A manually drawn high-grade glioma was positioned in a virtual brain simulator using a segmented brain MR image from the Internet Brain Segmentation Repository (IBSR) (16). The contour of the tumour bulk (T) is shown in Figure 1, together with the infiltration regions 1 and 2 calculated according to the mathematical model presented in the following subsection.
Mathematical model of infiltration. A mathematical model able to predict the cell infiltration beyond the borders of the tumour volume visible in CT images was developed and implemented in MatLab (MathWorks, version 2016a).
The model calculates the cell invasion in the voxels outside the tumour in a stepwise manner, by considering layers of voxels equidistant from the tumour border. Calculations are performed from the innermost to the outermost layer (Figure 2A). The stepwise model of infiltration mimics the natural behavior of diffusive infiltration (17), where the actual cell density of each voxel in the inner layers can actively affect the invasion to the outer layers. The invasion of the tumour cells in the voxels in each layer is calculated in a random order, so that the stochastic nature of the infiltration process is more closely reproduced (18, 19).
A 3×3×3 matrix of voxels containing invading cells (invading voxels) is centered around an invaded voxel in order to account for a potential infiltration from all possible directions. Naturally, the geometry of the tumour has a strong impact on the invasion pattern and the corresponding calculated density distribution of the infiltrated cells. Thus, voxels located at the same distance from the tumour border (i.e. within the same layer) might result in a non - homogenous density of infiltrated cells, even when the voxels under consideration correspond to the same type of tissue. As exemplified in Figure 2B, voxels A and B located at the same distance from the tumour border might be invaded differently because point A is located within the concave shape of the tumour and has a larger number of adjacent tumour (invading) voxels compared to point B.
A Total Infiltration Factor (TIF) is the cumulative result of cell invasions from all adjacent neighbours, where the contribution from each invading adjacent neighbor depends both on the density of cancerous cells in the invading voxel as well as the diffusion of cells with respect to their motility through the host tissue.
Contour of the virtual tumour (T) and of the corresponding infiltration regions 1 and 2 resulting from the mathematical model described in the following subsection.
The TIF value of an invaded voxel k is defined by equation 1:
A 3×3×3 voxel matrix, S, centered on voxel k is considered by the algorithm.
Ci is the concentration of cancerous cells in the voxel i adjacent to the voxel under consideration k, and it is defined as the ratio of the number of cancerous cells in the voxel i and the number of cancerous cells per voxel in the tumour volume (0≤Ci≤1). The cell density in the tumour mass is assigned to be 8000 cell/mm3, as provided by the detectability threshold in Ref. (7) and is assumed to be homogeneous in volume. The cell invasive flux is assumed to be directed from regions of high cell concentrations to regions of low cell concentrations, as reasonably expected from Fick's basic laws of flux flow.
Di is a weighting factor that accounts for the diffusion coefficient of the brain tissue in the voxel i relative to white matter. Di is equal to 1 for white matter (Dw=1), and to 0.1 for gray matter (Dg=1) (12), while it is assumed to be equal to 0 for other brain tissues, such as ventricles, skull, falx cerebri, which can be considered as natural barriers that prevent cell diffusion.
The normalization factor Nk is introduced to ensure that the TIF of the voxels outside the tumour boundaries is not higher than 1. In the first iteration of the algorithm, Nk is assumed to be constant and equal to the maximum number of cancerous neighbors for the voxels in the first layer outside the tumour border. For the following iterations, Nk is calculated separately for each invaded voxel, based on its surrounding voxels and the type of the host tissues, as:
where nw and ng are the numbers of neighboring voxels located in the white and gray matters, respectively, and Dw and Dg are the relative diffusion coefficients for the white and gray matter, respectively.
Illustration of a selected slice of the virtual high-grade glioma tumour. In (A), the layer-wise approach, adopted to perform the calculation of cell invasion outside the tumour border in the algorithm, is shown. At each step of the algorithm, the total infiltration factor (TIF) is calculated for each of the voxels within the layer under consideration. The tumour bulk is illustrated in white, the first and second layers encompassing the tumour are shown in black and dark grey, respectively. In (B), the voxels A and B are equidistant from the tumour border and are depicted together with their tumourous adjacent neighbors (indicated with circles) within the 3x3x3 calculation matrices. Assuming the same type of host tissue is surrounding the tumour, voxel A is expected to have a higher concentration of infiltrated cells given the higher number of invading cancerous neighbors.
Conformity index (CI), selectivity index (SI), gradient index (GI), volume receiving at least 12 Gy of dose (V12), maximum dose (Dmax), and tumour control probability (TCP) results for the different treatment plans under consideration.
The iterative algorithm that calculates the TIF values given by Equation 1 is repeated until a mean cell concentration, equal to 60% in the tumour region, is reached at 1 cm distance from the tumour border in the white matter, in agreement with the histopathological observations by Yamahara et al. (10).
Furthermore, the maximal extent of the tumour cell infiltration, as reported from the histopathological analyses by Silbergeld and Chicoine (9) to be about 4 cm away from the tumour border, is used as an additional constraint for the algorithm (i.e. the maximum cell density at 4 cm from the tumour border is set to 1 cell per voxel).
Consequently, three regions are identified in the cell density distribution, resulting from the mathematical model: i) the Gross Tumour Volume (GTV), with a cell density equal to 8000 cell/mm3, ii) zone 1, corresponding to the volume of tumour cells' densities in the interval 100-60% of the cell density in the tumour bulk (i.e. range equal to 7,999-4,800 cell/ mm3), and iii) zone 2, corresponding to the volume of tumour cells' densities in the interval 60-0% of the cell density in the tumour bulk (i.e. range equal to 4,799-1 cell/ mm3).
Treatment planning strategies. Several treatment plans that are different in prescription doses and target definition were created by an expert medical physicist using LeksellGammPlan® (version 10.1.1, Elekta Instrument AB, Sweden).
The first treatment plan was created by prescribing a dose of 20 Gy at 50% isodose to the GTV and accounting for the microscopic spread of the disease in the surrounding tissue, as advised by the clinical experience of the planner. This resulted in a spatial dose distribution that was not strictly conformal to the delineated GTV. In the second treatment plan, the same isodose level (i.e. 20 Gy at 50% isodose) was kept strictly conformal to the GTV. Another group of treatment plans (treatment plans 3 to 6) was generated by considering a target in zone 1 and prescribing a dose of 20 Gy to the 50%, 40%, 30%, and 20% isodose level. Finally, the last group of treatment plans (treatment plans 7 and 8) picked targets in zone 2 and provided dose prescriptions equal to 20 Gy at the 50% and 20% isodose level, respectively (Table I).
Density distribution of infiltrated cells outside the tumour mass. The grey scale shows the cell density ranging from 8,000 cells/mm3 (i.e. the detectability threshold of the tumour mass) to 0 cells/mm3 (darkest grey line in the figure).
Treatment plan quality assessment. The assessment of the quality of Gamma Knife™ Radiosurgery (GKRS) treatment plans was performed, from a physical point of view, by considering the standard conformity indexes available on the Leksell GammaPlan® treatment planning system. The Conformity Index (CI) and the Selectivity Index (SI) are complementary indexes that measure both the target coverage and the volume of irradiated normal tissues around the target, respectively. Specifically, CI is defined as the ratio between the Treated Target Volume (TTV), i.e. the portion of the contoured target that receives a dose equal to or greater than the prescribed dose, and GTV, while SI is defined as the ratio between TTV and the Prescribed Isodose Volume (PIV), i.e. the volume that receives a dose equal to or greater than the prescribed dose (14). The Gradient Index is a measure of the dose fall-off beyond the prescription isodose volume, and is defined as the ratio of the volume enclosed by half of the prescription isodose over the prescription isodose volume (20).
Illustration of the (A) dose distribution, (B) cell surviving fraction, (C) cell density distribution before treatment, and (D) cell density distribution following treatment for plan 1. The contour of the GTV, zones 1 (dark grey line) and 2 (bright grey line) are shown.
In conjunction with a physical assessment of the treatment plan's quality, a radiobiological evaluation of the treatment outcome was also performed.
The tumour cell survival fraction was calculated in each voxel following the linear-quadratic formalism (21), as follows:
where SFi and Di are the cell survival fraction and the absorbed dose in the voxel of interest i, respectively. Values 0.24 Gy−1 and 0.03 Gy−2 were used for the radiobiological parameters α and β for grade IV glioma cells, respectively, as reported by Malaise et al. (22).
The Poisson linear-quadratic model of tumour response to radiation (23) was chosen for the calculation of the Tumour Control Probability (TCP):
where n is the number of voxels in the considered volume, c0i is the density of tumour cells in the voxel under consideration before irradiation, Vi is the voxel volume, and Di is the absorbed dose in the voxel i.
Results
In Figure 3 isodensity lines of the cancerous cell distribution, as calculated by the proposed mathematical model of infiltration, are superimposed on the segmented MRI image of the brain. Nine representative slices out of the total 16 slices are shown in the figure. The tumour mass is visible in the figure in slices 38-40. The outermost isodensity line (dark grey contour) illustrates the maximal extent of infiltration, i.e. the outline of the cancerous region. As seen, the cancerous cells are distributed over the different brain tissues and they preferentially grow in the white matter compared to gray matter, while their infiltration is hindered from other brain tissues, such as the ventricles. The influence of the concave shape of the tumour resulting in a higher density of infiltrated tumour cells in the concavity of the tumour is reflected in the distance between the isodensity lines, which are spatially closer in the concavity of the tumour compared to the other areas.
Illustration of the (A) dose distribution, (B) cell surviving fraction, (C) cell density distribution before treatment, and (D) cell density distribution after treatment for the treatment plan 2. The contour of the GTV, zones 1 (dark grey line) and 2 (bright grey line) are shown.
Figures 4, 5, 6, 7, 8, 9, 10 and 11 show the: i) dose distribution, ii) the map of the cell surviving fraction, iii) the tumour cell density distribution before treatment, and iv) the tumour cell density distribution following treatment for the plans indicated in Table I. Results are shown in one representative slice of the volume under consideration.
Similarities are seen in the SF maps and the cell density distribution following treatment for plans 1 and 2, due to the rather similar dose distribution delivered to the target. The maximum density value of the survived cells is also of the same order of magnitude in plans 1 and 2. The group of plans with different dose prescription levels to zone 1 (plans 3 to 6) also show similar results in terms of SF maps, cell density distribution following treatment, and maximum density of surviving cells. A hot spot with a rather high density of surviving cells following treatment in zone 2 is seen in this group of treatment plans, because of the particular dose distribution planned by the user in this case, aiming at the coverage of zone 1. Similar results are also found in the plans with a dose prescription to zone 2 (plans 7 and 8), where the cancerous cells are almost completely eradicated.
Results of the physical and radiobiological assessment of the treatment outcome are shown in Table I. For each treatment plan, CI, SI, GI, the volume receiving at least 12 Gy of dose (V12) and the maximum dose (Dmax), are reported together with the TCP. With the exception of plan 1, where the selectivity of the plan was lowered on purpose to implicitly account for tumour cell infiltration beyond the contoured target region, and of plan 2, for which high selectivity was difficult to achieve because of the particular (horse-shoe) geometry of GTV, the quality of treatment is relatively similar for all the plans from a physical point of view (i.e. similar CI, SI, and GI values).
Illustration of the (A) dose distribution, (B) cell surviving fraction, (C) cell density distribution before treatment, and (D) cell density distribution after treatment for the treatment plan 3. The contour of the GTV, zones 1 (dark grey line) and 2 (bright grey line) are shown.
As predicted, V12 increased with size of the target volume and the maximum dose reflected the different dose prescriptions assumed in the planning process.
As expected from the rather high cell density values found following the completion of the treatment, the calculated TCP value for the treatment plans from 1 to 6 is zero. On the contrary, when accounting for the tumour cell distribution outside the target boundary, by prescribing either the 50% or the 20% isodose to the contour of zone 2 (plans 7 and 8), an almost complete eradication of the tumour cells was achieved, and so TCP values were close to 1.
Discussion
Even though the number of studies allowing a direct comparison between the obtained tumour cell infiltration pattern, as calculated by this novel algorithm, and other findings is rather limited, the cell density distribution presented in Figure 3 meets the expectations and agrees with the results of models of glioma growth proposed by Bondiau et al. (12) and Swanson et al. (24). The developed infiltration model described in this work is substantially detached from the other diffusion-reaction models of glioma growth (12, 24-26). While diffusion-reaction models provide information about the temporal evolution of the tumour, the approach proposed by this computational model is to predict the density distribution of infiltrated cells, solely based on the shape of the tumour mass delineated on the diagnostic image and based on the underlying brain tissue type. Thus, our model does not require time information as such, but only generic information regarding the infiltration of the specific histopathology, as has been previously provided by Yamahara et al. (10). In fact, the pattern of infiltration predicted by the proposed algorithm here could also serve as a more accurate initial condition for reaction-diffusion models in order to predict the growth pattern of high-grade gliomas.
The concave shape of the simulated tumour at an arbitrary location within the brain tissue might not be appreciated as a realistic one, given the fact that such concavity would most likely develop only in the presence of barriers or tissues with low motility in the vicinity of the cells where the tumour grew initially. However, the described tumour shape was intentionally chosen to illustrate the capability of the model to simulate the impact of invasion from different directions, a feature that would be particularly beneficial in case of multiple tumour bulks, such as in areas where the tumour masses are likely to be invaded from several directions.
Illustration of the (A) dose distribution, (B) cell surviving fraction, (C) cell density distribution before treatment, and (D) cell density distribution after treatment for the treatment plan 4. The contour of the GTV, zones 1 (dark grey line) and 2 (bright grey line) are shown.
With regards to the modelling of the diffusive infiltration flux of cancerous cells from high- to low-density regions, it should be noted that motility and diffusion are not limited to cancerous cells. In fact, we presume that normal cells are likely to diffuse from healthy areas towards the cancerous regions that have a lower concentration of normal cells. Nevertheless, when a higher proliferation rate of cancerous cells as compared to the normal cells is reasonably assumed, the diffusion of normal cells into the cancerous region may be neglected, as hypothesized in this study.
Results presented in Table I highlight the importance of accompanying a merely dosimetric assessment of the quality of the GKRS treatment plan through standard conformity indexes with radiobiological evaluation through TCP calculations in case of highly infiltrative tumours, such as high-grade gliomas. In fact, conformity indexes may fail in distinguishing between the plans with different treatment outcomes, as clearly seen when comparing results for treatment plans 1-6 having TCP equal to 0, and treatment plans 7-8 having TCP equal to almost 1. The low selectivity index scored for plan 1 in Table I (SI=0.37) indicates that a non-negligible portion of the normal tissues surrounding the target is irradiated by the experienced planners on purpose, as they are aware of the infiltration issues that accompany tumours of this type and stage. In fact, it is common practice in GKRS of high-grade gliomas in some clinics to implicitly account for the probability of the presence of infiltrated cells outside the delineated target, and deliberately decrease the selectivity of the plan (personal communication with Pierre Barsoum from Karolinska Hospital, Sweden). Given the current absence of standardization, a successful treatment outcome is strongly dependent on the planners, both with respect to their knowledge and experience, as well as their perception of the geometry and density distribution of infiltration. The use of mathematical models able to predict the actual geometry of the cancerous area seems therefore highly necessary for creating a standard practice in this type of treatment. Furthermore, the model of cell infiltration proposed in this study could be used to guide the contouring of high-grade glioma targets for GKRS, since the variability in contouring of high-grade gliomas appears to be very high (27-29). In their study, Sandström et al. (28) reported major and clinically significant differences in target delineation for an anaplastic astrocytoma case. This is especially troublesome in conformal techniques, such as in GKRS, where no margins are applied. By evaluating the extent of infiltration prior to contouring, this problem could be resolved to some degree.
Illustration of the (A) dose distribution, (B) cell surviving fraction, (C) cell density distribution before treatment, and (D) cell density distribution after treatment for the treatment plan 5. The contour of the GTV, zones 1 (dark grey line) and 2 (bright grey line) are shown.
Our method also provides information on the quality of the treatment plans in treating the cancerous areas from the radiobiological viewpoint. The remarkable differences between the outcomes of the plans suggest that, while changing the target definition may dramatically influence the number of survived cells, an increase in the maximum delivered dose level does not seem to impact the treatment outcome.
TCP calculation in conjunction with maps of predicted tumour cell density following treatment offer valuable information that can be used as a guide during the treatment planning session to increase the quality of the plan. While TCP gives an indication on how well the treatment plan covers the whole cancerous area (i.e. both GTV and regional infiltration), maps of tumour cell density following the treatment can be used to determine the regions where clonogenic cells should not have been eradicated at the end of the treatment, while they indicate areas where re-planning of the treatment might be necessary.
Illustration of the (A) dose distribution, (B) cell surviving fraction, (C) cell density distribution before treatment, and (D) cell density distribution after treatment for the treatment plan 6. The contour of the GTV, zones 1 (dark grey line) and 2 (bright grey line) are shown.
Nevertheless, the actual feasibility of delivering the treatment plan with a higher dose assigned to the periphery of the infiltrated zone may ensure the control of the tumour, but it should be carefully considered. In fact, such a plan might not be clinically acceptable to avoid prescribing a dose that exceeds the threshold of clinical acceptance in a large volume of the brain. V12, i.e. the volume of the brain that receives doses above 12 Gy, is one of the most important predictors of adverse radiation effects (ARE) following radiosurgery and is used as one of the criteria for clinically accepting the treatment plan. When V12 exceeds 8.5 cm3, the associated risk of ARE is greater than 10% and the delivery of such a plan in a single session is not recommended (30). As seen in Table I, V12 assumes values below 5 cm3 for the plans of group (I), about 40 cm3 for the plans of group (II) and about 150 cm3 for the plans of group (III). Thus, according to the mentioned constraints on V12, the treatment plans for target zones 1 or 2 result in being clinically unacceptable. In these cases, radiosurgery might not be an option and different treatment strategies, possibly, also different fractionation schemes might need to be envisaged. Alternatively, potential new ways of prescribing the dose of GKRS for infiltrative tumours could be investigated. One such example is to segment the target into different areas according to the predicted density of infiltrated cells and prescribe different dose levels, with the aim of increasing the TCP as much as possible while keeping the V12 below the threshold of clinical acceptance. This might offer a feasible solution to be pursued in future investigations.
From our results, it is evident that there are potential shortcomings from using conformity indexes alone for the evaluation of the treatment outcome. Radiobiological parameters, such as TCP and maps of survived cells following treatment that are able to reflect different and complementary aspects of the treatment outcome, are found to be crucial when assessing the quality of GKRS treatment plans in cases of infiltrative tumours.
Illustration of the (A) dose distribution, (B) cell surviving fraction, (C) cell density distribution before treatment, and (D) cell density distribution after treatment for the treatment plan 7. The contour of the GTV, zones 1 (dark grey line) and 2 (bright grey line) are shown.
Illustration of the (A) dose distribution, (B) cell surviving fraction, (C) cell density distribution before treatment, and (D) cell density distribution following treatment for the plan 8. The contour of the GTV, zones 1 (dark grey line) and 2 (bright grey line) are shown.
Acknowledgements
Financial support was obtained from the Cancer Research Funds of Radiumhemmet.
Footnotes
↵* These Authors contributed equally to this study.
Authors' Contributions
Research design: IT-D, ML; performing the computations: ZMK; development of the algorithm: ZMK, ML, IT-D, HS; treatment plans in Leksell GammaPlan®: PB; drafting of the manuscript: ML, ZMK. All authors read, critically revised and approved the final manuscript.
This article is freely accessible online.
Conflicts of Interest
The Authors report no conflicts of interest.
- Received January 31, 2019.
- Revision received March 12, 2019.
- Accepted March 13, 2019.
- Copyright© 2019, International Institute of Anticancer Research (Dr. George J. Delinasios), All rights reserved
