510 patients undergoing radiotherapy for H&N cancers at BC Cancer between November 2004 and July 2015 were enrolled in this study. Patients were treated with either intensity modulated radiation therapy or volumetric modulated arc therapy. Radiation dose for all radiotherapy plans were calculated using the analytical anisotropic algorithm using the same planning system, dose prescribing convention, and dosimetric grid size. DVHs for each patient’s parotid glands (both ipsilateral and contralateral to the tumor site) were extracted using DICOMautomation. Patients were excluded if: they were unable to follow written saliva collection procedures; they received atypical chemotherapy agents (i.e. an agent other than cetuximab, cisplatin, carboplatin, or gemcitabine); they received electron therapy; or they had previous interfering radiotherapy. In addition to routine clinical quality assurance procedures prior to delivery of radiotherapy plans, a single senior H&N Radiation Oncologist (JW) validated the consistency and accuracy of salivary contours of the parotid glands specifically for research quality assurance purposes after plan delivery. Although salivary function can also be impacted by radiation to the other major salivary glands, stimulated salivary function is most impacted by the parotid glands. In addition, the LKB model requires gland specific parameters and therefore cannot be fit or otherwise compared to the candidate approaches in a mixed-gland framework.

Stimulated saliva was measured prior to radiotherapy and one year following completion of treatment. Measurements comprise whole-mouth saliva collected with patients in prone or upright position over a five-minute period while chewing flavorless wax.

All clinical plans were created according to institutional guidelines using the Varian Eclipse treatment planning system. Dose-volume histograms for the parotid glands were extracted from clinical plans using DICOMautomaton [13], an open source toolkit for radiotherapy analysis.

The LKB model addresses the multidimensionality of a DVH as a predictor variable by reducing it to a single dose and volume in a process originally described by Lyman [6, 14]. This reduction of the multidimensional DVH to a single dose over an effective volume is justified by an assumed power law relationship, where \(i\) represents each step of the DVH.

$$V=\sum_\Delta __}_}\right]}^\frac$$

Transformed single-step histograms are assumed to have the same complication probability as the original one. The newly transformed dose and volume are then used to normalize the dose with \(m\times _(v)\) serving as an estimate of the standard deviation of dose, where V and D are the transformed values from the DVH and Vtotal is the total organ volume (in this case the parotid glands).

where \(_(v)\) is attained from another assumed power law.

$$_\left(1\right)=_(v)\times ^$$

The final model estimate is then calculated by plugging the estimated t into the cumulative distribution of a standard normal random variable.

Due to the LKB model requiring a binary definition of complication, patients were defined to have a complication if their post-radiation salivary flow rate was reduced to less than 25% of the preoperative rate (i.e., a ‘severe’ reduction). The model was fit using maximum likelihood. After transforming the DVHs to a single dose and volume, the ratio of post-treatment to pre-treatment whole salivary flow was dichotomized and the model fit resulting in the 3 organ specific parameters (TD50(1), m, n). Roesink et.al conducted a study of 93 patients in which these parameters were estimated to be 31 Gy, 0.54, and 1 respectively [7]. The fitted values as well as those reported by Roesink et al. were also used in the assessment of the candidate models fit with new methods.

The first alternative model addresses the high dimensionality of the DVH using a cubic spline basis. This procedure fits a polynomial function with a specified form to each of the DVHs in the dataset. For this application, it was decided that a spline function with 5 knots equally spaced across the range of observed doses imparted sufficient flexibility to adequately mimic the DVHs, Fig. 1. The resulting fitted model contains six coefficients which are then used as predictors in a logistic regression model, which also incorporates splines to improve model flexibility.

Three examples of spline approximation of dose-volume histograms and their approximation by cubic spline basis (red)

The second model extracts the volume recorded in the DVH at intervals of 1 Gy from 0 to 70 Gy. These values are then used directly as inputs into a neural network with a single hidden layer containing 12 nodes and a decay of 0.8, which is a form of regularization for the model. These model parameters were determined by using ten-fold cross validation to obtain optimal predictive performance.

Both of these candidate models contain tuning parameters, which were optimized using tenfold cross-validation of the AUC for predicting a decrease in salivary flow rate of 0.5*baseline. In the case of the regression-based approach, the cubic splines were applied to the model inputs with the number of knots being tuned by tenfold cross-validation. In the case of the neural network, the decay parameter was employed. The decay parameter regularizes the model penalizing the size of the weights to prevent overfitting and improves the performance of the model-fitting algorithm by reducing flat spots in the cost function by inducing a differential penalty between highly correlated inputs like those coming from the dose-volume histogram. Additionally, alternative architectures were tested in which the number of inputs were reduced to as few as 10-equally spaced readings from the DVH and hidden layer sizes ranging from 5 to 25. However, reducing the number of inputs did not improve performance with performance being negatively effected at the smallest number of inputs.

The data was partitioned into a training set containing 70% of the observations and a test set containing 30%. The predictive performances of the models were compared using area under the receiver operating characteristic curve in the test set. Sensitivity to the cutoff for reduction in salivary flow rate was examined by including a variety of other potential cutoffs. All analyses were conducted in the R statistical computing program [15].