Dose surface maps of the heart can identify regions associated with worse survival for lung cancer patients treated with radiotherapy.

BACKGROUND AND
PURPOSE: For lung cancer patients treated with radiotherapy, radiation dose to the heart has been associated with overall survival, with volumetric dose statistics widely presented. However, critical cardiac structures are present on the hearts surface, where this approach may be sub-optimal. In this work we present a methodology for creating cardiac surface dose maps and identify regions where excess dose is associated with in worse overall survival.
MATERIAL AND METHODS: A modified cylindrical coordinate system was implemented to map the cardiac surface dose for lung cancer patients. Validation was performed by mapping the cardiac chambers for 55 patients, fitting a point spread function (PSF) to the blurred edge. To account for this uncertainty, dose maps were blurred by a 2D-Gaussian with width described by the PSF. Permutation testing identified regions where excess dose was associated with worse patient survival. The 99th percentile of the max t-value then defined a cardiac surface region to extract dose, from each patient, to be analysed in a multivariable cox-proportional hazards survival model.
RESULTS: Cardiac surface maps were created for 648 lung cancer patients. Cardiac surface dose maps were blurred with a 2D- Gaussian filter of size σφ = 4.3° and σy = 1.3units to account for mapping uncertainties. Permutation testing identified significant differences across the surface of the right atria, p < 0.001, at all timepoints. The median dose to the region defined by the 99th percentile of the maximum t-value was 18.5 Gy. Multivariable analysis showed the dose to this region was significantly associated with survival, hazard ratio 1.01 Gy-1, p = 0.03, controlling for confounding variables.
CONCLUSIONS: Cardiac surface mapping was successfully implemented and identified a region where excess dose was associated with worse patient survival. This region extended over the right atria, potentially suggesting an interaction with the hearts electrical conduction system.

Introduction

In recent years, radiation induced cardiac toxicity (e.g. pericardial effusion, acute coronary syndrome, pericarditis, arrhythmia, and myocardial infraction) for lung cancer patients has been associated with poorer patient survival. In 2015, the results from RTOG 0617 showed radiation dose to the heart, in lung cancer patients, was associated with worse overall survival [1]. Since the publication of RTOG 0617 a number of papers have further analysed the impact of cardiac dose, with the majority of these papers analysing retrospective datasets from a single institution or data from clinical trials [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13]. These studies were summarised in the 2019 review by Zhang et al. [14], which concluded that no consistent cardiac dosimetric parameters yet exist. These papers analyse the dosimetric impact to either the whole heart or to individual cardiac sub-structures. The latter is arguably more appropriate as different cardiac sub-structures are likely to display different dose response relationships, with sub-structures in the base of the heart most commonly found [2], [3].

Additionally, the literature [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13] investigates dose to cardiac volumes, where some sub-structures are located primarily on the surface of the heart. The coronary arteries originate from the ascending aorta and are located across the surface of the heart. The electrical conduction system originates with the sinoatrial node which is located in the myocardium superior in the right atrium and the atrioventricular node located at the junction of the cardiac atria and ventricles. The electrical system then emerges from the apex of the heart and follows the surface. Additionally, the myocardium, the muscle of the heart, could be damaged due to the radiation exposure.

Dose surface maps have been implemented for the bladder and rectum in previous studies [15], [16], [17]. These studies showed the utility of these methods in identifying sub-regions of organs where excess dose was associated with toxicity. Dose surface maps have not yet been created for the heart. Therefore, we propose for the first time in this study, to analyse whether radiation dose on the surface of the heart has an impact on a patient’s overall survival. In this work, we propose a methodology to sample the dose across the hearts surface. We then use image-based data mining techniques to identify regions on the heart surface associated with overall survival.

Materials and methods

A total of 648 non-small cell lung cancer (NSCLC) patients were extracted from the treatment planning archive from a single institution. Institutional approval had been granted to use this data (research ethics committee reference: 17/NW/0060). All patients were stage 3 lung cancer patients, treated with 55 Gy in 20 fractions between 2011 and 2012. No further exclusion criteria were applied in patient selection, ensuring the heterogeneity of the patient population was maintained. Treatment plans were created in Pinnacle (Philips Radiation Oncology Systems, Andover, MA) using 6MV beams with 3D-conformal or IMRT delivery techniques and accounting for tissue inhomogeneity in the dose calculation. The number of treatment beams were variable across patients dependent on tumour volume and complexity of the target volume shape, beam angles were optimised to minimise dose to the contralateral lung.

