High rates of local recurrence in tobacco-related head and neck squamous cell carcinoma (HNSCC) are commonly attributed to unresected fields of precancerous tissue. Because they are not easily detectable at the time of surgery without additional biopsies, there is a need for noninvasive methods to predict the extent and dynamics of these fields. Here, we developed a spatial stochastic model of tobacco-related HNSCC at the tissue level and calibrated the model using a Bayesian framework and population-level incidence data from the Surveillance, Epidemiology, and End Results (SEER) registry. Probabilistic model analyses were performed to predict the field geometry at time of diagnosis, and model predictions of age-specific recurrence risks were tested against outcome data from SEER. The calibrated models predicted a strong dependence of the local field size on age at diagnosis, with a doubling of the expected field diameter between ages at diagnosis of 50 and 90 years, respectively. Similarly, the probability of harboring multiple, clonally unrelated fields at the time of diagnosis was found to increase substantially with patient age. On the basis of these findings, we hypothesized a higher recurrence risk in older than in younger patients when treated by surgery alone; we successfully tested this hypothesis using age-stratified outcome data. Further clinical studies are needed to validate the model predictions in a patient-specific setting. This work highlights the importance of spatial structure in models of epithelial carcinogenesis and suggests that patient age at diagnosis may be a critical predictor of the size and multiplicity of precancerous lesions. Cancer Res; 76(24); 7078–88. ©2016 AACR.

Major Findings

Patient age at diagnosis was found to be a critical predictor of the size and multiplicity of precancerous lesions. This finding challenges the current one-size-fits-all approach to surgical excision margins.

Head and neck squamous cell carcinoma (HNSCC) arises in the epithelial lining of the oral cavity, pharynx, and larynx. The annual incidence rate of HNSCC is estimated to be around 600,000 new cases worldwide (1), and in the United States alone, the death toll is approximately 11,500 cases per year (2). While a subgroup of HNSCC, including the oropharynx, is caused by infection with high-risk types of the human papillomavirus (HPV; ref. 3), the majority of HNSCC are HPV-negative and primarily associated with tobacco use and alcohol consumption (1). Despite a growing number of therapeutic strategies, survival in HPV-negative HNSCC has not improved significantly over the past decades, with a low median survival of approximately 20 months (4).

Poor prognosis in tobacco-related head and neck cancers is commonly attributed to the development of local recurrences and metastases after removal of the primary tumor (1). Field cancerization, or the presence of premalignant fields surrounding the primary tumor, has been shown to drive the high rate of local recurrence (5–11). In fact, molecular studies have shown that the majority of HPV-negative HNSCC develop within local fields of premalignant cells that are clonally related to the resected primary cancer (7, 9). These fields can be much larger than the actual carcinoma and are generally difficult to detect without genomic analyses due to their visually normal appearance (12). If such an invisible premalignant field extends beyond the surgical margins, the portion of the field left behind after resection of the primary tumor increases the risk of subsequent recurrence and contributes to poor prognosis (11). The occurrence of precancerous fields was first reported by Slaughter and colleagues (5) and has since been documented in most epithelial cancers (10, 13–15).

To date, it remains difficult to account for the phenomenon of field cancerization in clinical practice. The main reason for this translational barrier is a poor understanding of the dynamics and geometry of these invisible fields. Indeed, it is impossible to account for the risk factor of an unresected field of precancerous tissue in absence of any information on the extent of the suspected field. To overcome this barrier, we synthesized data and knowledge sources from the tissue, clinical and population scales to develop a quantitative model of HNSCC carcinogenesis that accounts for spatial features of the precancer field. On the basis of this model, we then sought to identify aspects of standard clinical practice that could be improved by means of patient-specific modeling tools.

Quick Guide to Equations and Assumptions
Model Assumptions

We developed and calibrated a spatial stochastic model of tobacco-related head and neck squamous cell carcinoma at the tissue level. The major model assumptions are:

• (Epi)genetic events lead to mutated cells with increased fitness advantage, and mutant clones can spread through the basal layer of the affected epithelium prior to onset of invasive cancer.

• Because of the time scales of carcinogenesis, only long-lived progenitor cells in the basal layer of the epithelium are relevant; this monolayer of progenitor cells is modeled as a two-dimensional lattice, where each node is occupied by a cell.

• Wild-type and mutant progenitor cells undergo evolutionary competition in the basal layer; only nearest-neighbor interactions are considered.

• Tissue architecture and homeostasis are maintained (constant population size) until the invasive cancer stage.

Key Model Parameters

Total population size (⁠$N$⁠); cellular transition rates from normal to precancer (⁠${u_1}$⁠) and from precancer to carcinoma in situ (⁠${u_2}$⁠), respectively; relative proliferative advantage of precancer (⁠${s_1}$⁠) and carcinoma in situ cells (⁠${s_2}$⁠); mean sojourn time from preclinical lesion to clinical diagnosis with cancer (⁠$1/\psi$⁠).

Key Equations

We were interested in quantifying the geometry of the field of precancerous cells surrounding the tumor at time of diagnosis (⁠$\sigma _3}$⁠). The key equations are:

