Sustained Coevolution in a Stochastic Model of Cancer–Immune Interaction Free

T-cell immunotherapy has revolutionized modern cancer therapy, delivering durable remission outcomes with higher rates of overall survival for many cancer subtypes (1, 2). The tumor–immune interaction is quite complex, determined by the recruitment of a variety of immune system effectors, and varies as a function of disease duration, progression, and subtype (3, 4). A majority of observed and therapeutic antitumor immunity is accomplished by the cytotoxic CD8⁺ T-cell repertoire, wherein immune cells recognize tumor-associated antigens (TAA) present on the surface of cancer cells (5). Adequate recognition leads to elimination of the target cells and expansion of the recognizing T-cell clone. Consistent with theoretical expectations (6, 7), effective immune repertoire detection of TAAs has been empirically validated during natural cancer progression as well as in the postimmunotherapeutic setting (8, 9).

One critical phenomenon that determines ultimate disease outcome is immunosurveillance, which relates to the degree to which tumor progression (typically in the absence of immunotherapy) is controlled by cancer–immune coevolution. To date, a majority of our understanding of this complex interaction has been experimentally driven (10–13) with immunoediting and TAA negative selection driving early-stage disease (14, 15). Recognition of TAAs relies on immune system discrimination between self and non-self epitopes, which is largely achieved through host-directed tolerance mechanisms. First, thymic negative selection functions by deleting immature T cells, which recognize common self-epitopes, thereby imparting central tolerance to the T-cell repertoire (16). This process is incomplete, however, and necessitates subsequent peripheral tolerance mechanisms (17, 18). Peripheral tolerance, in contrast, relies on T-cell inhibition at a systems-level and is achieved whenever a threat is perceived to be small based on a lower detection limit (19).

The “growth-threshold conjecture” is a key defining feature of systems-level immune detection wherein the total rate of change, not level, of a growing threat is the limiting determinant of immune system recognition (20, 21). Consequently, pathogens with small individual growth rates may slowly expand to large population sizes before risking immune system detection, consistent with the observed indolent course of intracellular viral illnesses (22) and escape of slow-growth tumors (23, 24). The corresponding evolutionary trajectories that arise as a result of immune detection based on threat dynamics are quite diverse and have implications for disease prognosis (15, 25). A mathematical framework that formalizes the underlying probabilistic dynamics governing the tumor–immune interaction promises to shed light on this complex process and could ultimately be implemented to better guide patient-specific immunotherapy.

The idea of modeling tumor progression and, in particular, treatment evasion is not new. Well-established stochastic models have revolutionized our understanding of acquired drug resistance in the setting of targeted therapy, in part, by highlighting the importance of random mutant acquisition and clonal selection on the population level (26, 27). Prior studies have also analyzed tumor escape using ordinary differential equations to model the tumor–immune interaction (28–30). A more recent analysis investigated systems-level cancer-immune dynamics in a deterministic setting without taking into account any evasion (23). Because the immune repertoire and evading cancer population are both evolving systems capable of random and repeated recognition and evasion, stochastic analysis is a useful framework for conceptualizing tumor recognition by the adaptive immune system. We previously considered an extreme version of this problem, namely an adaptive immune analogue susceptible to permanent defeat by a single but durable stochastic evasion event (24). Here we consider a more realistic scenario and use it to interpret data regarding age incidence curves as well as tumor clonal evolution.

Specifically, we introduce an analytically tractable mathematical framework for modeling the interaction between an adaptive T-cell repertoire capable of repeated recognition and a cancer population that may repeatedly evolve mechanisms of immune evasion. We solve for the cancer escape and elimination probabilities analytically and provide a statistical framework to study the mean-variance profiles of observed cancer coevolutionary trajectories. We apply our model to early cancer evolutionary data and predict “common initiation and rare progression,” wherein a majority of cancer-initiating events are controlled by an intact immune system, of which, only a small subset randomly escapes and leads to observed disease.