ADMIRE vr3.0 (Elekta AB, Stockholm, Sweden) was used to propagate heart contours using a set of ten atlases contoured by a clinical oncologist. Atlases consisted of previously treated lung cancer patients where a heart contour had been delineated on the planning CT scan. Heart contours were defined according to the UK stereotactic body radiotherapy (SABR) consortium guidelines [18]. Propagated heart contours from the 10 atlases were combined using the STAPLE algorithm [19] for each patient and visually checked to ensure they were suitable for the creation of cardiac surface dose maps.

A cylindrical coordinate system was devised with the origin positioned at the centre of mass for each slice in turn. This approach accounted for the asymmetric nature of the heart and can be described as a modified cylindrical coordinate system. Alternatively, this could be viewed as a cylindrical coordinate system where the heart has been transformed such that centre of the heart contour in each slice coincides with the polar axis. For each slice, 0° was defined in the anterior direction and the azimuthal angle defined from that point (Supplementary material, Fig. S1).

The dose on the surface of the heart was sampled by using an image processing technique called splatting [20]. Here, the desired resolution of the surface map was first specified and a line from the origin was traced for each angle until the surface of the heart was intersected (defined by creating a binary mask for the heart contour). The dose at this point was then sampled, a pictorial representation is included in the Supplementary material, Fig. S2. Step sizes were set to be 1° producing a surface map with dimensions 360 by number of slices containing the heart contour in the superior-inferior direction (z-direction). The z-direction was resampled using linear interpolation to be 50 units to ensure the surface maps across all patients were spatially normalised (50 being the mean number of slices containing the heart contour across patients). This can be considered as a form of non-rigid registration, ensuring the heart surface anatomy was matched into the same frame of reference across all patients.

To validate the accuracy of the surface mapping, the cardiac chambers, the left and right atria and ventricles, were segmented for 55 patients (atlas-based propagation, using 10 atlases, using ADMIRE vr3.0). Patients were selected to account for the range of heart surface doses, tumour volume and tumour location (left and right tumours, upper and lower lobes). Contours were visually inspected and amended by a clinical oncologist. These segmentations allowed the anatomical localisation of any dosimetric effects observed in the analysis. Each chamber was first expanded by 1 cm to ensure it intersected the heart surface. The same method as described was applied to map the intersections of each chamber and the heart surface (process described pictorially in the Supplementary materials, Fig. S3). The maps from the 55 patients were then summed, producing a distribution of chamber positions.

The image of the distribution of chamber positions could be interpreted as a map of the average chamber positions, convolved with a 2D Gaussian distribution with widths corresponding to the error in each direction. Therefore, the errors in each direction could be estimated by deconvolving the distribution using a semi-blind Richardson-Lucy algorithm described in [21] which estimates the point spread function (PSF) of a blurred image (Supplementary material S4 includes the mathematical formula). This iterative procedure first estimated the PSF using a number of iterations of the Richardson-Lucy algorithm. This PSF was then used in a second series of iterations to deconvolve the image. The images were represented by 2D discrete functions which took the value of the pixel intensity at each coordinate in the image. The convolutions were performed using fast Fourier transform methods from the SciPy library (Python vr3) and the PSF was fit to a Gaussian by first estimating the parameters by calculating the moments of the distribution and then optimising via a least squares routine. The PSF fit for the z-direction and azimuthal angle defined a 2D-Gaussian function and was used to blur the surface dose maps to account for mapping errors.

Surface maps were grouped to calculate the mean dose distributions of patients who survived or did not survive at a given timepoint post-treatment (6, 12, 18 and 24 months). Patients were censored for follow-up. To test the significance in the dose differences between the two groups permutation testing was applied. This procedure was first described by Chen et al. [22], with 1,000 permutations performed. Briefly, the test statistic, the maximum t-value, was calculated from the difference in the mean dose in each voxel between the two groups, divided by the standard deviation of the voxels. Permutations generated random samples, determining the distribution of the maximum t-values, testing the null hypothesis that there was no difference in the groups. To ensure permutation testing was not influenced by the standard deviation the map of standard deviation across the surface map was plotted.

The region of highest statistical significance was defined as voxels with a t-value greater that the 99th percentile of the max t-value. The dose to this region, across all patients, was then extracted for further analysis. The surface dose was included in univariable and multivariable cox-proportional hazards models, calculating hazard ratios for all variables, accounting for other confounding variables; tumour volume, age, gender, T-stage, N-stage, mean lung dose, performance status (defined using the Eastern Cooperative Oncology Group (ECOG) performance status definition: 0 fully active to 4 completely bed bound). Variables were selected into the multivariable model using forward selection where they displayed a p-value <0.05 on univariable analysis. Kaplan-Meier curves, for overall survival, were plotted for patients receiving greater or less than the median dose to this region. Log-rank was calculated to assess any significance difference in overall survival. All statistical analysis was performed in R version 3.6 [23].