• (i) The survival function (probability that cancer has not been diagnosed by time (⁠$t$⁠) of the model,

### Incidence data

To parameterize the tissue-level model and ensure its compatibility with population-level data, we calibrated the mesoscopic approximation of the microscopic model on the basis of age-specific incidence rates from the Surveillance, Epidemiology, and End Results (SEER) program of the National Cancer Institute (18 Registries, 2000–2012) in a Bayesian framework. As we focused on HPV-negative cancers, we restricted our search to males and females diagnosed with malignant tumors in the following categories of “Site Recode ICD-O-3/WHO 2008”: lip, tongue, floor of mouth, gum and other mouth, hypopharynx, and larynx (31). We further excluded cases with primary site labeled as “C01.9-Base of tongue” because cancers at this site are commonly HPV-positive (32). We derived the number of susceptible individuals and the number of cancer cases diagnosed in each of the following age groups: 15–19 years, 20–24 years,…, 80–84 years, and 85+ years. We note that the age-specific incidence curve for HPV-negative HNSCC has a peaked profile, see Fig. 2A. However, the model developed here predicts monotonically increasing incidence rates and thus cannot recapitulate peaked incidence patterns. Because changing baseline risks are more likely to cause the peaked profile than changes in biological parameters (see Supplementary Information S1), we truncated the incidence data used for parameter inference at the age of 74 years instead of introducing additional age-dependent model parameters. The vast majority of HPV-negative HNSCC are attributed to tobacco consumption (33, 34); thus, we reduced the pool of susceptible individuals to the approximately 20% of current smokers within the considered age groups (35). Furthermore, we assumed the exposure to start at the age of 15 years (36, 37).

Figure 2.

MCMC summary. MCMC Metropolis Hastings algorithm was used to infer the posterior distribution of the model parameters $\lambda$⁠, $\theta$⁠, and $\psi$⁠. A total of 500,000 iterations were computed; the first 50,000 burn-in iterations were discarded and then 1 in 100 samples was used to estimate the posterior distribution. A, Model fit to SEER data; age-specific incidence is restricted to smokers (see main text for details). The data point at 85 years reflects the incidence estimate for patients aged 85 years and older. B, Marginal distributions of $\lambda$⁠, $\theta$⁠, and $\psi$⁠.

Figure 2.

MCMC summary. MCMC Metropolis Hastings algorithm was used to infer the posterior distribution of the model parameters $\lambda$⁠, $\theta$⁠, and $\psi$⁠. A total of 500,000 iterations were computed; the first 50,000 burn-in iterations were discarded and then 1 in 100 samples was used to estimate the posterior distribution. A, Model fit to SEER data; age-specific incidence is restricted to smokers (see main text for details). The data point at 85 years reflects the incidence estimate for patients aged 85 years and older. B, Marginal distributions of $\lambda$⁠, $\theta$⁠, and $\psi$⁠.

Close modal

### Parameter inference

To compare the model predictions with SEER data, we first derived the age-specific incidence function under the evolutionary model at the tissue-level. To this end, we defined the random variable ${\sigma _3}$ as the time from mean age at smoking initiation to diagnosis with invasive cancer and calculated the survival function

where

$\gamma ( {a,x} ) = \mathop \int _0^x {t^{a - 1}}{e^{ - t}}dt$ is the lower incomplete ${\rm{\Gamma }}$-function, and $\psi$ is the transition rate from onset of CIS to diagnosis. The parameter groups (eq. C) have the following interpretation: $\lambda$ is the initiation rate of expanding type 1 clones, see also previous paragraph; $\theta$ can be rewritten as $( {\gamma t_2^{ - 3}} )/3$⁠, where ${t_2} = {( {c_2^2{u_2}{s_2}} )^{ - 1/3}}$ is approximately the time it takes for an expanding type 1 clone to give rise to an expanding type 2 clone (28). On the basis of $S( t ),$ we then computed the hazard function$h( t ) = - S^′( t )/S( t )$⁠, which provides a direct link to the SEER-derived age-specific incidence rates. The derivation of the likelihood function used for inference is outlined in Supplementary Information S1.

We then used a computational Markov Chain Monte Carlo (MCMC) Metropolis Hastings algorithm to derive posterior distributions for the parameters $\lambda$⁠, $\theta$⁠, and $\psi$⁠. We used improper priors over $[ {0,\infty } ]$ for all 3 parameters and ran a Markov chain of length 500,000. Computations were performed in MATLAB (Version 8.5.0, The MathWorks Inc. 2015). To determine plausible initial conditions and step lengths for the MCMC algorithm, we independently estimated the 3-parameter groups based on literature estimates of the underlying biologic parameters, see Supplementary Information S2. After discarding the first 50,000 burn-in steps, we sampled every 100th step to characterize the posterior distributions. Convergence to steady-state was visually ascertained. Credible intervals (CI; 95%) for the posterior distributions were computed by discarding the 2.5% largest and smallest values for each parameter. The marginals of the posterior distributions for $\lambda$⁠, $\theta ,$ and $\psi$ are summarized in Fig. 2B. Median estimates and 95% CIs of the posterior distributions, together with the literature-based order of magnitude estimates are summarized in Table 1. Goodness of fit was deemed satisfactory on the basis of visual inspection, see Fig. 2A.

Table 1.

Parameter estimates: MCMC and literature-based

ParameterUnitsPriorMedian (95% CI)Literature estimates
λ [y−1$[ {0,\infty } ]$ 4.3 × 10−3 (4.0 × 10−3–4.8 × 10−3∼10−3
θ [y−3${\rm{\ }}[ {0,\infty } ]$ 5.9 × 10−6 (4.4 × 10−6–7.9 × 10−6∼10−4
$\psi$ [y−1$[ {0,\infty } ]$ 9.8 × 10−2 (6.7 × 10−2–13.5 × 10−2∼10−1
ParameterUnitsPriorMedian (95% CI)Literature estimates
λ [y−1$[ {0,\infty } ]$ 4.3 × 10−3 (4.0 × 10−3–4.8 × 10−3∼10−3
θ [y−3${\rm{\ }}[ {0,\infty } ]$ 5.9 × 10−6 (4.4 × 10−6–7.9 × 10−6∼10−4
$\psi$ [y−1$[ {0,\infty } ]$ 9.8 × 10−2 (6.7 × 10−2–13.5 × 10−2∼10−1

NOTE: For each parameter, the median and 95% CIs for the MCMC Metropolis-Hastings–derived posterior distributions and the literature-based estimates are shown.

To compare model predictions with clinical recurrence patterns, SEER (13 Registries, 1991–2011) was queried for a retrospective analysis of the recurrence risk in patients diagnosed and treated for HNSCC lesions of size ≤ 3 cm. Data extraction was performed using an MP-SIR session with the software SEER*Stat (http://seer.cancer.gov/seerstat). Again, we focused on the predominantly HPV-negative sites as specified under “Model calibration.” Only patients aged 18 to 85 years at diagnosis with a known tumor size of ≤30 mm were included. Patients with unknown surgery status, unknown radiation status, and those identified at autopsy or on death certificate only were excluded. The following categorical variables were extracted from the dataset: age at diagnosis (younger: 18–49 years, older: 50–85 years) and treatment (surgery only, radiation only, surgery and radiation). For each patient, the time from diagnosis to recurrence (defined as second malignant event in any of the included head and neck sites), death, or censoring (defined as the minimum of end-of-study, loss to follow-up, and 10 years) was calculated. For each of the three treatment groups with known treatment status, a univariate Cox proportional hazards regression for recurrence was performed on the basis of the categorical age at diagnosis variable, stratified into younger (<50 years) and older (≥50 years) patients. Kaplan–Meier curves were used to visualize the results. Hazard ratios with 95% CIs were computed, and all P values were calculated as two-sided with significance declared for P values below 0.05. Computations were performed in R (http://www.R-project.org).

### Local field size at diagnosis

To characterize the field geometry at time of diagnosis with invasive cancer, we first derived the size distribution of the local field. More precisely, denoting the radius of the local field at time $t$ by ${R_l}( t )$⁠, we conditioned on $\{ {{\sigma _3} = t} \}$ to compute the probability density function of the local field radius as follows:

where $\zeta \equiv 3\theta /c_2^3$⁠. We refer to Supplementary Information S3 for details of the calculation. The corresponding density functions (eq. D) for ages at diagnosis of 40, 60, and 80 years are shown in Fig. 3A. Unfortunately, these results have an explicit dependence on the radial growth rate of the precancer field, ${c_2}( {{s_1}} )$⁠, which cannot be directly estimated due to identifiability constraints (we used literature-based estimates instead). To overcome this issue, we focused on the relative size of a precancerous field compared with the field size at the age of 50 years. More precisely, we introduced the relative mean field radius (RMFR), defined as the mean field radius at a given age divided by the mean field radius at age 50 years. Indeed, the RMFR is completely specified by the inferred parameters $\lambda$⁠, $\theta$ and $\psi$ and does not explicitly depend on ${c_2}$ (see Supplementary Information S3). On the basis of the posterior distributions from the Bayesian inference, we computed the RMFR as a function of patient age at diagnosis in Fig. 3B. The model predicted an approximately linear increase in RMFR between the ages of 20 and 90 years, and the median (95% CI) RFMR at age 90 was found to be 1.90 (1.88–1.93).

Figure 3.

Field characteristics as a function of age at diagnosis. A, Probability density function of the field radius at diagnosis is shown for patient ages 40, 60, and 80 years, respectively. Calculations were performed using the median posterior values of $\lambda$⁠, $\theta ,$ and $\psi$ as shown in Table 1, and ${c_2} = 0.4$ (see also Supplementary Information S2). B, On the basis of the posterior distributions, the statistics of the relative mean field radius (RMFR) for different ages at diagnosis were computed. Median and 95% CIs are shown. C, Probability distribution of the number of distant fields is shown for patient ages 40, 60, and 80 years, respectively. Parameter values as in A. D, Probability of multiple unrelated fields at time of diagnosis was computed on the basis of the posterior distributions. Median and 95% CIs are shown.

Figure 3.

Field characteristics as a function of age at diagnosis. A, Probability density function of the field radius at diagnosis is shown for patient ages 40, 60, and 80 years, respectively. Calculations were performed using the median posterior values of $\lambda$⁠, $\theta ,$ and $\psi$ as shown in Table 1, and ${c_2} = 0.4$ (see also Supplementary Information S2). B, On the basis of the posterior distributions, the statistics of the relative mean field radius (RMFR) for different ages at diagnosis were computed. Median and 95% CIs are shown. C, Probability distribution of the number of distant fields is shown for patient ages 40, 60, and 80 years, respectively. Parameter values as in A. D, Probability of multiple unrelated fields at time of diagnosis was computed on the basis of the posterior distributions. Median and 95% CIs are shown.

Close modal

### Multiple fields at diagnosis

In addition to the local field size at diagnosis, we estimated the probability of harboring distant fields in addition to the local field. The exact distribution of the total number $M( t )$ of clonally independent fields, including the local and all distant fields, is found in Supplementary Information S3. In Fig. 3C, the probability distribution for the number of multiple fields is shown for different ages at diagnosis (40, 60, and 80 years), based on the median posterior values of $\lambda$⁠, $\theta ,$ and $\psi$⁠. Using this distribution, we then computed the probability of harboring at least 2 clonally unrelated fields in the head and neck region at time of diagnosis,

Accounting for the posterior distributions of the parameter groups, the probability of harboring at least one distant field is shown as a function of age at diagnosis in Fig. 3D. The model predicts that the probability (95% CI) of harboring at least one distant field increases by one order of magnitude from 2.1% (2.0%–2.4%) at 20 years of age to 20.0% (18.6%–22.7%) at 80 years of age.

### Age-specific recurrence patterns

To test the validity of our model, we made model-based predictions and tested them against outcome data from SEER. The first prediction was based on the fact that current clinical practice recommends an age-independent excision margin width of ∼1 cm (38, 39). Considering the predicted increase of local field size with age at diagnosis (Fig. 3A), an age-independent margin width implies that, for the same tumor size, the area of precancerous tissue left behind after resection of a primary tumor is bigger in older patients (with larger precancer fields) than in younger patients (with smaller precancer fields), see Fig. 4. This increase in recurrence risk in older patients is further increased by the elevated probability of harboring distant fields (Fig. 3B) that may not be affected by local excision of the primary tumor. To test the resulting hypothesis of an age-dependent recurrence risk, we performed a retrospective time-to-event analysis for patients with small ≤ 3 cm) tumors who received surgery but no radiation (n = 9,665). In alignment with our model-based prediction, we found a 58% increase in recurrence risk in older patients (≥50 years, n = 7,774) than in younger patients (<50 years, n = 1,891; HR, 1.58; 95% CI, 1.27–2.10; P < 0.0001). The corresponding Kaplan–Meier curve for recurrence-free survival is shown in Fig. 5A.

Figure 4.

Surgical margins and residual field. A, Illustration of age-related differences in local field size and number of unrelated fields. Before surgery, only one local field may be present in a younger patient (left), whereas a larger local field and additional distant fields may be present in an older patient (right). B, During surgery, the local field is removed in the younger patient (left) but only partially resected in the older patient (right), where the residual field portions elevate the risk of recurrence.

Figure 4.

Surgical margins and residual field. A, Illustration of age-related differences in local field size and number of unrelated fields. Before surgery, only one local field may be present in a younger patient (left), whereas a larger local field and additional distant fields may be present in an older patient (right). B, During surgery, the local field is removed in the younger patient (left) but only partially resected in the older patient (right), where the residual field portions elevate the risk of recurrence.

Close modal
Figure 5.

Recurrence-free survival after treatment of primary tumor. Recurrence-free survival among younger (<50 years) and older (≥50 years) patients, stratified by treatment regime. A, surgery only; B, surgery and radiation; C, radiation only. Analyses restricted to primary tumors of size ≤ 3 cm. Significant difference in recurrence risk for surgery only, P < 0.0001.

Figure 5.

Recurrence-free survival after treatment of primary tumor. Recurrence-free survival among younger (<50 years) and older (≥50 years) patients, stratified by treatment regime. A, surgery only; B, surgery and radiation; C, radiation only. Analyses restricted to primary tumors of size ≤ 3 cm. Significant difference in recurrence risk for surgery only, P < 0.0001.

Close modal

Because adjuvant radiation therapy after surgery targets the margins of the resected tissue portion, the combination of surgery and radiation is more likely to eradicate the entire local field and hence diminish the recurrence risk from a residual field portion. Therefore, we hypothesized a reduction in the age-related difference in recurrence risk among patients who received radiation in addition to surgery. This prediction was corroborated by the analysis of outcome data. Repeating the above time-to-event analysis for patients who received both surgery and radiation, we found a smaller, nonsignificant difference (HR, 1.29; 95% CI, 0.91–1.84) between older (≥50 years, n = 4,537) and younger patients (<50 years, n = 993), see Fig. 5B.

Finally, we note that for the patient group with radiation but without surgery (Fig. 5C), no significant difference in recurrence risk between the two age groups was observed. Interestingly, in comparison to younger patients (<50 years, n = 266), older patients (≥50 years, n = 2,622) did have a slightly smaller recurrence risk (HR, 0.78; CI, 0.40–1.52). However, we were not able to draw any definite conclusions because of a small sample size and patient attrition in the younger age group.

Leftover fields of precancerous tissue that are not removed during resection of HPV-negative, tobacco-related HNSCC lead to an increased risk of local recurrence. Because of their visually normal appearance, epithelial precancer fields are generally not detectable at the time of surgery without harmful biopsies of the surrounding tissue portions. Thus, there is a need for noninvasive methods to predict the extent of the field and to make patient-tailored decisions with respect to optimal treatment modality, size of excision margins, and frequency of postoperative surveillance. In this work, we developed and calibrated a mechanistic model of spatial tumorigenesis in tobacco-related HNSCC. Following the lines of previous multiscale modeling work (40–43), we linked biologic mechanism at the tissue scale to data at the population scale and used a Bayesian framework to ensure compatibility of the model dynamics across scales. Using the model, we made predictions about the geometry of precancerous fields at the time of diagnosis with invasive cancer and tested them against outcome data.

The models predicted a strong dependence of the local field size on age at diagnosis. More precisely, we found the expected field size of the local precancerous field surrounding the primary tumor to double between the ages of 50 and 90 years. The model further predicted a substantial increase in the probability of harboring independent synchronous fields at the time of diagnosis, increasing to 20% at age 80. Together, these findings indicate that the current one-size-fits-all approach to the width of surgical excision margins may need to be critically re-evaluated.

Our study has several limitations. First, inherent limitations of the SEER database such as ascertainment biases and incomplete recurrence records may impact the validity of our results (46). Second, HPV status is not recorded in SEER, and despite careful selection of sites that are predominantly HPV-negative, a portion of recorded cases may have been misclassified. Third, although the incidence data used for model calibration were restricted to a relatively short period (2000–2012), it is likely subject to secular trends that are at least partially due to changing smoking patterns in the population (4, 47). In addition, it has been shown that smoking cessation leads to a slow decrease in head and neck cancer risk (48), which may result in differential field sizes between former and current tobacco users. In future work, the use of more granular smoking prevalence and cancer incidence data, adjusted for secular trends, is expected to address these issues and improve the model predictions. Fourth, a limitation shared with most multistage modeling analyses is the assumption of identical parameters for all individuals. This issue is partially mitigated by the Bayesian approach, which provides posterior distribution of parameters rather than point estimates. However, incorporating patient-level heterogeneity into the modeling framework constitutes a critical next step toward the long-term goal of developing personalized approaches to head and neck cancer care. Finally, we did not explicitly account for the role of the immune system as a first line of defense against neoplastic progression. Although immune effects could be incorporated into the model, current knowledge about the exact mechanisms of immune response to neoplastic transformation seems insufficient to develop meaningful models.

Historically, different types of multistage models have been used to infer the nature of cancer-causing mechanisms on the basis of incidence and mortality data. The first such models (49) were based on the assumption that cancer arises as the product of an organ-specific number of rare mutations. These models assumed a well-mixed population of cells and neglected cellular dynamics and spatial tissue structure. Later, these models were extended to account for clonal expansions of precancerous cells (40, 41, 50) and used to analyze the number and size of premalignant clones in nonspatial populations, both for exponential mean growth (51) and for more general growth dynamics (52, 53). In parallel to multistage models, population dynamic models such as the Wright–Fisher and Moran processes have also been used to model cancer initiation under well-mixed assumptions (54, 55). In these models, the expansion of premalignant clones is constrained by competition with healthy cells or premalignant cells at other stages. Our current work constitutes an extension of both multistage and population dynamic models. In particular, we accounted for the spatial structure of the epithelial lining of head and neck sites, and we developed a mechanistic model on the basis of the current understanding of tissue homeostasis and molecular biology of HPV-negative HNSCC. Using the spatially explicit model, we achieved a very good fit to the SEER-based incidence data, Fig. 2A. For the sake of comparison, we also fitted the classical multistage model (49) to the incidence data, see Supplementary Information S4. While the inferred number of 3 stages corresponds well to the current understanding of HNSCC etiology, the model fit was found to be poor. More recently, A study (47) showed that nonspatial 2-stage clonal expansion models yield good fits to the data when accounting for period and cohort effects in biologic parameters. Nevertheless, our findings suggest that for solid cancers subject to field cancerization, spatial effects may play an important role in shaping population-level incidence patterns.

While our predictions for age-dependent recurrence risks and field sizes were successfully corroborated by SEER-based outcome data, pathologic studies are needed to provide a definite confirmation of our models. Our predictions of an increase in field size in older patients will hopefully spur a critical re-evaluation of the one-size-fits-all approach to surgical excision margin width, as well as the effectiveness of surgery without radiation therapy in older patients with extended exposure and potentially very large fields. In addition, appropriate follow-up intervals for the monitoring of local recurrences may be optimized according to patient-specific predicted field geometries, and a quantitative understanding of the field dynamics combined with a set of reliable biomarkers for premalignant fields may enable physicians to perform targeted risk-assessments in high-risk groups, such as heavy smokers. Finally, the proposed modeling framework can be applied to other 2-dimensional epithelial sites affected by field cancerization, such as bladder, esophagus and skin, and it may provide valuable insights into observed differences in field extent and outcome between HPV-positive and HPV-negative head and neck cancers.

No potential conflicts of interest were disclosed.

The views expressed in this article are those of the author and do not necessarily represent the views of the Department of Veterans Affairs or the United States government.

Conception and design: M.D. Ryser, W.T. Lee, K.Z. Leder, J. Foo

Development of methodology: M.D. Ryser, W.T. Lee, K.Z. Leder, J. Foo

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): W.T. Lee

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M.D. Ryser, N.E. Ready, J. Foo

Writing, review, and/or revision of the manuscript: M.D. Ryser, W.T. Lee, N.E. Ready, K.Z. Leder, J. Foo

Study supervision: W.T. Lee, J. Foo

The authors thank Prof. F. Michor (Dana-Farber Cancer Institute) and Prof. R. Durrett (Duke University) for valuable feedback on model design and analysis.

The study was supported by the following grants: NIH R01-GM096190 and SNSF P300P-154583 (M.D. Ryser); NSF CMMI-1362236 (K.Z. Leder); and NSF DMS-1224362 and NSF DMS-1349724 (J. Foo).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Leemans
CR
,
Braakhuis
BJ
,
Brakenhoff
RH
.
The molecular biology of head and neck cancer
.
Nat Rev Cancer
2011
;
11
:
9
22
.
2.
Siegel
R
,
D
,
Jemal
A
.
Cancer Statistics, 2013
.
CA Cancer J Clin
2013
;
63
:
11
30
.
3.
Marur
S
,
D'Souza
G
,
Westra
WH
,
Forastiere
AA
.
Hpv-associated head and neck cancer: a virus-related cancer epidemic
.
Lancet Oncol
2010
;
11
:
781
9
.
4.
Chaturvedi
AK
,
Engels
EA
,
Pfeiffer
RM
,
Hernandez
BY
,
Xiao
W
,
Kim
E
, et al
Human papillomavirus and rising oropharyngeal cancer incidence in the United States
.
J Clin Oncol
2011
;
29
:
4294
301
.
5.
Slaughter
DP
,
Southwick
HW
,
Smejkal
W
.
Field cancerization in oral stratified squamous epithelium; clinical implications of multicentric origin
.
Cancer
1953
;
6
:
963
8
.
6.
Califano
J
,
van der Riet
P
,
Westra
W
,
Nawroz
H
,
Clayman
G
,
S
, et al
Genetic progression model for head and neck cancer: implications for field cancerization
.
Cancer Res
1996
;
56
:
2488
92
.
7.
Tabor
MP
,
Brakenhoff
RH
,
van Houten
VM
,
Kummer
JA
,
Snel
MH
,
Snijders
PJ
, et al
Persistence of genetically altered fields in head and neck cancer patients: biological and clinical implications
.
Clin Cancer Res
2001
;
7
:
1523
32
.
8.
van Houten
VM
,
Tabor
MP
,
van den Brekel
MW
,
Kummer
JA
,
Denkers
F
,
Dijkstra
J
, et al
Mutated P53 as a molecular marker for the diagnosis of head and neck cancer
.
J Pathol
2002
;
198
:
476
86
.
9.
Tabor
MP
,
Brakenhoff
RH
,
Ruijter-Schippers
HJ
,
Van Der Wal
JE
,
Snow
GB
,
Leemans
CR
, et al
Multiple head and neck tumors frequently originate from a single preneoplastic lesion
.
Am J Pathol
2002
;
161
:
1051
60
.
10.
Braakhuis
BJ
,
Tabor
MP
,
Kummer
JA
,
Leemans
CR
,
Brakenhoff
RH
.
A Genetic explanation of slaughter's concept of field cancerization: evidence and clinical implications
.
Cancer Res
2003
;
63
:
1727
30
.
11.
Oliveira
MVMd
,
Fraga
,
Pereira
CS
,
Barros
LO
,
Oliveira
ES
,
Guimarães
ALS
, et al
Field cancerization in head and neck squamous cell carcinoma: immunohistochemical expression of P53 and Ki67 proteins: Clinicopathological Study
.
Rev Clín Pesq Odontol(Impr)
2010
;
6
:
17
27
.
12.
Poh
CF
,
Zhang
L
,
Anderson
DW
,
Durham
JS
,
Williams
PM
,
Priddy
RW
, et al
Fluorescence visualization detection of field alterations in tumor margins of oral cancer patients
.
Clin Cancer Res
2006
;
12
:
6716
22
.
13.
Heaphy
CM
,
Griffith
JK
,
Bisoffi
M
.
Mammary field cancerization: molecular evidence and clinical importance
.
Breast Cancer Res Treat
2009
;
118
:
229
39
.
14.
Nonn
L
,
Ananthanarayanan
V
,
Gann
PH
.
Evidence for field cancerization of the prostate
.
Prostate
2009
;
69
:
1470
9
.
15.
Graham
TA
,
McDonald
SA
,
Wright
NA
.
Field cancerization in the GI tract
.
Future Oncol
2011
;
7
:
981
93
.
16.
Squier
CA
,
Kremer
MJ
.
Biology of oral mucosa and esophagus
.
J Natl Cancer Inst Monogr
2001
:
7
15
.
17.
Clayton
E
,
Doupe
DP
,
Klein
AM
,
Winton
DJ
,
Simons
BD
,
Jones
PH
.
A single type of progenitor cell maintains normal epidermis
.
Nature
2007
;
446
:
185
9
.
18.
Doupe
DP
,
Alcolea
MP
,
Roshan
A
,
Zhang
G
,
Klein
AM
,
Simons
BD
, et al
A single progenitor population switches behavior to maintain and repair esophageal epithelium
.
Science
2012
;
337
:
1091
3
.
19.
Martincorena
I
,
Roshan
A
,
Gerstung
M
,
Ellis
P
,
Van Loo
P
,
McLaren
S
, et al
Tumor evolution. High burden and pervasive positive selection of somatic mutations in normal human skin
.
Science
2015
;
348
:
880
6
.
20.
Sprouffske
K
,
Pepper
JW
,
Maley
CC
.
Accurate reconstruction of the temporal order of mutations in neoplastic progression
.
Cancer Prev Res
2011
;
4
:
1135
44
.
21.
Takeda
T
,
Sugihara
K
,
Hirayama
Y
,
Hirano
M
,
Tanuma
JI
,
Semba
I
.
Immunohistological evaluation of Ki-67, P63, Ck19 and P53 expression in oral epithelial dysplasias
.
J Oral Pathol Med
2006
;
35
:
369
75
.
22.
Kushner
J
,
G
,
Jordan
RC
.
Patterns of P53 and Ki-67 protein expression in epithelial dysplasia from the floor of the mouth
.
J Pathol
1997
;
183
:
418
23
.
23.
Gale
N
,
Zidar
N
,
Poljak
M
,
Cardesa
A
.
Current views and perspectives on classification of squamous intraepithelial lesions of the head and neck
.
2014
;
8
:
16
23
.
24.
Reed
AL
,
Califano
J
,
Cairns
P
,
Westra
WH
,
Jones
RM
,
Koch
W
, et al
High frequency of P16 (Cdkn2/Mts-1/Ink4a) inactivation in head and neck squamous cell carcinoma
.
Cancer Res
1996
;
56
:
3630
3
.
25.
Smeets
SJ
,
Braakhuis
BJ
,
Abbas
S
,
Snijders
PJ
,
Ylstra
B
,
van de Wiel
MA
, et al
Genome-wide DNA copy number alterations in head and neck squamous cell carcinomas with or without oncogene-expressing human papillomavirus
.
Oncogene
2006
;
25
:
2558
64
.
26.
Smeets
SJ
,
van der Plas
M
,
Schaaij-Visser
TB
,
van Veen
EA
,
van Meerloo
J
,
Braakhuis
BJ
, et al
Immortalization of oral keratinocytes by functional inactivation of the P53 and pRb pathways
.
Int J Cancer
2011
;
128
:
1596
605
.
27.
Durrett
R
,
Moseley
S
.
Spatial Moran models I. Stochastic tunneling in the neutral case
.
Ann Appl Probab
2015
;
25
:
104
15
.
28.
Durrett
R
,
Foo
J
,
Leder
K
.
Spatial Moran models, II: cancer initiation in spatially structured tissue
.
J Math Biol
2016
;
72
:
1369
400
.
29.
Foo
J
,
Leder
K
,
Ryser
MD
.
Multifocality and recurrence risk: a quantitative model of field cancerization
.
J Theor Biol
2014
;
355
:
170
84
.
30.
Bramson
M
,
Griffeath
D
.
On the Williams-Bjerknes Tumor-Growth Model .2
.
Math Proc Cambridge
1980
;
88
:
339
57
.
31.
Elrefaey
S
,
Massaro
MA
,
Chiocca
S
,
Chiesa
F
,
Ansarin
M
.
HPV in oropharyngeal cancer: the basics to know in clinical practice
.
Acta Otorhinolaryngol Ital
2014
;
34
:
299
309
.
32.
Attner
P
,
Du
J
,
Nasman
A
,
Hammarstedt
L
,
Ramqvist
T
,
Lindholm
J
, et al
The role of human papillomavirus in the increased incidence of base of tongue cancer
.
Int J Cancer
2010
;
126
:
2879
84
.
33.
Boyle
P
,
Macfarlane
GJ
,
Blot
WJ
,
Chiesa
F
,
Lefebvre
JL
,
Azul
AM
, et al
European school of oncology advisory report to the European Commission for the Europe against Cancer Programme: Oral Carcinogenesis in Europe
.
Eur J Cancer B Oral Oncol
1995
;
31B
:
75
85
.
34.
Hashibe
M
,
Brennan
P
,
Chuang
SC
,
Boccia
S
,
Castellsague
X
,
Chen
C
, et al
Interaction between tobacco and alcohol use and the risk of head and neck cancer: pooled analysis in the International Head and Neck Cancer Epidemiology Consortium
.
Cancer Epidemiol Biomarkers Prev
2009
;
18
:
541
50
.
35.
Agaku
IT
,
King
BA
,
Dube
SR
,
Control CfD, Prevention
.
Current cigarette smoking among adults—United States, 2005–2012
.
MMWR Morb Mortal Wkly Rep
2014
;
63
:
29
34
.
36.
Freedman
KS
,
Nelson
NM
,
Feldman
LL
.
Smoking Initiation among young adults in the United States and Canada, 1998–2010: a systematic review
.
Prev Chronic Dis
2012
;
9
:
E05
.
37.
Goel
RK
,
Nelson
MA
.
The effectiveness of anti-smoking legislation: a review
.
J Econ Surv
2006
;
20
:
325
55
.
38.
Weinstock
YE
,
Alava
I
III
,
Dierks
EJ
.
Pitfalls in determining head and neck surgical margins
.
Oral Maxillofac Surg Clin North Am
2014
;
26
:
151
62
.
39.
Varvares
MA
,
Poti
S
,
Kenyon
B
,
Christopher
K
,
Walker
RJ
.
Surgical margins and primary site resection in achieving local control in oral cancer resections
.
Laryngoscope
2015
;
125
:
2298
307
.
40.
Luebeck
EG
,
Moolgavkar
SH
.
Multistage carcinogenesis and the incidence of colorectal cancer
.
Proc Natl Acad Sci U S A
2002
;
99
:
15095
100
.
41.
Meza
R
,
Jeon
J
,
Moolgavkar
SH
,
Luebeck
EG
.
Age-specific incidence of cancer: phases, transitions, and biological implications
.
Proc Natl Acad Sci U S A
2008
;
105
:
16284
9
.
42.
Gallaher
J
,
Babu
A
,
Plevritis
S
,
Anderson
AR
.
Bridging population and tissue scale tumor dynamics: a new paradigm for understanding differences in tumor growth and metastatic disease
.
Cancer Res
2014
;
74
:
426
35
.
43.
Luebeck
EG
,
Curtius
K
,
Jeon
J
,
Hazelton
WD
.
Impact of tumor progression on cancer incidence curves
.
Cancer Res
2013
;
73
:
1086
96
.
44.
Holmberg
L
,
Garmo
H
,
Granstrand
B
,
Ringberg
A
,
Arnesson
LG
,
Sandelin
K
, et al
Absolute risk reductions for local recurrence after postoperative radiotherapy after sector resection for ductal carcinoma in situ of the breast
.
J Clin Oncol
2008
;
26
:
1247
52
.
45.
Kong
I
,
Narod
S
,
Taylor
C
,
Paszat
L
,
R
,
Nofech-Moses
S
, et al
Age at diagnosis predicts local recurrence in women treated with breast-conserving surgery and postoperative radiation therapy for ductal carcinoma in situ: a population-based outcomes analysis
.
Curr Oncol
2013
;
21
:
96
104
.
46.
James
BYM
,
Gross
CP
,
Wilson
LD
,
Smith
BD
.
Nci seer public-use data: applications and limitations in oncology research
.
Oncology
2009
;
23
:
288
.
47.
Brouwer
AF
,
Eisenberg
MC
,
Meza
R
.
Age effects and temporal trends in HPV-related and HPV-unrelated oral cancer in the United States: a multistage carcinogenesis modeling analysis
.
PLoS One
2016
;
11
:
e0151098
.
48.
Marron
M
,
Boffetta
P
,
Zhang
ZF
,
Zaridze
D
,
Wunsch-Filho
V
,
Winn
DM
, et al
Cessation of alcohol drinking, tobacco smoking and the reversal of head and neck cancer risk
.
Int J Epidemiol
2010
;
39
:
182
96
.
49.
Armitage
P
,
Doll
R
.
The age distribution of cancer and a multi-stage theory of carcinogenesis
.
Br J Cancer
1954
;
8
:
1
12
.
50.
Moolgavkar
SH
,
Dewanji
A
,
Venzon
DJ
.
A stochastic two-stage model for cancer risk assessment. I. The hazard function and the probability of tumor
.
Risk Anal
1988
;
8
:
383
92
.
51.
Dewanji
A
,
Venzon
DJ
,
Moolgavkar
SH
.
A stochastic two-stage model for cancer risk assessment. II. The number and size of premalignant clones
.
Risk Anal
1989
;
9
:
179
87
.
52.
Luebeck
EG
,
Moolgavkar
SH
.
Stochastic analysis of intermediate lesions in carcinogenesis experiments
.
Risk Anal
1991
;
11
:
149
57
.
53.
Curtius
K
,
Hazelton
WD
,
Jeon
J
,
Luebeck
EG
.
A multiscale model evaluates screening for neoplasia in Barrett's Esophagus
.
PLoS Comput Biol
2015
;
11
:
e1004272
.
54.
Beerenwinkel
N
,
Antal
T
,
Dingli
D
,
Traulsen
A
,
Kinzler
KW
,
Velculescu
VE
, et al
Genetic progression and the waiting time to cancer
.
PLoS Comput Biol
2007
;
3
:
e225
.
55.
Foo
J
,
Leder
K
,
Michor
F
.
Stochastic dynamics of cancer initiation
.
Phys Biol
2011
;
8
:
015002
.