Our basic tumor–immune model consists of clones of cancer cells with pure birth at rate r until they are detected by the immune system. Recognition may occur once the tumor has reached a minimal detection size m₀ and we work assuming either all clones are detected at size m₀ with a one-time detection probability q (referred to as deterministic detection) or an ongoing chance of detection at sizes above m₀ (adaptive detection). Once detection occurs, cells are killed at rate d > r, giving rise to a finite lifetime for that clone. There is also a probability per birth event μ << 1 that the clone will give rise to an evasive tumor subclone. We previously investigated a simpler model that assumed evasion events to be complete and durable against recognition by all current and future T cells (24). Effectively, evading clones were rare but catastrophic cells that were invisible to the immune system, modeling an extreme event such as MHC-I downregulation. In general, it is more reasonable to assume that most evasion events impart transient escape via cancer antigen downregulation, for example, from recognition by the current dominant effector T-cell clone, with future detection of less-frequent or de novo TAAs possible by other T-cell clones.

We study here the clonal evolution that results from repeated recognition and evasion events. Elimination of the parent clone follows from its recognition and eventual elimination by the immune system. However, the cancer population evolves if it successfully acquires an immune-evasive daughter clone prior to elimination. This process repeats for each subsequent population until either elimination or escape occurs, and we refer henceforth to the periodic growth followed by detection and elimination as a recognition cycle. To accomplish this, we assume that any evasive clonal population has a chance of being detected upon maturity to detection size. A typical trace of the population dynamics generated by our model is shown in Fig. 1A–C where clones are created by evasion and die by detection over repeated cycles. Note that when the tumor population reaches some higher threshold M the cancer becomes clinically observable and “escapes” immune surveillance; in this simulation this occurs at around t = 50. Here the color code reflects subclones created during a specific cycle of the process. In general, a given clone can give rise to multiple evasive lineages during its lifetime. For the deterministic detection scheme, it is straightforward to argue (24) that the number of such daughter clones may be approximated by a Poisson distribution with parameter

Thus, the total evader “intensity,” that is, the parameter governing the total number of clones produced at period n, is simply γ_n,j = jλ, where j is the number of evasive clones at period n. A generalization of this formula for the adaptive assumption is presented in the Supplementary Data. The model can therefore be framed as generalization of the Galton–Watson branching process, well-established in modeling biological systems (31–33), by the addition of possible cancer escape. The process thus eventually culminates in either escape to size M or complete eradication of the tumor. We shall refer to the expected per-capita progeny as the branching parameter. In the deterministic case, this value is simply λ from Eq. A.

Instead of tracking the population size at a given time t along with the distribution of arriving clones (Supplementary Fig. S1), we focus instead on the number of evading clones, Z_n, per recognition cycle, n, starting from one initial clone (Fig. 2A). The clearance probability, or likelihood that a clone is recognized at cycle n, is denoted by q_n. This framework can in general handle any reasonable assumptions on the dependence of detection on cycle number. For example, it may be the case that the immune system becomes less capable over time due to exacerbation by an increasing number of unique clones. One could also generalize the above treatment to allow for the recognition rate to depend explicitly on the number of clones present at cycle n or even on the individual evasive clone as a function of its appearance order. These assumptions lead to an inhomogeneous Markov process relating the number of clones at cycle n to the number of clones at cycle n + 1 (Fig. 2B and C). For simplicity, we will mostly assume that q_n is a constant, which will be denoted as q_c, and in this case the Markov chain is homogeneous in cycle number. The probability of generating k clones in the next cycle given the existence of j clones in the current cycle is

with

The immune system must recognize every subclone at each period or else the threat escapes. Thus, for j clones at period n, the escape probability is

Finally, the probability of tumor elimination is just the probability that no evasive clones are generated and no escape occurs, namely

All simulations, data analysis, and data visualization were performed using MATLAB version R2017b.

Derivations of the escape and elimination probabilities, together with the branching parameter, were solved using standard probabilistic analysis.

All simulations of the evolutionary process were performed using a modified Gillespie algorithm. Analyses involving escape were simulated for all surviving populations until the emergence of a dominant clone that exceeded the upper clinical escape threshold.

Age-specific cancer incidence data for breast, bladder, kidney, pancreas, prostate, head and neck, ovary, lung, melanoma, acute myeloid leukemia, and chronic lymphocytic leukemia were obtained from Cancer Research UK. Age incidence for triple-negative breast cancer was obtained by applying Web Plot Digitizer (34) to Fig. 1 of ref. 35 and taking a weighted average based on the empirical frequency of reported groups (Asian, n = 368; Black, n = 509; Hispanic, n = 920; White, n = 2521). Least-squares linear regression was performed on all available incidence data between ages 0 and 40 and in a similar manner between the age at which incidence is first measured and 40. The corresponding least-squares slope parameter is plotted for each cancer subtype against the mutation rate.