Results

Patient demographics for the 648 patients are included in Table 1. Fig. 1a shows the count in each bin from the surface mapping of the cardiac chambers, demonstrating the agreement of the delineated chambers. The PSF was fitted, to the distribution of chamber positions, using the semi-blind Richardson-Lucy algorithm, with a width of 4.3° in the angular coordinate and 1.3 unit in the z-direction.

Table 1

Patient demographics for the 648 patients included in the analysis are described.

Variable	Sub-variable	Sub-total	Percentage complete
Gender	Male	360	100%
	Female	288
Age (years) (median)		73 (39–95)	100%
Tumour volume (cm3) (median)		42 (20 – 79)	100%
T Stage	T1	76	94%
	T2	244
	T3	170
	T4	117
N stage	N0	290	94%
	N1	87
	N2	195
	N3	40
Performance Status	0	85	93%
	1	272
	2	202
	3	42

Fig. 1

Results from mapping the cardiac chambers for 55 patients. (a) Shows the distribution in position for each chamber, showing the number of patients in each position. (b) Defines the median chamber position with the dashed line defining one standard deviation of the distribution.

Permutation testing on the blurred dose maps showed significant results (p < 0.001) for all timepoints tested. Significance maps from 12 months are included as Fig. 2 (patients grouped as those who survived or did not survive at 12 months). After censoring patients for follow-up, the number of patients who survived and did not survive were most balanced at 12 months, providing the most robust results in the permutation testing. However, the map was representative of other timepoints tested; patients grouped on survival at 6, 12, 18, 24 months. Fig. 2a shows the significance map with the 99th and 95th percentiles of the maximum T-value highlighted. Fig. 2b shows the significance map with the chamber positions overlaid to allow regions of significance to be localised to the hearts surface anatomy. The region of significance was located across the base of the heart, with the highest significance extending between the right atrium and right ventricle. Interestingly, there was no significance across the left atrium or left ventricle. Importantly, Fig. 2c shows the standard deviation in dose across all patients. Note that the standard deviation was fairly homogenous across the majority of the 99th percentile region and displayed a different pattern than the t-map

Fig. 2

(a) Shows the results of the permutation testing with the 95th and 99th percentiles highlighted. (b) Shows the map overlaid with the cardiac chambers to aid localisation of the dosimetric effect. (c) Shows a FECOG of the standard deviation in dose across the surface map for all patients. Colour bars are included for each figure to describe the value of the T statistic in (a) and (b) (higher values show increased significance) and for the standard deviation in (c).

Mean dose to the region defined at the 99th percentile was extracted and used in a time-to-event survival analysis including other available clinical variables. Median surface dose across all patients was 18.5 Gy (range 0 – 55 Gy), compared to the mean whole heart dose of 12.7 Gy. Table 2 shows the univariable and multivariable models. On multivariable analysis the mean dose to the heart surface region was significantly associated with survival (Hazard Ratio (HR) 1.01, p = 0.03 as continuous variable). Tumour volume was also significant (p < 0.001, continuous), age (p = 0.01, continuous) and nodal stage (p < 0.05, categorical, reference N0). Mean lung dose was not significantly associated with survival in this model (p = 0.08). Interactions between these variables are included in the Supplementary material (S5).

Table 2

Univariable and multivariable analysis showing significance of the dose to the defined surface region of the heart and controlling for other tumour and clinical covariates. Performance status was not significant for this group of patients on univariable analysis and was not brought forward into the multivariable model.

	Univariate			Multivariate
	HR	(95% CI)	p	HR	(95% CI)	p