Estimates of the per-cell mutation rates were obtained by using median rates from the available estimates from ref. 36 for acute myeloid leukemia (mutation rate μ = 3.33 × 10⁻⁷; n = 134), bladder (μ = 5.17 × 10⁻⁶; n = 35), breast cancer (μ = 9.33 × 10⁻⁷; n = 121), chronic lymphocytic leukemia (μ = 7.67 × 10⁻⁷; n = 91), head and neck cancer (μ = 3.17 × 10⁻⁶; n = 181), lung cancer (μ = 7.20 × 10⁻⁶; n = 514), melanoma (μ = 1.32 × 10⁻⁵; n = 121), ovarian cancer (μ = 1.65 × 10⁻⁶; n = 394), pancreatic cancer (μ = 1.07 × 10⁻⁶; n = 13), prostate cancer (μ = 7.33 × 10⁻⁷; n = 227), and renal cancer (μ = 1.53 × 10⁻⁶; n = 225).

Correlations between the linear slope parameter of early cancer age incidence and per-cell evasion rates were assessed using Pearson correlation coefficient, and reported P values were obtained. For this, the null hypothesis was that there existed no correlation between the slope of early cancer age incidence and the per-cell evasion rate. Rejection of the null hypothesis and the presence of significant (nonzero) correlation was assessed at a significance value of α = 0.05.

The time to most recent common ancestor relative to cancer initiating time for renal clear cell carcinoma was estimated from the TRACERx renal consortium using previously published data (37). The data were separated into two categories: evasion dominated by branching (samples K59, K167, K146, K065, K003, K150, K104, K163, K153, K126, K108, K135, K169, K139, K145, K156, K162, K180) and primarily linear evolution (samples K113, K165, K130, K023, K176, K037, K143, K124, K136, K096, K027, K158, K105) based on samples having their most recent common ancestor earlier or within 10 years of clinical detection, respectively. Mean and variance estimates of the mutational burden in patients with lung cancer were calculated from data reported in (38). Mean-variance frontiers were plotted assuming the quadratic relationship of Eq. J parameterized by some unknown fraction of relevant evasion events α, where we considered α up to 10%. Immune-enhanced, immune-compromised, and immune-neutral regions were defined by calculating all allowable mean-variance frontiers based on trajectories generated by the domain of interest.

The following sections present the main findings of our analysis. Full mathematical details are provided in the Supplementary Data.

Let E_n (resp. F_n) be the event that the threat escapes (resp. is eliminated) at period n. One can directly derive a backwards discrete-time master equation that determines the long-time limit of these splitting probabilities. Specifically focusing on the elimination and taking q = q_c to be constant, we have

where P_k (F_n) is the probability that a tumor trajectory having k clones at cycle n eventually leads to complete tumor elimination. Clearly, each clone is independent from its future trajectory and hence P_k (F_n) = P₁ (F_n)^k. Thus, Eq. F becomes a closed-form equation for P₁ (F_n)

Taking the steady-state limit of this equation yields an equation for the asymptotic elimination probability p*

One can similarly find the asymptotic value of the escape probability and verify that it is just 1- p*; in other words, the system always chooses one of the absorbing states in the long-time limit. A more complete treatment presented in the Supplementary Data allows this result to be extended to evaluate the full time-dependent value of P₁ (F_n) along with more general assumptions for the q_n.

The above outlines the model for the deterministic case since detection, should it occur, must do so at size m₀. Adaptive detection incorporates immune decline via an immunomodulatory parameter ν > 0 that results in larger average detection sizes and yields a similar elimination probability |${\tilde{p}^*}$| (See Supplementary Data for full details). This analytic framework above agrees well with results obtained from simulating the full coevolutionary process (Fig. 3A and B). “Escape by underwhelming” is an important experimentally observed feature wherein threats of intermediate growth rates are at a survival disadvantage relative to their faster and slower growing counterparts (39). Our model recapitulates this behavior for minimal detection thresholds m₀ limited by the total cancer growth rate (Fig. 3C), in agreement with the growth-threshold conjecture (22). This simple theoretical framework generates a significant degree of diversity in dynamical behavior, governed completely by the clearance probability and the branching parameter (Fig. 3D; Supplementary Figs. S2–S10).

Our foundational model may also be adapted to investigate a variety of assumptions on the nature of immune escape, including declines in immune clearance as a function of elapsed time and total surviving clone numbers (Sec. S4.4, Intertemporal immune decline; Sec. S4.5, Clone frequency-dependent recognition) along with their effects on the likelihood of cancer elimination and escape (Supplementary Figs. S2 and S3). Prior empirical studies have proposed the existence of an equilibrium state, wherein a cancer population is neither fully eliminated nor grows to overwhelm the host (40, 41). Surprisingly, our framework would suggest that ultimate coexistence is not reasonably predicted, occurring only when unreasonable assumptions are imposed on the clearance probabilities (Sec. S4.3).

As proof-of-principle that local physiologic changes in immune function may affect observed cancer frequency, even among healthy individuals, we study the effects of physiologic declines in immune turnover on increases in early age incidence for various cancer types using publicly available large population datasets (Supplementary Fig. S6; ref. 42). Previous studies have shown that age-incidence data can be fit to simple models assuming a declining immune system with age (6, 43). Our dynamical model predicts, all else being equal, that tumor escape likelihoods vary directly as a function of evasion rate μ implicitly through the branching parameter. Although evasion need not be limited exclusively to genomic mechanisms, we use experimentally derived per-cell mutation rates obtained from a large pan-cancer analysis as an indirect measure of evasion (36).

We calculate the slope parameter of linear regression for early age incidence (Fig. 4A; Supplementary Fig. S6) and find a strong correlation between these values and the evasion rate. Aggregate breast cancer incidence appears to be an exception, perhaps a result of early screening detection and the incidence following puberty of hormone-sensitive disease. Restricting the comparison to triple negative breast cancer (35) reveals statistically significant correlation between incidence increases and evasion rate (Fig. 4B; Supplementary Fig. S7). Decreases in immune function modeled via an increased adaptive parameter ν (see Supplementary Data), assumed directly proportional to early age increases, result in increased ultimate escape probability, |$1 - {\tilde{p}^*}$|⁠. We find excellent agreement in the predicted incidence slope-μ profile for the best-fit scaled range of immune function parameter ν (Fig. 4C). Moreover, using this normalized range for ν in the model predicts a linear increase in incidence for all observed mutation rates (Fig. 4D). Together, these findings demonstrate that cancer progression and escape are well-explained by variability in immune function across nearly all cancer types.

Given the ability of our mathematical model to link observed cancer evasion and recognition cycle timing to underlying tumor–immune coevolutionary dynamics, we next assessed the extent of immune evasion and clearance observed in cancer evolutionary data. The TRACERx renal clear cell carcinoma (RCCC) dataset (37) was particularly useful, providing estimates on the timing of landmark evolutionary events. Surprisingly, and consistent with model behavior under sustained control (Fig. 3D, q_c = 0.9, λ = 0.96), the authors found extended periods of clonal evolution sustained by small population sizes (estimated to be of order 10² cells). Using the arrival time of the most recent common ancestor (MRCA) as a representative of escape, we estimated the expected number of recognition cycles given the observed intervening time between MRCA detection and disease initiation. If t_m is the time it takes a clone to grow to minimal detection size X(t_m) ∼ 10² and T_M is the times it takes to grow to ultimate detection size X(t_m) ∼ 10⁹, then the intervening time, assuming n recognition cycles occurs, is given by Δt = t_M − nt_m. From this, we can estimate the number of recognition cycles via

where nt_m and t_M are the observed times (see Supplementary Data for details). Our calculations predict that, on average, RCCC undergoes approximately 27 distinct, immunologically relevant recognition cycles prior to escape (Fig. 5A and B; Supplementary Fig. S11A).