Tumour volume (log) (continuous)	1.42	(1.31–1.55)	<0.001	1.35	(1.21–1.50)	<0.001
Heart surface region mean dose (continuous)	1.02	(1.01–1.02)	<0.001	1.01	(1.00–1.03)	0.03
Age (continuous)	1.01	(1.00–1.02)	0.02	1.02	(1.00–1.03)	0.01
Gender (male vs. female)	1.35	(1.15–1.59)	<0.001	1.17	(0.96–1.43)	0.13
Mean lung dose (continuous)	1.07	(1.04–1.09)	<0.001	0.95	(0.90–1.01)	0.08
T-stage (T1 ref)
T2	1.56	(1.17–2.07)	<0.001	1.05	(0.74–1.49)	0.78
T3	1.98	(1.47–2.66)	<0.001	1.01	(0.69–1.48)	0.96
T4	2.35	(1.72–3.21)	<0.001	1.14	(0.77–1.71)	0.51
N-stage (N0 ref)
N1	0.87	(0.68–1.13)	0.3	0.73	(0.55–0.98)	0.04
N2	1.47	(1.22–1.78)	<0.001	1.25	(0.97–1.60)	0.09
N3	1.60	(1.13–2.26)	<0.001	1.23	(0.79–1.91)	0.36
Performance status (PS0 ref)
PS1	1.23	(0.93–1.62)	0.15	–	–	–
PS2	1.27	(0.96–1.70)	0.10	–	–	–
PS3	1.13	(0.80–1.67)	0.53	–	–	–
PS4	0.78	(0.19–3.19)	0.73	–	–	–

Kaplan Meier curves were plotted with patients’ grouped on those receiving greater than or less than the median dose, 18.5 Gy, to the identified surface region, Fig. 3, with log-rank p < 0.001. Median survival times for patients receiving lower than 18.5 Gy was 22 months (95% confidence interval 20 – 26 months) compared to patients receiving greater than 18.5 Gy median survival 13 months (95% confidence interval 12 – 15 months).

Fig. 3

Kaplan Meier curves showing the survival of patients who received greater than 18.5 Gy to the identified region on the surface of the heart versus those that received <18.5 Gy. Median survival times for the two groups were 12.9 months and 22.4 months respectively.

Discussion

In this work we have devised a method to map the radiotherapy dose on the surface of the heart for the first time. We applied our novel methodology to mapping the cardiac surface dose for 648 stage III NSCLC patients. In analysing the surface dose maps we identified regions in the base of the heart, overlapping with the right atrium, which were significantly associated with patient’s overall survival (p = 0.03). Importantly, apart from being treated with 55 Gy in 20 fractions, no exclusion criteria were applied in the patient selection. Therefore, we have maintained the heterogeneity of the ‘real world’ patient population with the majority of patients over 70 years old and performance status 2 and 3.

Dose surface maps have previously been created for the bladder and rectum to perform toxicity assessments for radiotherapy patients. Due to their shape the bladder lends itself to a spherical coordinate system while the rectum to a cylindrical coordinate system. Palorini et al performed a pixel-wise analysis of bladder dose surface maps to investigate localised effects resulting in toxicities for prostate cancer patients [15]. Significant differences in the bladder surface maps were found, with different spatial patterns for different toxicities, including urinary frequency being associated with higher dose on the trigone [16]. Wortel et al employed a cylindrical coordinate system for creating surface dose maps of the rectum for prostate cancer patients [17]. Their approach is similar to our modified cylindrical coordinate system implemented in this work. Average surface dose maps for reported rectal toxicities were created and permutation testing used to identify significant differences. For all toxicities reported significant differences in the surface dose maps were identified.

Dose surface maps have never been created for the heart before, for any patient group. The heart is cone-shaped, with the base of the heart broader than the apex. However, on an individual slice, the heart displays a circular nature. Hence our decision to implement a modified cylindrical coordinate system where each slice is handled independently of its neighbours allowing the radius to modified. This is a simplistic choice, however, the approach worked well as evidenced by the agreement of the position of the atria and ventricles with uncertainties of 4.3° in the azimuth direction and 1.3 units in the z-direction. Additionally, the simplicity of this approach makes independent validation easier.

To quantify the potential uncertainty in the process we fitted a PSF to the distribution of the mapped chamber positions. Identification of cardiac anatomy on CT is difficult, chambers were automatically contoured by ADMIRE with an independent observer performing a manual check. Certainly, the uncertainty measured will be an over-estimate capturing both the uncertainty in the mapping process and any uncertainty in the contour propagation. However, despite this over-estimation regions of significance were identified in the analysis, with the region defined by the 99th percentile larger than the measured uncertainty. Performing the statistical analysis without blurring the surface dose maps identified a similar region showing that these uncertainties have minimal effect on the results.

There is further potential for uncertainty introduced by respiration or cardiac motion. This will primarily act to blur the boundary of the heart introducing an uncertainty in the dose sampled. However, as the dose across neighbouring pixels is highly correlated this is not believed that these uncertainties will affect our results. Additionally, there is the potential for the heart to change shape during radiotherapy, particularly where chemotherapy is being delivered concurrently (due to the extra fluid delivered during this procedure). Patients in this study received induction chemotherapy minimising this possibility. These uncertainties can be better estimated by utilising 4DCT scans and incorporating on-treatment imaging, such imaging will be explored in further studies.