The parameter range that generates this behavior is near criticality |$( {\lambda \undertilde{ \lt } 1} )$| with high clearance probabilities (q_c approaching 1), suggesting effective immune handling of renal cell threats (Fig. 5C; Supplementary Fig. S11B). To further investigate the extent of immune protection, we used parameter estimates consistent with these findings to assess the frequency of observed cancer incidence relative to the total number of (unobserved) cancer-initiating events, and predict that initiation is 18 times as frequent as measured incidence. These results suggest that the adaptive immune system actively filters many potential threats manifesting as de novo initiating cancer clones, but may occasionally miss due to statistical chance, in support of the “bad-luck hypothesis” (44). The resulting population dynamics closely resemble those observed in the cancer incidence data, with a large disease period consisting of clonal evolution with low tumor burden, followed by either elimination or escape with extensive subclonal evolution (Fig. 5D). When separating samples into categories based on extensive linear or branched evolution and repeating this analysis, we find that patients with disease dominated by linear evolution are predicted to have significantly larger cycle numbers (57 ± 35) with a concomitant increase in the number of predicted initiating events relative to observed incidence. This behavior contrasts with trajectories dominated by branched evolution, which have a reduced cycle number (5 ± 2) with increased escape probability (Supplementary Fig. S12).

Given the predicted implications of immune performance on eliminating early and indolent disease, we wanted to investigate the relationship between disease subtype and predicted immune status. Available data on large cohort patients with lung cancer (38) provided enough samples for each studied subtype to estimate mean and variance profiles for the distribution of clonal mutational events. This analysis was repeated for each lung cancer subtype and then compared with all predicted values in the model. Because only a small (unknown) fraction α << 1 of total mutations are expected to impart immune evasion, we consider the resulting allowable mean-variance frontier assuming a range of small α. Very briefly, for X being the total number of clonal nonsynonymous mutational events with respective mean and variance given by μ_X and |$\sigma _X^2$|⁠, the mean and variance for the fraction αX that are relevant to immune evasion relate to the observed parameters μ_X and |$\sigma _X^2$| via

Equation J defines a quadratic mean–variance curve parameterized by α that can be applied to compare the observed curves against all allowable simulated parameters in the domain considered (Fig. 5A and C) to determine model conditions that are consistent with empirical data. Our above model defines a baseline relationship that well characterizes the reported mean-variance profiles across all lung cancer evolutionary data subtypes (Fig. 6; Supplementary Table S1). However, several variations in the mean-variance frontier are only achievable in our generalized model assuming either immune suppression or heightened immune surveillance (see Supplementary Data for full details). Interestingly, squamous cell carcinoma, and to a lesser extent positive smoking status, are more consistent with a process that is significantly immunosuppressed, suggesting that the observed pattern of mutational variability emerges in the setting of reduced immune targeting efficiency. Intriguingly, the predicted mean-variance frontier requires enhancements to immune clearance rates over successive periods to explain the nonsmoker signature (Eq. S123).

The underlying early clonal architecture of an evolving cancer population contains a significant amount of information related to cancer progression and the extent of immunosurveillance. Understanding this complex behavior is essential for better predicting the extent of tumor immunoediting, measuring the frequency of disease progression relative to incidence, estimating immune targeting efficacy, and proposing optimized treatment strategies conditioned on post-escape clonal distributions. Here, we developed a stochastic population dynamical model of the battle between an evolving cancer population and the CD8⁺ T-cell immune repertoire. Prior models of the cancer–immune interaction have explicitly represented and tracked immune effectors as a distinct population (22, 29). While useful for describing population control in a variety of contexts, including nonimmunogenetic tumors where immune evasion is a nondominant contributor to cancer progression, this approach falls short in characterizing the phylogenetic architecture and elimination probabilities of evasive cancer subclones following recognition via the adaptive immune system. Our alternative formulation explicitly models recognition on the population level assuming that the cancer compartment may randomly acquire evasive clones, each of which requires distinct T-cell clones in the immune repertoire for recognition and control. We made the implicit assumptions that cancer adaptivity is represented by an evasion rate μ, while immune system adaptivity depends on the repertoire diversity and turnover rates. The resulting tumor escape or elimination behavior is thus reduced to an effective branching parameter, λ. This framework is sufficiently general for studying the coevolution of tumor and immune system, in addition to formulating the key issue of immune escape and elimination with regards to the ability of the immune repertoire to keep pace with a heterogeneous and evading cancer population.

Previous empirical studies have proposed the existence of an equilibrium state wherein a growing cancer population can be maintained indefinitely (40, 41). In contrast, our model only predicts ultimate equilibrium for the unlikely case in which recognition rates of new clones increase dramatically as an increasing function of total clone number. More realistic assumptions on the decline in immune recognition, occurring either over time or when overwhelmed by a greater number of competing clones, prevents the cancer population from reaching ultimate equilibrium and implies that the immune recognition process ultimately results in either escape or elimination with certainty, given sufficient time.

The model's underlying dynamical behavior subsequently encodes information on both the likelihood and timing of ultimate cancer escape or elimination. Our framework highlights the importance of immunosurveillance in early disease progression as evidenced by differential increases in cancer incidence that scale with evasion rate across nearly all cancer types. We utilized per-cell mutation rates as an indirect measure of evasion capacity. Our findings were consistent despite changing the age range used as input for the linear regression. Indeed, many mutations are independent of immune evasion events and may increase the antigenic burden of the cancer population. This framework only distinguishes progeny based on an evasion event, which results in alteration or removal of an immune target recognized on the parent clone. Evaluating the competing effects of advantageous antigen down-regulation and undesirable tumor neoantigen generation on overall recognition by the immune system is an important consideration and warrants further investigation.

Using the excellent data available on the timing of early evolutionary events in RCCC, together with the observation that early clonal evolution is estimated to be sustained at small cell sizes, we show that the mean number of recognition cycles is quite large, suggestive of a drawn-out competition between immune and cancer compartments. Our results suggest that cancer progression from early initiation to appreciable disease only occurs a small fraction of the time, in support of our proposed “common initiation and rare progression” postulate, wherein most cancer initiating events are effectively controlled by host recognition mechanisms. In the minority of cases where the cancer population escapes, our analysis predicts this to be due to random chance, consistent with the “bad-luck” hypothesis (44). Moreover, partitioning of samples into groups based on trajectories dominated by linear versus branched evolution revealed differences in mean recognition cycle and relative escape likelihood, the latter group characterized by punctuated evolution and early disease. We remark that these results are most sensitive to the assumed minimal population size (here taken to be m₀ = 10² as estimated empirically); the resulting mean cycle number estimate scales as a decreasing function of m₀ as 2/log₁₀ m₀. Despite this, significant and prolonged coevolution is still predicted (15 recognition cycles in renal cancer patients) assuming an order of magnitude increase in the minimal detection size. Although empirically challenging, further investigation into the early tumor–immune interaction following initiation could refine this estimate and shed light on the underlying evolution of early disease for a variety of cancers. When comparing the mean-variance frontiers in the clonal distributions of TRACERx LCA subtypes, we found that observed evolutionary trajectories for squamous cell carcinoma and positive smoking status were consistent with simulated profiles assuming immune impairment, while nonsmokers and adenocarcinoma followed an opposite trend. Our results demonstrate the utility of further experimental investigation into the benefit of distinct treatment strategies based on cancer subtype and predicted immune function of each patient.

Cancer evolution is complex and highly variable owing to the dynamic interaction that occurs between tumor cells and the adaptive immune system. Improvement in our understanding of cancer progression therefore requires a detailed description of individual disease evolutionary trajectories. Toward this end, we have developed a mathematical framework that models a heterogeneous population of cancer cells subjected to adaptive immunosurveillance. When applied to empirical data, our model is a useful tool for studying the tumor–immune interaction as it establishes a fundamental relationship between evasion and the observed branching structure that emerges during disease progression. Although we have focused our analysis on describing the control and progression of cancer, our model is broadly applicable in understanding similar phenomena, such as intracellular infection by pathogens that evolve mechanisms of immune evasion. In HIV, for example, the delayed minimal disease burden produced by repeated recognition cycles in our model shares some general features of the latent period (45). Future efforts to apply this model in a broader context may provide insights regarding the observed dynamics of slow-progression intracellular threats.

No potential conflicts of interest were disclosed.

Conception and design: J.T. George, H. Levine

Development of methodology: J.T. George, H. Levine

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

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J.T. George, H. Levine

Writing, review, and/or revision of the manuscript: J.T. George, H. Levine

Study supervision: J.T. George

We would like to thank Jeffrey Molldrem for fruitful discussions on immunotherapy and cancer-immune coevolution.

J.T. George is supported by the NCI of the NIH (F30CA213878). H. Levine is supported by Physics Frontiers Center NSF Grant PHY-1427654 and by NSF PHY-1605817.

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.