As summarised in the review article by Zhang et al. [14] no consensus yet exists as to the optimal cardiac region and dose threshold for radiotherapy planning. The literature analysing whole heart volumetric dosimetric parameters. In RTOG 0617 [1], patients were randomised to radiotherapy with 74 Gy versus 60 Gy. The higher dose arm showing worse overall survival, 20.3 months versus 28.7 months. The multivariable analysis showed that higher dose to 5% of the heart volume and higher dose to 30% of the heart volume were associated with increased mortality. Work exploring cardiac sub-structures has predominately identified the base of the heart as most dose sensitive. McWilliam et al. found the base of the heart as most dose sensitive in 1101 NSCLC treated with 55 Gy in 20 fractions, hazard ratio 1.25 for patients receiving higher than median dose to the heart (16.3 Gy) [2]. Similarly, in 803 patients treated with SABR the left atrium and the superior vena cava were associated with non-cancer patient death [3]. Johnson-Hart et al. analysed residual shifts from image guided radiotherapy of NSCLC patients [24]. This latter work showed that in patients with a 1 mm or greater residual shift, which pointed towards the base of the heart, showed worse overall survival. The results in this study are in agreement, highlighting a surface region located in the base of the heart as most dose sensitive. Importantly, in this patient cohort, whole heart dosimetric parameters were not significantly associated with overall survival on univariable analysis (Mean heart dose, p = 0.2, heart V30, p = 0.4 and heart V5, p = 0.5).

With any analysis exploring the association of dosimetric parameters with overall survival there exists issues of co-linearity. The test of correlations between the cardiac surface dose and clinical variables included in the multivariable model are included in the Supplementary material (S5). Importantly there is no correlation between tumour volume and heart surface dose, this is not unexpected as tumours can arise in either lung. The heart surface region identified exists primarily on the right-hand side of the heart, therefore, for a tumour of a given volume located in the left lung would contribute less dose to this region than if located in the right lung. Interestingly, the T-stage for the patients shows an association with the heart surface dose. N-stage shows a strong association with cardiac surface dose, higher N-stage patients will have more extensive disease in the mediastinum and therefore higher dose delivered across the heart. Lung mean dose showed a correlation of 0.7 with the cardiac surface dose. There is suggestion of an interplay effects between lung and heart dose and associated toxicities [25]. Such interactions need further investigation to determine the relevant importance of dose to the lungs versus dose to the heart. No other clinical variables showed a strong correlation with the heart surface dose. Despite these correlations, heart surface dose remained significant when these variables were controlled for in multivariable analysis.

The sinoatrial (SA) node is located in the wall of the myocardium, superior in the right atrium, located lateral to the superior vena cava. This region is potentially encompassed with our identified region of significance. This suggests that the observed poorer survival for patients receiving a higher dose to this region of the surface of the heart may be caused by damage to the SA node. Indeed, work by Vivekanandan et al identified changes in patient’s electrocardiogram (ECG) measurements pre- and post-radiotherapy [6]. However, they found that dose to the left atrial wall was most strongly associate with a measured change in a patient’s ECG, found in 38% of patients analysed. Patients were treated in a dose escalated, isotoxic clinical trial, with doses thresholds, on the left atrial wall, of 63–73 Gy identified. A different patient population to those included in this work who received 55 Gy, who represent ‘real world’ patients, with complex multi-morbidities and polypharmacy.

The patient database lacks information on multi-morbidities and cause of death, a common problem with retrospective datasets. Ideally, prospective studies are required to better capture patient’s burden and severity of multi-morbidities to better understand the interaction between underlying cardiac health and radiation dose. To better build this understanding we are currently recruiting patients to a prospective cardiac biomarker study funded by Yorkshire Cancer Research (Research ethics committee: 18/NW/0706). In this study cardiac comorbidity data will be extensively collected alongside prospective ultrasound, CT angiogram and circulating cardiac biomarkers. Such datasets will allow a more robust analysis to be performed, allow baseline cardiac health to be accounted for directly in the analysis. Additionally, cardiac dose surface maps may show utility in further patient groups where cardiac toxicities impact patient outcomes, for example breast cancer patients [26]. However, breast cancer treatments are likely to show less variability in dose than in lung cancer patients, potentially making the permutation testing less powerful.

In summary, we have devised a methodology for mapping the radiation dose across the surface of the heart for radiotherapy patients. Increased radiation dose to the surface of the heart, particularly across the right atrium, is significantly associated with poorer patient’s survival. Results should be validated in external datasets and prospective studies with detailed cardiac comorbidity information.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.