Maria Bilous

Cindy Serdjebi

Arnaud Boyer

Claudia Pouypoudat

Pascale Tomasini

Dominique Barbolosi

Fabrice Barlesi

François Chomy

Sébastien Benzekry
email: sebastien.benzekry@inria.fr

Quantitative mathematical modeling of clinical brain metastasis dynamics in non-small cell lung cancer

described [START_REF] Gui | A prospective evaluation of whole brain volume loss and neurocognitive decline following hippocampal-sparing prophylactic cranial irradiation for limited-stage small-cell lung cancer[END_REF] . This highlights the need for a better prediction of patients with highest risk of BM development or relapse, in order to select patients with the higher benefit/risk ratio for WBRT. Several phase III trials were conducted but no firm conclusion applicable for the entire patient population could be drawn [START_REF] Mulvenna | Dexamethasone and supportive care with or without whole brain radiotherapy in treating patients with nonsmall cell lung cancer with brain metastases unsuitable for resection or stereotactic radiotherapy (QUARTZ): results from a phase 3, non-inferiority, randomised trial[END_REF][START_REF] Pechoux | Whole brain radiotherapy in patients with NSCLC and brain metastases[END_REF] . The prediction of the risk of BM development and indication for WBRT to prevent brain dissemination is particularly interesting for specific populations such as EGFR or ALK positive NSCLC patients, for whose incidence of BMs is high [START_REF] Soffietti | Diagnosis and treatment of brain metastases from solid tumors: guidelines from the European Association of Neuro-Oncology (EANO)[END_REF][START_REF] Doherty | Treatment options for patients with brain metastases from EGFR/ALK -driven lung cancer[END_REF][START_REF] Jiang | EGFR TKIs plus WBRT Demonstrated No Survival Benefit Other Than That of TKIs Alone in Patients with NSCLC and EGFR Mutation and Brain Metastases[END_REF] . This points to the need of rational tools to decide therapeutic action in a patient-specific way.

For such aims, quantitative mathematical modeling may be of considerable help, by inferring underlying processes from observed data or providing useful numerical tools in the era of personalized medicine [START_REF] Altrock | The mathematics of cancer: integrating quantitative models[END_REF][START_REF] Barbolosi | Computational oncology -mathematical modelling of drug regimens for precision medicine[END_REF] . However, despite numerous studies in mathematical and computational oncology, the majority of the efforts have to date focused on mathematical models at the scale of a single tumor [START_REF] Byrne | Dissecting cancer through mathematics: from the cell to the animal model[END_REF] . Moreover, the comparison of the models to empirical data remains infrequent. Historically, modeling efforts in the field of metastasis were first initiated by statistical models phenomenologically describing relapse hazards [START_REF] Bross | Statistical testing of a deep mathematical model for human breast cancer[END_REF][START_REF] Guiguet | Size distribution of metastases during detection and adjuvant treatment: biomathematical approach[END_REF][START_REF] Koscielny | A simulation model of the natural history of human breast cancer[END_REF][START_REF] Retsky | Computer simulation of a breast cancer metastasis model[END_REF] . In the 1970's, Liotta et al. pioneered the development of biologically-based and experimentally-validated models for all the main steps of the metastatic process [START_REF] Liotta | Stochastic model of metastases formation[END_REF] . Since then, only a small number of studies addressed this topic [START_REF] Iwata | A Dynamical Model for the Growth and Size Distribution of Multiple Metastatic Tumors[END_REF][START_REF] Hanin | A stochastic model for the sizes of detectable metastases[END_REF][START_REF] Michor | Stochastic dynamics of metastasis formation[END_REF][START_REF] Haeno | Computational Modeling of Pancreatic Cancer Reveals Kinetics of Metastasis Suggesting Optimum Treatment Strategies[END_REF][START_REF] Newton | A Stochastic Markov Chain Model to Describe Lung Cancer Growth and Metastasis[END_REF][START_REF] Scott | Mathematical modeling of the metastatic process[END_REF][START_REF] Hartung | Mathematical modeling of tumor growth and metastatic spreading: validation in tumor-bearing mice[END_REF][START_REF] Baratchart | Computational Modelling of Metastasis Development in Renal Cell Carcinoma[END_REF][START_REF] Benzekry | Modeling Spontaneous Metastasis following Surgery: An In Vivo-In Silico Approach[END_REF][START_REF] Poleszczuk | Abscopal Benefits of Localized Radiotherapy Depend on Activated T-cell Trafficking and Distribution between Metastatic Lesions[END_REF] . Of specific relevance to the current work, the Iwata model [START_REF] Iwata | A Dynamical Model for the Growth and Size Distribution of Multiple Metastatic Tumors[END_REF] introduced a size-structured population approach to capture the time development of a population of regional metastases from a hepatocellular carcinoma. It was further extended and studied from the mathematical and numerical viewpoints [START_REF] Barbolosi | Mathematical and numerical analysis for a model of growing metastatic tumors[END_REF][START_REF] Benzekry | Mathematical analysis of a two-dimensional population model of metastatic growth including angiogenesis[END_REF][START_REF] Benzekry | Mathematical and numerical analysis of a model for anti-angiogenic therapy in metastatic cancers[END_REF] , and was found able to reproduce experimental in vivo data [START_REF] Hartung | Mathematical modeling of tumor growth and metastatic spreading: validation in tumor-bearing mice[END_REF][START_REF] Baratchart | Computational Modelling of Metastasis Development in Renal Cell Carcinoma[END_REF][START_REF] Benzekry | Modeling Spontaneous Metastasis following Surgery: An In Vivo-In Silico Approach[END_REF] .

To understand biological dynamics of metastasis, animal studies allowing to capture the natural course of the disease [START_REF] Hartung | Mathematical modeling of tumor growth and metastatic spreading: validation in tumor-bearing mice[END_REF][START_REF] Baratchart | Computational Modelling of Metastasis Development in Renal Cell Carcinoma[END_REF] and the effect of therapeutic interventions in controlled environments [START_REF] Benzekry | Modeling Spontaneous Metastasis following Surgery: An In Vivo-In Silico Approach[END_REF] provide valuable data for quantitative analysis, although with limited clinical relevance. Experimental procedures are tedious and often only provide access to the total metastatic burden, neglecting its distribution into distinct metastatic lesions when no anatomical imaging modalities are employed [START_REF] Baratchart | Computational Modelling of Metastasis Development in Renal Cell Carcinoma[END_REF][START_REF] Benzekry | Modeling Spontaneous Metastasis following Surgery: An In Vivo-In Silico Approach[END_REF] .

Individual-level clinical data of metastatic development are also challenging to obtain because precise number and size of existing lesions are not routinely collected. Apart from the landmark work of Iwata et al. [START_REF] Iwata | A Dynamical Model for the Growth and Size Distribution of Multiple Metastatic Tumors[END_REF] (which only included one untreated patient with hepatocellular carcinoma), no study has yet provided a quantitative modeling analysis of longitudinal data of individual-level number and size of metastases. The purpose of the current work is to develop this approach for BMs dynamics in NSCLC, in order to address the biological and clinical questions stated above.

Results

pre-diagnosis natural course of lung primary tumors. We investigated two possible growth models for the natural course of lung PTs: exponential and Gompertz. Exponential growth is the simplest model expressing uncontrolled proliferation and is often well adapted to describe tumor growth kinetics during limited observation periods [START_REF] Friberg | On the growth rates of human malignant tumors: implications for medical decision making[END_REF] . However, it has been shown that on longer timeframes (typically under volumes increases of 100 to 1000-fold), the specific growth rate of tumors decreases [START_REF] Laird | Dynamics of Tumour Growth: Comparison of Growth Rates and Extrapolation of Growth Curve To One Cell[END_REF] , which is well captured by the Gompertz model [START_REF] Laird | Dynamics of Tumour Growth: Comparison of Growth Rates and Extrapolation of Growth Curve To One Cell[END_REF][START_REF] Norton | A Gompertzian model of human breast cancer growth[END_REF][START_REF] Benzekry | Classical mathematical models for description and prediction of experimental tumor growth[END_REF] . For calibration of the parameters, we used the data of the PT size and mean doubling time at diagnosis (see Methods). The latter was assumed to be histology-dependent and retrieved from a meta-analysis comprising a total of 388 adenocarcinomas and 377 squamous cell carcinomas [START_REF] Friberg | On the growth rates of human malignant tumors: implications for medical decision making[END_REF][START_REF] Detterbeck | Turning gray: the natural history of lung cancer over time[END_REF] (see Supplementary Table S1). For the Gompertz model (two parameters), we additionally assumed a carrying capacity of 10 12 cells [START_REF] Klein | Parallel progression of primary tumours and metastases[END_REF] . Under such assumptions, we obtained predictions of pre-diagnosis phases of 19 years versus 5.4 years in the exponential and Gompertz cases, respectively, for an adenocarcinoma with diameter 35 mm (median value from ref. [START_REF] Mujoomdar | Clinical Predictors of Metastatic Disease to the Brain from Non-Small Cell Lung Carcinoma: Primary Tumor Size, Cell Type, and Lymph Node Metastases1[END_REF] ), see Fig. S1. The first prediction seems highly unrealistic compared with previous reports estimating the age of lung primary tumors to be 3-4 years old [START_REF] Patrone | How long have I had my cancer, doctor? Estimating tumor age via Collins' law[END_REF] , based on a different method using time to recurrence [START_REF] Collins | Observations on growth rates of human tumors[END_REF] . The Gompertz prediction on the other hand is rather consistent with the literature range. Moreover, the obtained value of the cellular proliferation rate was realistic, corresponding to a length of the cell cycle

    =      α log(2) p 0,
of 24.4 hours [START_REF] Norton | Holland-Frei Cancer Medicine[END_REF] .

We therefore concluded that the Gompertz model was best adapted to describe the pre-diagnosis natural history of the PT. probability of brain metastasis as a function of primary tumor size can be described by a mechanistic model of metastasis. PT size was reported by Mujoomdar et al. [START_REF] Mujoomdar | Clinical Predictors of Metastatic Disease to the Brain from Non-Small Cell Lung Carcinoma: Primary Tumor Size, Cell Type, and Lymph Node Metastases1[END_REF] as a major predictive factor of BM in NSCLC. We investigated whether our mechanistic model of BM was able to fit the data from these authors, at the population level (see Methods for the definition of the model). To do so, we assumed a statistical distribution associated with the inter-individual variability for the metastatic dissemination parameter μ (see Methods). The growth parameters were determined as above, i.e. from histology only. This approach yielded size-dependent population probabilities of BM occurrence, which could be compared to the data. The model fits were reasonably accurate, for both histological types (Fig. 1). Importantly, this demonstrates the descriptive power of our model, which is based on mechanistic simulations of the underlying biological process. In this regard, this is superior to a mere fit of a biologically agnostic statistical model such as a linear regression. Notably, this result was obtained with a minimal number of parameters to describe the inter-patient heterogeneity in metastatic potential. Namely, only two parameters defining the lognormal statistical distributions of μ in the populations. The coefficients of variation were very large (>6,000%). The median value μ pop gives a quantitative way to measure median BM aggressiveness. It was found two log 10 orders of magnitude larger for adenocarcinomas than for squamous cell carcinomas (μ = . × - 1 39 10 ).

Model-based comparison of biological scenarios with individual longitudinal data suggests dormancy.

Data. These consisted of 10 and 11 PT size measurements, and 47 and 16 BM size measurements in two patients, respectively (one with an adenocarcinoma, the other with a squamous cell carcinoma). We first analyzed the adenocarcinoma patient (Fig. 2) to develop the model and used the second patient to investigate the generalizability of our findings. In the adenocarcinoma patient, the PT first responded to systemic therapy before slowly regrowing (Fig. 2A). A first distant BM was detected 20 months after diagnosis, and kept growing uncontrolled afterwards (Fig. 2B). Other BMs appeared during follow-up (Fig. 2C), reaching a total of 20 BMs at 48 months, date of last examination (Fig. 2D).

Growth kinetics. To model the effect of systemic treatment on the PT, we found that a tumor growth inhibition model [START_REF] Claret | Model-based prediction of phase III overall survival in colorectal cancer on the basis of phase II tumor dynamics[END_REF] (equation ( 1)) was able to adequately fit the data (Fig. 2A). Interestingly, for the BM kinetics, we found that the untreated PT Gompertz model predicted surprisingly well the observed data (Figs 2B andS2). Specifically, when employing the Gompertz parameters used for the pre-diagnosis PT history, and using initial conditions on each BM, the future growth curves were accurately predicted. This indicated that: 1) for this patient, the BMs did not respond to systemic therapy, at least from the time of first observation, 2) it is reasonable to assume that all BMs grew at the same growth rate, 3) the BMs growth rate might be similar to the PT growth rate, at least during the clinically overt phase.

Quantitative assessment of five biologically-based models of metastatic dissemination and colonization. However, this mere description of the BMs growth is not satisfactory as a model of systemic disease, since the dissemination part is absent. In particular, no model is given for the initiation times of the BMs. To include the dissemination component of the metastatic process, we relied on a modeling framework first initiated by Iwata et al. [START_REF] Iwata | A Dynamical Model for the Growth and Size Distribution of Multiple Metastatic Tumors[END_REF]

(see Methods).

Five models, corresponding to five possible biological scenarios, were investigated (Fig. 3). These included: (A) a "basic" model with same growth model and parameters for the PT and the BMs and no secondary dissemination (i.e. BMs spread only by the PT); (B) a model allowing different growth parameters for the PT and the BMs; (C) a model with secondary dissemination, i.e. the ability of BMs to spread BMs themselves 31 ; (D) a model allowing a delay before initiation of metastatic ability of the PT (the so-called linear model in which dissemination occurs at a late stage, opposed to the parallel model where dissemination is an early event [START_REF] Klein | Parallel progression of primary tumours and metastases[END_REF] ) and (E) a model allowing BM dormancy phases, i.e. time periods with null growth rate [START_REF] Chambers | Dissemination and growth of cancer cells in metastatic sites[END_REF][START_REF] Nguyen | Metastasis: from dissemination to organ-specific colonization[END_REF][START_REF] Aguirre-Ghiso | How dormant cancer persists and reawakens[END_REF] . Each of these scenario was translated into corresponding mathematical equations (see Methods).

The models exhibited different descriptive power, as quantified by the best-fit value of the objective function (Table 1). The best description of the data was obtained by the dormancy model (Fig. 4). Importantly, this model achieved accurate goodness-of-fit for 47 BM size measurements over 6 time points (Fig. 4A,B) with only 3 parameters fitted from the BM sizes data (μ, τ and γ). Consequently, parameter identifiability was found excellent (relative standard errors <5%, Table 2). Apparition of the first BM above the visibility threshold was also accurately described (19 vs 20 months, model vs data, Fig. 4C). The model overestimated the number of small visible BMs at 40 months, possibly because a lot of these were present but too small to be visible at MRI. Indeed, a lot of BMs were present at 43 months -in accordance with the model (Fig. 4A) -which thus emphasizes its predictive potential.

The basic (A), secondary (C) and delay models (D) exhibited substantially worse goodness-of-fit (Table 1 and Figs S3-5), suggesting that their underlying hypotheses are insufficient to provide valid explanations of BM dynamics. In the delay model, the fit was improved compared to models (A) and (B) but the value of the delay t d was found unrealistically large (4.8 years after the first cancer cell, corresponding to 6 months before diagnosis).

Finally, model (D) yielded a significant improvement of the goodness-of-fit for both the dynamics of the number of visible BMs and cumulative size distributions, while not deteriorating the practical identifiability of the parameters (Fig. S6). Of note, the estimated values of the BM growth parameters were close to the PT ones. However, given the previous observation that during the clinical phase the BMs grew at a growth rate consistent with the one predicted for the PT, combined with the facts that the optimal objective value was achieved for the dormancy model and that it is more parsimonious (one parameter less), we selected the dormancy model for further analysis. Clinically relevant simulations of the disease course reveal times of BM initiation. From the quantitative calibration of the dormancy model to the data, several predictions of clinical interest can be made. The value of γ that generated the best-fit was 1, suggesting equal probability for all the cells within the PT to disseminate. In turn, this could indicate of a well-vascularized tumor, which might be prone to anti-angiogenic therapy. Interestingly, the value of μ inferred from this patient-level data was in the same range as the one inferred from the above population analysis (μ = 2 × 10 -12 cell -1 .day -1 versus μ pop = 1.39 × 10 -11 cell -1 .day -1 ), giving further support to our approach. Once calibrated from the data, our model allowed to simulate the predicted natural history of the disease. The Supplementary Movie S1 shows a simulation of the PT growth together with the apparition of the entire population of BMs (visible + invisible). Stars indicate tumors that are invisible (<5 mm), and the BMs size distribution time course is also simulated. BMs represented in gray were born before the diagnosis, while BMs in white are the ones born after.

Apart from the age of the PT, prediction of birth times and invisible BM burden at any time point can be performed. Interestingly, we found that the first BM -which clinically appeared 19 months after diagnosis -was already present 14 months before clinical detection of the PT (Fig. 5A). In fact, at the time of diagnosis, our model predicted the presence of 10 occult BMs, representing a total burden of 1,167 cells mostly distributed into the largest (first) BM (size 1,046 cells  0.126 mm), see Fig. 5B. Notably, at the the PT had reached the imaging detection limit, no metastasis had occurred yet. This suggests that if the disease had been detected earlier (e.g., by systematic screening), BM occurrence might have been prevented provided the tumor had been operable at this time. The amount of predicted BMs present at T 2 was 28 tumors, with an overall BM burden of 1.7 × 10 7 cells (Fig. 5C). Therefore, provided that neuro-cognitive risks would be acceptable, the model would recommend WBRT or systemic treatment rather than localized intervention only.

Together, these results demonstrate the potential clinical utility of the model for prediction of the invisible BM, in order to inform therapeutic decision.

Generalization to a second patient.

To test whether the dormancy model was generalizable, we used data from a second patient, which was not employed during the model development phase. Given the different histological type of the lung PT (squamous cell carcinoma), we adapted the doubling time accordingly (Table S1) and found a younger age of the PT of 2.1 years. Estimation of the log-kill parameter from the TGI PT model during therapy suffered from lack of identifiability, due to an estimated short duration of treatment effectiveness (see parameter κ in Table 2). Response of the PT to therapy was indeed characterized by a faster relapse (parameter α 1 ), compared with the first patient. The qualitative structure of the model confirmed its descriptive power by where se is the standard error (square root of the covariance diagonal entry, see Materials an methods) and est is the parameter estimate. *For value of the parameter γ, due to identifiability issues, only five possible values were considered during the fit procedure (see Methods).

being able to give an accurate description the BMs size dynamics (Figs S7 and8), while the "basic" and "secondary dissemination" models were not (data not shown). Several parameters appeared to be patient-dependent, such as the dormancy duration τ, estimated to 171 days versus 133 days for the first patient. Due to the lower number of data points available for this patient, parameters identifiability was found worse, but still acceptable (Table 2).

Resulting clinical predictions were distinct (Fig. S8), emphasizing the patient-specific nature of BM dynamics. The first BM was clinically detected 14 months after diagnosis, but was predicted to have been disseminated 45 months prior to diagnosis. While for both patients cerebral dissemination had already occurred at the time of diagnosis, its extent was different with a much larger mass (>1,000 cells) for the first patient than for the second (8 cells). This is due to the long period of dormancy for patient 2, resulting in all 8 predicted BMs being still dormant at the time of diagnosis.

These results demonstrate the flexibility of our model, which is able to describe and predict distinct situations in terms of repartition of the metastatic burden into individual lesions.

Discussion

Using both population-and individual-level data of BM development in NSCLC, we have developed a general method based on biologically grounded computational models that allowed: 1) to infer the disease age from data on PT size and histology only, 2) to test multiple scenarios of metastatic dissemination and colonization against macroscopic data available in the clinic and 3) to infer times of BM initiations and number and size of invisible lesions during the clinical course of the disease. Estimation of the duration of the pre-diagnosis PT growth phase has important implications, since BMs are more likely to have occurred for an old PT compared to a young one. Our results showed a significant difference whether considering exponential or Gompertz growth, which was consistent with previous findings [START_REF] Detterbeck | Turning gray: the natural history of lung cancer over time[END_REF] . Of note, our age estimates of 5.3 and 2.1 years-old are in relative agreement with the 3-4 years range found by others relying on a different method [START_REF] Patrone | How long have I had my cancer, doctor? Estimating tumor age via Collins' law[END_REF] .

Importantly, we were then able to include this description of PT past growth into a mechanistic model to describe the probability of BM occurrence. We found that the cell-scale dissemination parameter μ was significantly larger in adenocarcinomas than in squamous cell carcinomas, which offers a quantitative theory for the reported differences in BM occurrence between these two histological types [START_REF] Mujoomdar | Clinical Predictors of Metastatic Disease to the Brain from Non-Small Cell Lung Carcinoma: Primary Tumor Size, Cell Type, and Lymph Node Metastases1[END_REF] .

More generally, our computational platform provided a way to translate biological findings into clinically useful numerical tools. However, in order to provide robust inference, the complexity of the models had to remain commensurate with the available data. Thus, several higher order phenomena relevant to the metastatic process were ignored or aggregated into mesoscopic parameters. For instance, the metastatic potential μ is the product of several cell-scale probabilities relating to the multiple steps of the metastatic cascade [START_REF] Gupta | Cancer metastasis: building a framework[END_REF] . The median values inferred from population analyses based on probabilities of BM (μ = .

× - 1 39 10 pop 11 and . × - 1 76 10 13 ) matched relatively well the ones that were found from analysis of patient-specific data (μ = × - 2 10 12 and . × - 1 02 10 12 ), given the standard deviations. The larger value found in patient 1 could be due to histology or to the EGFR mutational status, known to impact on the BM aggressiveness [START_REF] Tomasini | EGFR and KRAS Mutations Predict the Incidence and Outcome of Brain Metastases in Non-Small Cell Lung Cancer[END_REF][START_REF] Jiang | EGFR TKIs plus WBRT Demonstrated No Survival Benefit Other Than That of TKIs Alone in Patients with NSCLC and EGFR Mutation and Brain Metastases[END_REF][START_REF] Wang | Impacts of EGFR mutation and EGFR-TKIs on incidence of brain metastases in advanced non-squamous NSCLC[END_REF][START_REF] Han | A retrospective analysis in patients with EGFR-mutant lung adenocarcinoma: is EGFR mutation associated with a higher incidence of brain metastasis?[END_REF] . Interestingly, while relying on distinct modeling techniques, the value of μ that we obtained was in the same range as obtained by others using stochastic evolutionary modeling [START_REF] Haeno | Computational Modeling of Pancreatic Cancer Reveals Kinetics of Metastasis Suggesting Optimum Treatment Strategies[END_REF][START_REF] Reiter | Minimal functional driver gene heterogeneity among untreated metastases[END_REF] .

When evaluating the descriptive power of multiple biologically-based models against longitudinal data of BM sizes, our results suggested possible phases of dormancy. This finding is supported by several preclinical and clinical reports [START_REF] Nguyen | Metastasis: from dissemination to organ-specific colonization[END_REF][START_REF] Kienast | Real-time imaging reveals the single steps of brain metastasis formation[END_REF] . For instance, Nguyen et al. mention possible dormancy phases for cerebral metastases from lung cancer [START_REF] Nguyen | Metastasis: from dissemination to organ-specific colonization[END_REF] . Durations of these latent periods are mentioned to be of the order of months, which is consistent with our quantitative results (4.4 and 5.7 months for patients 1 and 2, respectively). BM dormancy could possibly be explained by the fact that colonizing lung carcinoma cells face a foreign micro-environment when reaching the brain, which hampers their establishment and growth [START_REF] Winkler | The brain metastatic niche[END_REF] . This could cause single cells to remain quiescent or tumors to remain at a small, pre-angiogenic size [START_REF] Chambers | Dissemination and growth of cancer cells in metastatic sites[END_REF][START_REF] Holmgren | Dormancy of micrometastases: balanced proliferation and apoptosis in the presence of angiogenesis suppression[END_REF] . Indeed, vascular endothelial growth factor-A (VEGF-A) has been shown to mediate BM dormancy in a mouse model of lung cancer [START_REF] Kienast | Real-time imaging reveals the single steps of brain metastasis formation[END_REF] .

Several second-order phenomena were ignored here for the sake of identifiability, but could nevertheless impact BM dynamics and dormancy. These include tumor-tumor interactions, either through soluble circulating factors [START_REF] Dewys | A Quantitative Model for the Study of the Growth and Treatment of a Tumor and Its Metastases with Correlation between Proliferative State and Sensitivity to Cyclophosphamide[END_REF] distantly inhibiting growth and possibly maintaining dormancy [START_REF] Holmgren | Dormancy of micrometastases: balanced proliferation and apoptosis in the presence of angiogenesis suppression[END_REF] , or by exchange of tumor cells between established lesions [START_REF] Gundem | The evolutionary history of lethal metastatic prostate cancer[END_REF][START_REF] Norton | Is cancer a disease of self-seeding?[END_REF] . We have recently proposed a model for distant interactions that was validated in a two-tumors experimental system and could be incorporated into the current modeling platform [START_REF] Benzekry | Mathematical Modeling of Tumor-Tumor Distant Interactions Supports a Systemic Control of Tumor Growth[END_REF][START_REF] Benzekry | Global Dormancy of Metastases Due to Systemic Inhibition of Angiogenesis[END_REF] .

While not problematic for the two patients that were investigated here since BMs did not respond to the systemic treatment -possibly because of the blood-brain barrier hampering delivery of the anti-cancer agents -a major limitation of our model is that the effect of therapy on the BMs was neglected. We intend to address this in future work, in particular for optimizing and personalizing combination therapies [START_REF] Barbolosi | Computational oncology -mathematical modelling of drug regimens for precision medicine[END_REF][START_REF] Imbs | Revisiting Bevacizumab + Cytotoxics Scheduling Using Mathematical Modeling: Proof of Concept Study in Experimental Non-Small Cell Lung Carcinoma[END_REF][START_REF] Serre | Mathematical Modeling of Cancer Immunotherapy and Its Synergy with Radiotherapy[END_REF] . Moreover, the current study needs to be extended to a larger number of patients.

Our study demonstrates the potential clinical utility of our computational platform as a personalized predictive tool of BM that could be helpful in the decision of WBRT versus localized treatment only. For instance for patient 1, at the time of detection of the first BM or at 40 months when only 6 BMs were visible, large amounts of invisible BMs were predicted by the model, which would advise in favor of WBRT. In contrast, a smaller metastatic burden was predicted in patient 2, thus advising for localized intervention only. Thus, our quantitative modeling approach could complement existing predictive metrics such as the brain metastasis velocity index [START_REF] Farris | Brain Metastasis Velocity: A Novel Prognostic Metric Predictive of Overall Survival and Freedom From Whole-Brain Radiation Therapy After Distant Brain Failure Following Upfront Radiosurgery Alone[END_REF][START_REF] Yamamoto | Validity of a Recently Proposed Prognostic Grading Index, Brain Metastasis Velocity, for Patients With Brain Metastasis Undergoing Multiple Radiosurgical Procedures[END_REF] , for patients undergoing stereotactic radiosurgery (SRS), i.e. patients presenting with a small (≤3-5) number of BMs. While the latter requires data at a second time point after a first SRS, our model could predict the invisible burden of BM already at the first SRS, and consequently inform the risk and timing of future relapse.

In order to translate our findings into a clinically usable tool, further methods need to be developed to calibrate the parameters of the model from data already available at diagnosis or at the first BM occurrence. By relying on a small number (3) of free parameters for the BM dynamics, our model represents a realistic candidate for such a purpose. To this respect, in addition to routine clinico-demographic features, molecular gene expression signatures [START_REF] Shedden | Gene expression-based survival prediction in lung adenocarcinoma: a multi-site, blinded validation study[END_REF] (for parameter μ, for instance) as well as radiomics features [START_REF] Aerts | Decoding tumour phenotype by noninvasive imaging using a quantitative radiomics approach[END_REF] (for parameter γ) might represent valuable resources.

Materials and Methods

patient data. The data used in this study concerned patients with non-small cell lung cancer (NSCLC) and were of two distinct natures: 1) population data of probability of BM as a function of PT size retrieved from the literature [START_REF] Mujoomdar | Clinical Predictors of Metastatic Disease to the Brain from Non-Small Cell Lung Carcinoma: Primary Tumor Size, Cell Type, and Lymph Node Metastases1[END_REF] and 2) longitudinal measurements of PT and BM diameters in two patients with NSCLC retrieved from imaging data (CT scans for lung lesions, MRI for brain tumors). Both patients had unresectable PT at diagnosis. The first patient (used for model development) was extracted from an EGFR mutated cohort from Institut Bergonié (Bordeaux, France). The second patient (used for model validation) had an EGFR wild type squamous cell lung carcinoma and was issued from routinely treated patients in the thoracic oncology service of the University Hospital of Marseille. The data comprised 10 PT sizes and 47 BM sizes (spanning 6 time points) for the first patient and 11 PT sizes and 16 BM sizes (spanning 4 time points) for second patient.

Data use from the Marseille patient was approved by a national ethics committee (CEPRO, Comité d'Evaluation des Protocoles de Recherche Observationnelle, reference number 2015-041), according to French law. Data use from the Marseille patient was approved by the institutional ethics committee of Institut Bergonié (Comité de Protection des Personnes Sud-Ouest et Outre Mer III) and IRB approval was obtained for use of the images.

The need for written informed consent of data collection was waived for this patient, in accordance with the related policy of Institut Bergonié. The study was performed in accordance with the Declaration of Helsinki, Good Clinical Practices, and local ethical and legal requirements.

Mathematical modeling of primary tumor growth and metastatic development. Primary tumor

growth. The pre-diagnosis natural history of the primary tumor size S p (t) for times < t T d (diagnosis time)expressed in number of cells in the model computations -was assumed to follow the Gompertz growth model [START_REF] Norton | A Gompertzian model of human breast cancer growth[END_REF][START_REF] Benzekry | Classical mathematical models for description and prediction of experimental tumor growth[END_REF] , i.e.

α β            = - = = ⇒ = ∀ < α β -β - dS dt S S S t S t e t T ( l n( )) ( 0) 1 ( ) ,
p p p t 0, (1 )
where time t = 0 corresponds to the first cancer cell, parameter α 0,p is the specific growth rate (i.e. S dS dt 1 p p ) at this time and β p is the exponential rate of decrease of the specific growth rate. Conversions from diameter measurements to number of cells were performed assuming spherical shape and the classical assumption 1 mm 3 = 10 6 cells [START_REF] Spratt | Rates of growth of human solid neoplasms: Part I[END_REF] . After treatment start (T d ), the primary tumor size was assumed to follow a tumor growth inhibition model [START_REF] Claret | Model-based prediction of phase III overall survival in colorectal cancer on the basis of phase II tumor dynamics[END_REF] consisting of: 1) exponential growth (rate α 1 ), 2) log-kill effect of the therapy (efficacy parameter κ) [START_REF] Benzekry | Metronomic reloaded: Theoretical models bringing chemotherapy into the era of precision medicine[END_REF] and 3) exponential decrease of the treatment effect due to resistance, with half-life t res . The equation is:

α κ = - ∀ ≥ . - - dS dt S e S t T , ( 1
)
p p t T t p d 1 ln2 d res
Metastatic development. Basic and different growth models. The general modeling framework we employed was derived from the work of Iwata et al. [START_REF] Iwata | A Dynamical Model for the Growth and Size Distribution of Multiple Metastatic Tumors[END_REF] . It consists in modeling the population of metastases by means of a size-structured density ρ t s ( , ), of use to distinguish between visible and invisible tumors. Metastatic development of the disease is reduced to two main phases: dissemination and colonization [START_REF] Chaffer | A perspective on cancer cell metastasis[END_REF] . The multiple steps of the metastatic cascade [START_REF] Gupta | Cancer metastasis: building a framework[END_REF] are aggregated into a dissemination rate with expression:

μ = γ d S S ( ) , (2)
p p which corresponds to the number of successfully born BM per unit of time. In this expression, the geometric parameter γ corresponds to the intra tumor repartition of the cells susceptible to yield metastasis and μ is the per day per cell probability for such cells to overcome all the steps of the metastatic cascade (acquisition of metastatic-specific mutations, epithelial-to-mesenchymal transition, invasion of the surrounding tissue, intravasation, survival in transit, extravasation and survival in the brain). For γ = 1 all cells in the PT have equal probability to give a BM whereas a value of γ = 0 indicates a constant pool of cells having metastatic ability (cancer stem cells). Intermediate values 0 < γ < 1 can be interpreted as the geometric disposition of the metastatic-able cells, including the surface of the tumor (γ = 2/3) or a fractional dimension linked to the fractal nature of the tumor vasculature [START_REF] Baish | Cancer, angiogenesis and fractals[END_REF] . Assuming further that the growth of the metastasis follows a gompertzian growth rate

α β = - g s s s ( ) ( l n( ))
0 with growth parameters α 0 and β possibly equal (basic model (A)) or distinct (different growth model (B)) compared to the PT ones, the density ρ satisfies the following transport partial differential equation [START_REF] Iwata | A Dynamical Model for the Growth and Size Distribution of Multiple Metastatic Tumors[END_REF] : where S 0 is the size of a BM at birth (here assumed to be one cell). From the solution to this equation, the main quantity of interest for comparison to the empirical data is the number of metastasis larger than a given size s (cumulative size distribution):

ρ ρ ρ ρ            ∂ +∂ = ∈ +∞ ∈ +∞ = ∈ +∞ = ∈ +∞ t s t s g s t s S g S t S d S t t s s S ( , ) ( ( , ) ( )) 0 (0, ), ( , ) ( ) ( , ) ( ( )) ( 0, ) (0, ) 0 (, ) , (3)
∫ ρ = ′ ′. +∞ f t s ts ds ( , ) (, )
s

The total number of metastases -denoted N(t) -is obtained by using s = V 0 above and its expression can be directly computed without solving the entire problem (3) as it is given by:

∫ = ′ ′. N t d S t dt ( ) ( ( )) (4) t p 0
Using the method of characteristics, one can derive the following relationship between N and f:

= - f t s N t t s ( , ) ( ( )), (5)
where t(s) is the time for a tumor growing at rate g to reach the size s. In the case of Gompertz growth one has:

β β α = -       -       . t s s ( ) 1 ln 1 l n( ) 0
Of particular interest is the number of visible BMs f t S ( , ) vis with S vis the minimal visible size at CT scan taken here to be 5 mm in diameter.

Secondary dissemination. In the previous model formulation, all BMs were assumed to have been seeded by the primary tumor. When BMs are also allowed to spread metastases themselves, this results in a second term in the boundary condition of (3) and the model becomes [START_REF] Iwata | A Dynamical Model for the Growth and Size Distribution of Multiple Metastatic Tumors[END_REF] : In this case, formula (4) is not valid anymore, which complicates substantially the computation of the cumulative size distribution. A dedicated scheme based on the method of characteristics was employed [START_REF] Benzekry | Mathematical and numerical analysis of a model for anti-angiogenic therapy in metastatic cancers[END_REF] . The cumulative size distribution is then given by:

∫ ρ ρ ρ ρ ρ              ∂ +∂ = ∈ +∞ ∈ +∞ = + ∈ +∞ = ∈
τ = - = - -. f t s N t t s N t t s ( , ) ( ( )) ( ( ) )
dorm dorm Discrete versions of the models. While continuous versions of the models were used for fitting the model to the data because they allow computations to be tractable, discrete versions were implemented for forward calculations, because of the small number of BMs. Briefly, in a discrete simulation, the appearance time of the i-th BM T i was defined by

= > ≥ T t N t i inf{ 0; ( ) }, i
with N(t) the number of metastases as defined above. The size of the i-th BM S i (t) was then defined, for t > T i , by:

         = = = . dS dt g S S t T S ( ) ( ) i i i i 0
This discrete version corresponds to a Poisson process where the number of events (metastatic births)  t ( ) is replaced by its expectation:

 = N t t ( ) [ ( )]  .
For further links between stochastic (Poisson process) and continuous versions of the Iwata model, the reader is referred to [START_REF] Gomez | Stochastic and Deterministic Models for the Metastatic Emission Process: Formalisms and Crosslinks[END_REF] . Simulations reported in Fig. 5B,C and 6 were generated using discrete versions of the models.

Models fitting and estimation of the parameters.

Pre-clinical primary tumor growth. To parameterize the Gompertz function defining the PT growth, two parameters needed to be defined (α 0,p and β p ). In the absence of longitudinal measurements of the PT size without treatment, these two parameters were determined from two considerations: 1) the maximal reachable size (carrying capacity, equal to ) of a human tumor is 10 12 cells [START_REF] Retsky | Computer simulation of a breast cancer metastasis model[END_REF][START_REF] Klein | Parallel progression of primary tumours and metastases[END_REF] and 2) the histology-dependent value of the doubling time at diagnosis, retrieved from a meta-analysis of published literature about the natural growth of lung PTs (see Table S1, extended from [START_REF] Friberg | On the growth rates of human malignant tumors: implications for medical decision making[END_REF][START_REF] Detterbeck | Turning gray: the natural history of lung cancer over time[END_REF] ). The latter yielded values of 201 days for an adenocarcinoma and 104 days for an undifferentiated carcinoma. For the Gompertz model, the doubling time is size-dependent and its value for the PT diagnosis size DT(S d ) is given by:

β β β α = -       -        . DT S S S ( ) 1 ln ln(2 ) ln d p p d p d p 0,
Using the formula linking α 0,p and β to the carrying capacity, this nonlinear equation was numerically solved.

Population level: probability of BM apparition. To fit the data from [START_REF] Mujoomdar | Clinical Predictors of Metastatic Disease to the Brain from Non-Small Cell Lung Carcinoma: Primary Tumor Size, Cell Type, and Lymph Node Metastases1[END_REF] describing the probability of BM in a population of lung adenocarcinoma patients, we employed a previously described methodology [START_REF] Benzekry | Modeling Spontaneous Metastasis following Surgery: An In Vivo-In Silico Approach[END_REF][START_REF] Barbolosi | Modélisation du risque d' évolution métastatique chez les patients supposés avoir une maladie localisée[END_REF] . Briefly, we considered that the probability of developing BM after diagnosis was the probability of having already BM at diagnosis, i.e.  > N T ( ( ) 1) d , with N defined by equation ( 4). Of note, since this quantity only depends on the dissemination rate, this probability is equal for all the model scenarios introduced above. We fixed the value of the PT growth parameters as described above from the cohort histology and set γ = 1 as the simplest dissemination model. The inter-individual variability was then minimally modeled as resulting from a lognormal population distribution of the parameter

D  ∫ µ µ µ µ µ =      > |      . σ σ
The best-fit of these probabilities -evaluated by Monte Carlo simulations -to the empirical data was then determined by least squares minimization performed using the function fminsearch of Matlab (Nelder-Mead algorithm) [START_REF] Mathworks | Matlab with statistics and optimization toolboxes[END_REF] .

Individual level: description of longitudinal data of number and size of BM growth. Maximum likelihood estimation. Due to the discrete nature of the data at the individual level (diameters of a small number of BMs at discrete time points), a direct comparison between the size distribution ρ solution of the problem (3) was not possible. Instead, we compared the data to the model by means of the cumulative size distribution. Denoting by t j the observation times, x j i the sorted BM sizes at time t j and y j i the number of metastases larger than size x j i at time t j , we considered the following nonlinear regression problem: regroups all the parameters of the model. Note that, except for the "secondary dissemination" model, all models can be viewed as submodels of a general model including all these parameters (the "basic model" consisting of the case α α β β = = = t , , 0 p p 0, 0 0 and τ = 0, for instance). Classical maximum likelihood estimation then leads to the following estimate: where θ Ĵ( ) is the jacobian matrix of the model (with respect to the parameter vector θ) at all time points and all sizes, evaluated at the optimal parameter θ ˆ and σ ˆ is the a posteriori estimate of σ given by

θ σε ε = + ∼ y f t
σ θ = - ˆN P l 1 ( )
with N the total number of data points and P the number of free parameters. Using the standard errors as an identifiability metric, we repeatedly observed a lack of identifiability of parameters μ and γ when fitted together. Indeed, standard errors for μ and γ were larger than 200% when fitting the basic model to the data. Further investigation of the shape of the objective function confirmed this lack of idenfiability (Fig. S7). To address this issue, we only considered a finite set of relevant possible values for γ and only optimized the value of μ. These values were . . 1 2 . When more parameters were let free in the model (delay t d or dormancy period τ ), we generated 4 × 4 parameters grids for initial conditions (with ∈ t {0, 500, 1700, 2000} d or τ ∈ {0, 50, 100, 350}). Of these 16 optimization problems, the one with the minimal value of l at convergence was selected.

Figure 1 .

 1 Figure 1. Probability of brain metastasis as a function of the primary tumor size. Fit of the mechanistic model for probability of BM to data from Mujoomdar et al. 50 . (A) Adenocarcinoma data (n = 136). Inferred values for the distribution of μ: μ = . × - --1 39 10 cell day pop

Figure 2 .

 2 Figure 2. Data of primary tumor and brain metastases in a patient with non-small cell lung cancer. (A) Postdiagnosis kinetics of the primary tumor largest diameter, measured on follow-up computerized tomography images (inlet). This EGFR mutated patient was treated first with a tyrosine kinase inhibitor (erlotinib) then with several rounds of additional systemic treatments upon relapse (cytotoxic chemotherapy (docetaxel), rechallenge with erlotinib and a second cytotoxic agent (gemcitabine)), as indicated by arrows and dashed lines in the figure. The solid line corresponds to the model fit during treatment (tumor growth inhibition (TGI) model). (B) Growth kinetics of the brain metastases. The solid line corresponds to Gompertz growth predictions based on parameters estimated from the primary tumor size at diagnosis and histological type, using as initial condition the first observation of each BM. (C) Time course of the apparition of visible metastases. (D) Cerebral magnetic resonance image showing two brain metastases 48 months post-diagnosis (other brain lesions not visible on this slice). EGFR = Epithelial Growth Factor Receptor, TKI = Tyrosine Kinase Inhibitor, mTKI = maintenance Tyrosine Kinase Inhibitor, CT = (cytotoxic) chemotherapy. Time is in months from diagnosis.

Figure 3 .

 3 Figure 3. Schematic of the models and investigated hypotheses. Several biological hypotheses about metastatic dynamics were formalized as mathematical models and tested against the data. These included: same (A) or different (B) growth rates between primary and secondary tumors; only primary dissemination from the primary tumor or (C) secondary dissemination from brain metastases themselves; (D) possibility of a delay between cancer initiation and onset of metastatic ability; and (E) the possibility of a dormancy lag time between metastatic spread and growth initiation.

Figure 4 .

 4 Figure 4. Fit of the dormancy model. (A) Time course of the visible brain metastases (BMs) size distributions during follow-up. Comparison between model calibration and data. (B) Time course of the number of visible BMs. T 1 = time of diagnosis. T 2 = detection time of the first brain metastasis. (C) Comparison of the BM size distribution between the model fit and the data at last follow-up (48 months post-diagnosis). Time is in months from diagnosis.

Figure 5 .

 5 Figure 5. Model simulated predictions of the pre-and post-diagnosis time course of the primary and cerebral disease. (A) Model-inferred growth kinetics of the primary tumor (in blue) and the brain metastases (in red), compared with data (circles). T 0 = time of first cancer cell, T 1 = time of diagnosis, T 2 = detection time of the first brain metastasis. Only brain metastases that will become visible are shown. (B) Predicted size distribution of the brain metastases at diagnosis. (C) Predicted size distribution of the brain metastases at the time of clinical occurrence of the first one. Simulations performed with the discrete version of the model. PT = primary tumor, BM = brain metastasis. Time is in months from diagnosis.

2

 2 Parameters identifiability. Standard errors can be computed from this statistics' covariance matrix, given by 80 :

Table 1 .

 1 Minimal value of the objective function obtained when fitting each of the models to the data. The objective function was the sum of squared residuals (see expression(7)).

	Model	Patient 1	Patient 2
	Basic	5.51	2.53
	Secondary	5.43	2.3
	Delay	5.23	1.53
	Dormancy	4.93	1.71
	Diff. growth	4.95	1.79

Table 2 .

 2 Patient-specific parameters. Values of the parameters from either the data (Data), the model assumptions

	Parameter	Meaning	Unit	Patient 1	Patient 2	Determination
	S d	PT size at diagnosis	mm	36.0	53.7	Data
	Histology			Adenocarcinoma	Squamous cell carcinoma	Data
	T d	Age of the PT at diagnosis	year	5.3	2.1	Predicted
	Visibility threshold		mm	3	3	Fixed
	S 0,p	Initial size of the PT	cell	1	1	Fixed
	S 0	Initial size of a BM	cell	1	1	Fixed
	α p	Proliferation rate at one cell (PT)	day -1	0.0284	0.0858	Computed
	β p	Exponential decrease of growth rate (PT)	day -1	1.03 × 10 -3	3.10 × 10 -3	Computed
	DT PT	PT doubling time at 3 cm	day	173	57	Predicted
	α 1	PT growth rate during relapse	day -1	5.72 × 10 -4 (23.6)	0.0101 (5.28 × 10 -3 )	Fit PT
	κ	PT log-kill effect	day -1	4.46 × 10 -3 (20.6)	0.0439 (1803)	Fit PT
	t res	Half-life of treatment response	day	149 (1.05 × 10 -3 )	35.4 (0.998)	Fit PT
	α 0	Proliferation rate at one cell (BM)	day -1	0.0284	0.0858	Computed
	β	Exponential decrease of growth rate (BM)	day -1	1.03 × 10 -3	3.10 × 10 -3	Computed
	τ	Dormancy duration	day	133 (4.23)	171 (33.3)	Fit BM
	DT BM	BM doubling time at 3 cm	day	173	57	Predicted
	μ	Cellular metastatic potential	cell -1 × day -1 2.00 × 10 -12 (2.83)	1.02 × 10 -12 (71.9)	Fit BM
	γ	Fractal (Hausdorff) dimension of the vasculature (x1/3)	-	1	1	Fit BM*
	d	Dissemination rate at 3 cm	day -1	0.0282	0.0144	Predicted

(Fixed)

, estimated from direct computations from diagnosis data (Computed), fit of the TGI model to the PT data (Fit PT) or fit of the metastatic dormancy model to the brain metastases data (Fit BM), or computed from these estimations (Predicted). Values in parenthesis for fitted parameters correspond to the relative standard errors expressed in percent, i.e. × 100 se est

.

  Delay model. Consideration of a delay t 0 before onset of metastatic dissemination in the model can be taken into account by remarking that Dormancy model. For inclusion of dormancy in the model -defined as a period of duration τ during which a newborn metastasis remains at size S 0 -the time to reach any given size >

		N delay	t ( )	=	∫ t	t 0	′ ′ = d S t dt ( ( ) p	N t ()	-	N t ( ) 0	.
	Thus
	f delay	t s ( , )	=	N delay	(	-t t s ( ))	=	-N t t s ( ( ))	-	N t ( ) 0	.
												s S 0 becomes	t dorm ( ) s	=	t s ( )	+	τ

Scientific RepoRtS | (2019) 9:13018 | https://doi.org/10.1038/s41598-019-49407-3

Scientific RepoRtS |

Acknowledgements

The authors would like to thank Pr A. Iliadis for valuable comments and advice.

Data Availability

The datasets analyzed during the current study are available from the corresponding author on reasonable request.

Author Contributions

Additional Information

Supplementary information accompanies this paper at https://doi.org/10.1038/s41598-019-49407-3.

Competing Interests:

The authors declare no competing interests.

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Mario Barbatti
email: mario.barbatti@univ-amu.fr

Kakali Sen

Effects of different initial condition samplings on photodynamics and spectrum of pyrrole

Semi-classical simulations of spectrum and dynamics of complex molecules require statistical sampling of coordinates and momenta. The effects of using thermal and quantum samplings are analyzed taking pyrrole as a test case. It is shown that there are significant differences in the results obtained with each of these two approaches. Overall, quantum sampling based on a Wigner distribution renders superior results, comparing well to the experiments. Dynamics simulations based on surface hopping and ADC(2) reveal that pyrrole internal conversion to the ground state occurs not only through H-elimination path, but also through ring-distortion paths, which have been systematically neglected by diverse experimental setups. The analysis of the reaction paths also shows that the ionization potential varies by more than 5 eV between ionization of the excited state at the Franck-Condon region and at the intersections with the ground state. This feature may have major implications for time-resolved photoelectron spectroscopy.

Introduction

In diverse types of simulations, one needs to sample nuclear conformations of a molecule considering the kinetic energy deposited in the internal vibrational modes. This is the case, for instance, of the sampling of the phase space (geometries and velocities) to initiate excited-state dynamics simulations or the sampling of the configurational space (geometries only) for spectrum simulations. There are two basic approaches to carry out these samplings: one may either sample the space by considering that the molecule is at its zero-point level (we will refer to it as the quantum or Q sampling); or, alternatively, one may sample the space supposing that the molecule is in thermal equilibrium with the environment (thermal or T sampling).

Although we can find examples of both approaches in the literature, they are far from equivalents. While in the Q sampling each vibrational degree of freedom j has energy / 2 Top: Illustration of the classical (T = 300 K) and quantum energies and amplitudes in a onedimensional harmonic potential. Bottom: Ratio between the classical and quantum amplitudes as a function of the wavenumber (Eq. ( 1)). The proportions in the upper figure correspond to 1600 cm -1 . This difference of energy per degree of freedom reflects on the vibrational amplitude.

Considering that the amplitude of a harmonic oscillator with energy  is  

1/ 2 2 2 / a Q   
, the ratio between the amplitudes in the thermal and quantum samplings is

1/2 2 aT aQ j Q kT Q            . (1)
In particular, Eq. (1) implies that the deviation between the thermal and quantum amplitudes tends to increase for larger wavenumbers (Figure 1-bottom). In this work, we take pyrrole (Scheme 1) as a test case for investigating the differences between the Q and T samplings. Pyrrole is a small molecule with ultrafast photodynamics. 1 This makes it ideal for benchmarking [2][3][4][5][6][7][8][9][10] as well as for concept and method developments. [11][12][13][14] Moreover, after high-energy excitation into the bright * state, pyrrole dynamics branches in several distinct reaction paths, 15 facilitating the comparison and evaluation of the sampling effects. The ultrafast dynamics of photoexcited pyrrole has also been the subject of a large number of experimental (time 1,[16][17][18][19][20][21][22] and frequency [23][24][25] domains) and theoretical investigations (statics [26][27][28] and dynamics [29][30][31][32][33][34][35] ).

We have ourselves taken pyrrole before to investigate dynamics branching, 36,37 to test methodology, 38 and to understand the role of the conical intersection branching space on the photoproducts. 39 Here, we will use it for simulating spectrum and excited-sate dynamics starting from Q and T samplings, aiming at understanding their impact on the simulations.

The difference between the Q and T samplings has not gone unnoticed, although it has been largely neglected in the field. As far as we know, the only work comparing them in the context of photochemistry is in the paper by Klaffki et al., 40 where they analyzed the results of excited-state dynamics simulations for a protonated Schiff base (PSB3) based on initial conditions generated with both quantum and thermal samplings. They observed a slight trend towards faster reactions when using Q sampling, although the differences were not really statistically significant. They finally prescribed the use of thermal sampling as the preferred methodology, as the Q sampling would provide "artificially high kinetic energies" and lower numerical stability.

We tend to disagree with this analysis. The high kinetic energies of the Q sampling are not artificial, but correspond to the kinetic energies expected for a quantum system. In conceptual terms, the Q sampling through a Wigner distribution should be the proper way of mapping quantum densities on the classical phase space. 41 From a practical perspective, we have observed before 42 that while spectra simulated with Q sampling rendered satisfactory results, those based on T sampling were too narrow, needing much larger nominal temperatures to reach proper band width description. We have even proposed 43,44 a hybrid method to sample initial conditions for QM/MM surface hopping dynamics based on Q sampling for the QM and T sampling for the MM regions. The present results for spectrum and dynamics will add to these previous evidences of the superiority of Q over T samplings.

The main contribution of the present work is at providing a more complete and systematic investigation of the effects of sampling on electronic spectrum and nonadiabatic excited-state dynamics. Although the effects of sampling will most likely depend on each particular molecule, we try to figure out features that may reflect general trends of the different sampling methods. We also aim at establishing protocols for each step of the simulations for evaluating the sampling effects, which may be used in other cases as well. As a side-product of this methodological investigation, we also report a new set of excited-state dynamics for pyrrole in the gas phase.

Computational details

Electronic Structure and Dynamics

Excited states were computed with the algebraic diagrammatic construction to the second order [ADC(2)], 45,46 using the resolution-of-the-identity (RI) approximation 47 and frozen 1s orbitals. The reference ground state was computed at the second-order Møller-Plesset Perturbation Theory (MP2). 48 Calculations were done with the aug-cc-pVDZ basis set assigned to for N and C atoms, and with the cc-pVDZ assigned to H atoms. 49 The absorption spectrum and initial conditions for dynamics were simulated with the nuclear ensemble method 50 using different types of sampling as discussed below. On-the-fly dynamics simulations were carried out in the excited states computed with ADC(2). Nonadiabatic effects were taken into account by the surface hopping approach. 51 Classical equations were integrated with 0.5 fs time step, while quantum equations were integrated with 0.025 fs, using interpolated quantities between classical steps. The maximum simulation time was 300 fs.

Hopping probabilities between excited states were computed with the fewest switches approach 52 including decoherence corrections ( = 0.1 Hartree). 53 Nonadiabatic couplings with ADC(2) were computed by finite differences with the method discussed in Ref. 54 , which is based on the Hammes-Schiffer-Tully approach. 55 As a single-reference and low-order method, ADC(2) cannot provide reliable nonadiabatic couplings for crossings with the ground state. 54 For this reason, when a trajectory reached an S 1 -S 0 energy gap smaller than 0.15 eV before the maximum simulation time, it was stopped. This crossing time was taken as an estimate of the internal conversion time to the ground state. ADC(2) calculations were carried out with the Turbomole program. 56 The spectrum and dynamics simulations were performed with Newton-X 57,[START_REF] Barbatti | NEWTON-X: a package for Newtonian dynamics close to the crossing seam[END_REF] interfaced with Turbomole. Surface intersections (within 0.01 eV) were optimized with the penalty Lagrange multiplier technique ( = 0.02 hartree) implemented in the Ciopt program, [START_REF] Levine | [END_REF] which we have adapted to work with ADC(2) and Turbomole. Cremer-Pople parameters were computed with Platon. 60

Ensemble Sampling

Quantum sampling

The ground-state phase space of a molecule with N at atoms can be conveniently sampled by a harmonic-oscillator Wigner distribution    

N i i W q N i qi pi q p P                           q p  , ( 2
)
where 2 2 2 2 qi i i i i pi           (3)
In these equations, q i and p i are the coordinate and momentum for each normal mode i with reduced mass  i and angular frequency  i . N p = 2000 nuclear geometries and momenta were sampled using Eq. ( 2) and converted back to Cartesian coordinates. The normal coordinates for the electronic ground state minimum were computed at MP2/(aug)-cc-pVDZ. This set will be referred to as the Q sampling throughout the paper.

Thermal sampling

A Born-Openheimmer (BO) ground-state dynamics trajectory was run for 10 ps with 1 fs time step using MP2/(aug)-cc-pVDZ. Canonical ensemble (300 K) was enforced via the Andersen thermostat 61 with 0.2 fs -1 collision frequency. A total of 20 snapshots were randomly selected between 5 and 10 ps of this trajectory to spawn 20 new trajectories. These trajectories run for 25 ps with the same parameters as the original one. Combined, these new trajectories corresponded to 500 ps ground-state dynamics. N p = 2000 snapshots (geometries and velocities) were randomly selected from this set and used to compute the absorption spectrum and also to sample the initial conditions for the excited state dynamics. This set will be referred to as the T 300 sampling. The same procedure was also adopted to build a second thermal ensemble at T = 912 K (T 912 sampling).

The thermal sampling was also prepared based on molecular mechanics for pyrrole in the gas phase using AMBER 11 suite of programs. 62 The charges and parameters of pyrrole were taken from Charmm General Force Field for drug like molecules. 63,64 Pyrrole was first minimized, then slowly heated from 0 to 300 K for 30 ps, and finally equilibrated for 20 ns at 300 K using Langevin thermostat with a collision frequency of 1.0 ps -1 . After heating, the production run of 500 ps was performed without electrostatic cutoffs and using the Andersen thermostat with 0.01 fs -1 collision frequency. The production run was done with 1 fs time step for 500 ps. 5000 snapshots where used for analysis of the ensemble, which will be referred to as the MM 300 T sampling.

Spectrum Simulations

Ensembles computed with quantum and thermal samplings were used to simulate the absorption cross section based on the nuclear ensemble approach. For each of the N P geometries (R l ) in the ensemble, transition energies (E 0n ) and oscillator strengths (f 0n ) for N fs = 25 excited states were computed at ADC(2)/(aug)-cc-pVDZ level. The absorption cross section is given by 50

          2 0 0 0 0 1 , 2 fs p N N n l n l n l n l p e E E f gE E mc E N           R R R  , ( 4
)
where E is the photon energy of the radiation,  0 is the vacuum permittivity, c is the speed of light, and e and m are the electron charge and mass. g is a normalized Gaussian line shape

          2 0 0 1/ 2 2 2 1 , e x p 2 /2 2 /2 n n E E g E E                  , (5)
centered at E 0n and with width  = 0.05 eV.

Results and Discussion

Sampling characterization

We start by analyzing the energy and geometry distributions arising from each type of sampling. In total, we produced four ensembles: one from quantum sampling (Q) based on the harmonic-oscillator Wigner distribution, and three from thermal samplings. Among the thermal sampling, two ensembles were built at 300 K (T 300 from a BO trajectory and MM 300 T from a MM trajectory) and the last one at 912 K (T 912 also from a BO trajectory), as explained later.

The distribution of kinetic and total ground-state energies for the quantum (Q) and thermal samplings at 300 K (T 300 and MM We may also build the thermal ensemble with a temperature T ZP set to correspond to the zero point energy E ZP :

3 2 2 ZP at ZP E N kT  . ( 6
)
For pyrrole, T ZP = 912 K. In this case, the mean value of the T 912 kinetic-energy distribution is displaced to E ZP /2 by construction.

Figure 2 also shows that the width of the energy distributions (kinetic and total) is much narrower in the T 300 and MM 300 T samplings than in the Q sampling. To understand this difference, notice that the total energy variance in the quantum sampling based on a harmonic-oscillator Wigner

distribution is   2 2 2 2 3 6 3 6 2 2 1 1 1 2 2 2 at at N N pi i i qi EQ i i i i                                        , (7)
yielding a quantum distribution with width 2 1 .12 eV EQ   for pyrrole. In the case of the thermal sampling, supposing that the molecular vibration can be described as Boltzmann-distributed independent harmonic oscillators, the energy variance is

  2 2 2 3 6 ET at N k T    , (8)
Thus, the width of the classical distribution is 2 0.25 eV

ET

 

, much narrower than the width of the quantum distribution. Once more, we may check the performance of the T 912 sampling. As shown in Figure 2, the result is in good agreement with the Q sampling.

In the literature, it is also common to sample the geometries using a harmonic oscillator distribution, but then constraining velocity magnitude such that the total energy is the E ZP . This procedure, used for instance in Ref. 40 , results in a total energy distribution for the quantum sampling with zero variance.

In terms of coordinates, the deviations between the distributions of the T-and Q-samplings depend on the specific vibrational frequencies, as expressed in Eq. ( 1) (see also Figure 1-bottom). In Figure 3, we plot the distributions of three internal coordinates in the ground state for the thermal and quantum samplings.

The first internal coordinate (Figure 3-top) is the NH distance. In the Q sampling, the main normal mode contributing to this coordinate is the NH stretching vibrating at 3683 cm -1 . As expected (Eq. ( 1)), this distribution is much wider than that obtained for both T samplings at 300 K. The deviations between T 300 and Q are much smaller for the CN-distance distributions (~1200 cm -1 , see Figure 3-middle), and almost negligible for the CNCC-dihedral distributions (~650 cm -1 ). The overheated T 912 sampling produces distributions that are equivalent to those of the Q sampling for middle-range wavenumbers, like in the case of the CN-distance distribution (Figure 3middle). However, it produces too narrow distributions for high-frequency cases (Figure 3-top) and too broad distributions for low-frequency cases (Figure 3-bottom). We will see in the next section that this unbalanced distribution of energy among the modes leads to artifacts in the spectrum simulations.

Once more, we observe that there is no significant difference between using a BO or MM trajectories to sample the geometries in T 300 and MM 300 T . For this reason, in the next section we will proceed only with BO samplings, T 300 and T 912 , to simulate the spectrum.

Absorption Spectrum

We have used the nuclear ensemble approach to simulate the absorption spectrum of pyrrole based on ADC(2) excitation energies and oscillator strengths for the Q, T 300 , and T 912 samplings.

The nuclear ensemble approach is possibly the simplest approach to simulate electronic spectrum including band shapes. In this method, energies and transition probabilities (proportional to the oscillator strengths) are computed for each geometry in an ensemble representing the initial state distribution, and then summed over (see Eq. ( 4)) to provide the spectrum. 50,65 The nuclear ensemble approach is intuitive, straightforward to implement, and delivers absolute widths and heights even for dark vibronic bands. However, as it does not incorporate any information on the final states wave functions, it completely misses vibrationally resolved features, and tends to predict too symmetric band shapes. (This and more involved approaches for spectrum simulations are discussed, for instance, in Ref. 66 .) The nuclear ensemble approach is also very handy as the initial step for excitedstate dynamics simulations. The same set of geometries used for the spectrum, together with their corresponding momentum, can be used as the initial conditions to start semiclassical trajectories.

The vertical excitations of pyrrole computed with ADC(2) are shown in Table 1. The rootmean-square deviation (RMSD) to the vertical energies computed with CC3 (Ref. 10 ) is 0.2 eV. Most of vertical transitions are slightly overestimated compared to CC3, with exception of the 1 1 A 1 state, the only state in the series with a dominant valence character. The triple- basis marginally stabilizes the double- results by 0.1 eV (RMSD).

-Singlet vertical excitations of pyrrole at ADC(2)/(aug)-cc-pVDZ. Velocity-gauge oscillator strengths. Table 1 The experimental 25 oscillator strength between 5.5 and 6.5 eV is given in the parenthesis. Results at ADC(2)/aug-cc-pVTZ are given within brackets.

ADC(2)

CC3 10 State E (eV) f E (eV) Assignment Expt. (eV) 25 The experimental spectrum 25 of pyrrole below 6.5 eV is dominated by a broad and intense transition into the 1 1 B 2 state, which is peaked at 5.95 eV (Figure 4). A narrow intense peak due to the origin of the 2 1 B 1 transition features at 5.86 eV. A recent assignment of pyrrole spectrum based on quantum-dynamics simulations is discussed in Ref. 33 .

The band maximum is usually taken as the reference value to check computed vertical excitations. Nevertheless, to properly compare theory and experiment, especially in the case of a broad band as that of the 1 1 B 2 transition, we should take into account that the experimental band maximum does not correspond to the vertical excitation. We can estimate the deviation between the band maximum and the vertical excitation from the simulated spectrum. The simulated band maximum (Q sampling, Figure 4) is at 6.08 eV. The vertical excitation into the 1 1 B 2 state is at 6.25 eV (Table 1), 0.17 eV higher than the maximum. Assuming that the same shift applies to the experimental data, whose band maximum is at 5.95 eV, the "experimental vertical excitation" may be taken as 6.12 eV. The agreement between the ADC(2) vertical excitation into the 1 1 B 2 state (6.25 eV for both basis sets) and the "experimental vertical excitation" (6.12 eV) is excellent. Photoabsorption cross section computed with the quantum and thermal (T = 300 and 912 K) samplings compared to the experimental results. 25 The absorption spectrum computed with the thermal and quantum samplings based on ADC(2)/(aug)-cc-pVDZ states are shown in Figure 4. The experimental result from Ref. 25 is also shown in the figure. At 300 K, the T 300 sampling produces a spectrum that has qualitatively wrong band shape, when compared to the experimental result. The spectrum based on the Q sampling provides a much better result for the band shape, especially in the low-energy edge.

The failure of the T 300 sampling in the description of the low-energy edge is related to the geometry of the bright * minimum. This minimum, the S 4 state with ADC(2), is reached by shortening the C3-C4 bond (by 0.06 Å at ADC(2) level) and stretching the C2-C3 and C5-C4 bonds (both by 0.07 Å). These bond deformations, which lead to small energy gaps of about 5.7 eV, are not well described by the T 300 sampling.

The spectrum computed with the T 912 sampling is also shown in Figure 4. The result is significantly better than at 300 K, with the band shape much closer to the experimental one. There is, however, an overestimation of the absorption cross section in the low-energy edge. This is caused by the excess of energy deposited in ring-deformation modes, which leads to artificially too broad distributions as shown in Figure 3-bottom.

Dynamics simulations

To simulate semiclassical dynamics, we need an ensemble of initial geometries and velocities to start each trajectory. Rigorously, surface hopping should be based on a double ensemble strategy: 51 trajectories are started from different geometries and velocities representing the initial distribution of the molecule in the phase space; and for each of these points several trajectories should also be started to sample the random nonadiabatic processes. (Different from classical dynamics, in surface hopping, trajectories initiated with the same conditions may have different fates due to the stochastic nature of the hops.) Commonly, to reduce the computational costs, this double-ensemble approach is replaced by a single ensemble, where only a single trajectory is started from each phase-space point. This is the approach we adopt here as well.

Two sets of excited-state dynamics simulations were run, one starting from initial conditions generated with the Q sampling and another with the T 300 sampling. In both cases, initial conditions were restricted to the 6.1 ± 0.1 eV energy window, exciting the →3p x state near the center of the first absorption band. 100 trajectories were computed for each sampling. The number of trajectories starting in each adiabatic state (see Table 2) was taken as proportional to the transition probabilities into each state within the energy window. Thus, in the Q sampling, for instance, 19 trajectories were started in S 3 , 30 in S 4 , 36 in S 5 , and 15 in S 6 . This procedure ensures that the initial conditions in both samplings are equivalent in terms of excitation energy and diabatic character of the initial state.

-Distribution of ensemble points among the electronic states according to their oscillator strengths, in The general time-evolution of the trajectories is illustrated in Figure 5. Qualitatively, the results for Q and T 300 samplings are similar. The occupation of the initially populated states (S 3 to S 6 ) quickly reduces; first populating S 2 and S 1 with time constant  1 and then S 0 with  0 :

      0 12 1/ 1/ 3 7 1 2 0 ... S S S S S           . ( 9
)
With such a partition (which includes S 7 in the first group to account for the up-hops from S 6 ), we can fit the [S 1 +S 2 ] occupation analogously to a unimolecular reaction with an intermediary step

  0 1 2 12 0 12 0 exp exp S S t t F t                                . ( 10
)
The result of this fitting is given in Table 3. Clearly, the initial transfer into S 2 and S 1 is independent of the sampling, and about 27 fs. Nevertheless, the transfer from S 2 and S 1 into S 0 with Q sampling is significantly faster than that with the T 300 sampling. While with the Q sampling S 0 is populated within approximately  0 = 90 fs, with the T 300 sampling it takes about 148 fs.

-Characterization of the T-sampled (300 K) and Q-sampled dynamics in the 6.1±0.1 eV excitation Table 3 window. Time constants:  12 for exponential transfer from the high excited states into S 1 and S 2 ;  1 for the S 2 transfer to S 1 ;  0 for the transfer from S 1 to S 0 ;  H for S 1 to S 0 transfer via H elimination;  R for S 1 to S 0 transfer via ring distortion;  * for crossing from 3s to * (H fragments formation). There is also an important difference in the activation of the different reaction paths. Pyrrole internal conversion may occur either through a */ground-sate intersection induced by the NH stretching (H-elimination pathway) or through */ground-state intersection induced by simultaneous ring opening and puckering (ring-distortion pathway). 37 Dynamics starting from the Q sampling features 59% of H-elimination, against 70% when starting from the T 300 sampling (Figure 6-top). For 100 trajectories, the uncertainty in these values is ±8% for 90% confidence interval. Thus, in spite of some overlap in the extreme of the error bars, it seems that the T 300 sampling really leads to a preference for H elimination in comparison to the Q sampling.

Sampling

There is an extensive patch of ring-distortion crossing seam that can give rise to internal conversion to the ground state. 67 It varies in terms of degree and type of ring puckering. The distribution of ring-distortions at the S 1 /S 0 crossing point given by the Cremer-Pople parameters 68  and Q are shown in Figure 6-bottom. For both samplings, internal conversion to the ground state tends to occur with a substantial puckering degree, between 0.3 and 0.5 Å. Also for both, there is a large variation of types of puckering specially from envelope conformations involving the N atom (E 1 ) to twisted conformations involving C2 and C3 ( 2 T 3 ). It is not possible, however, to distinguish any special trends in either sampling.

Figure 6. Top: Number of S 1 →S 0 hops as a function of the NH distance. For 90% confidence interval, the errors in the fractions of ring distortion and H-elimination are ±8%. Bottom: Cremer-Pople parameters at the S 1 →S 0 hop time for the ring-distortion cases. The corresponding puckering conformations are indicated at the right axis. E n : envelope puckering with atom n below the ring plane (in n E, n is above the plane); m T n : twisted puckering with m above and n below the plane. C2↔C5 and C3↔C4 equivalence was taken into account, such as that, for instance, E 5 is counted as E 2 .

The bias of T 300 towards H elimination is somewhat counterintuitive. We could expect that with a narrower NH-distance distribution than in the Q sampling, the activation of the H elimination would also reduce. The reason for having the opposite output is related to the nonadiabatic events bringing pyrrole from S 2 to S 1 . If we apply Eq. ( 10) to fit the S 1 occupation instead of S 1 +S 2 , we see that the S 2 →S 1 transfer time is  1 = 55 fs in the Q and 89 fs in the T 300 sampling. This reflects the smaller momentum of T 300 trajectories when they reach the S 2 /S 1 crossing seam, compared to the momentum of Q trajectories in the same region. The Landau-Zener model 69 of nonadiabatic transitions helps to understand what happens then: in a two-state crossing, everything else being the same, the larger the momentum, the larger the probability of remaining in the same diabatic pathway.

In this way, as the S 2 /S 1 crossing seam is reached along ring-distortion modes, 37 trajectories with larger momentum are more likely to remain moving along the ring distortion pathway than those with Ring distortion Q: 41% T: 30%

Q T  (°) Q ( Å ) 4 T 3 E 3 2 T 3 2 E 2 T 1 E 1
smaller momentum, which have larger probabilities of changing to another diabatic pathway, resulting in H-elimination.

With the exception of the time constant for transfer the populations from high excited states into S 1 and S 2 ( 12 ), dynamics is clearly faster with Q than with T 300 sampling (see Table 3). The S 1 →S 0 time constant for trajectories following the ring-distortion path is  R = 98 fs with Q and somewhat slower with T 300 , 124 fs. For H elimination, this time increases from  H = 87 to 154 fs between the Q and the T 300 samplings.

Comparison to experimental data

We can summarize the dynamics of pyrrole excited into the * band as schematically illustrated in Figure 7. This figure shows the singlet (neutral) and doublet (cation) states of pyrrole along the Helimination and ring distortion computed by linear interpolation of internal coordinates between the ground state minimum and the optimized S 1 /S 0 intersections of each path. These calculations were done at the same level as the spectrum and dynamics simulations, ADC(2)/(aug)-cc-pVDZ. In the following description, we will only mention the Q-sampling time constants and populations fractions.

Qualitatively, the T 300 description is exactly the same, but with the values provided in Table 2, Table 3, and Figure 6.

Excitation into the 6.1 ± 0.1 eV (Figure 7a) populates the S 3 -S 6 manifold, especially S 4 and S 5 .

After that, pyrrole quickly relaxes to the S 2 state in only 27 fs. S 1 is populated through a S 2 /S 1 crossing induced by out-of-plane ring modes (Figure 7b) within 55 fs. The split between the two deactivation paths happens at this crossing: part of the population continues moving along the ring distortion, and the remaining population planarizes and starts the H elimination. The population following the ring distortion (41±8%) reaches the S 1 /S 0 intersection (Figure 7c) within 98 fs. The population following H-elimination (59%) reaches the S 1 /S 0 intersection (Figure 7d) within 87 fs. Both paths considered, pyrrole returns to the ground state within 90 fs.

The split of population between the two pathways has been systematically neglected in diverse experimental works, which usually take into account only the H elimination. However, photodissociation experiments at 193 nm (6.42 eV) 15 showed that pyrrole splits half-half between C 4 NH 4 + H and C 3 H 4 N + HCN channels. Moreover, the H fragments also split between fast and slow fractions, corresponding to two different mechanisms. Fast H fragments arise from the H-elimination pathway, while slow H fragments arise from statistical dissociation in the hot ground state competing with HCN production. 24 Therefore, the measured systematic increase of slow H fragments between 254 and 193.3 nm 24 is an indication that another relaxation path but H-elimination is becoming more and more important at higher excitation energies. The only known alternative for this relaxation pathway is the ring distortion. HCN and slow H fragments should result from deactivation through ring-distorted pathways.

For this reason, we expect that the ring-distortion/H-elimination branching in the dynamics simulations of pyrrole excited into the * band should also be large, about 50% to cope with the experimental results. The results from Q sampling (41% ring distortion against 59% H elimination; see Figure 6-top) are doing better in this analysis than those from T 300 (30% ring distortion). In any case, either with Q or T 300 samplings, the substantial split between the two reaction paths is a very positive point for ADC (2). Previous full dimensional dynamics of pyrrole carried out with MRCI and TDDFT failed to show such split. 36,38 The comparison between experiment and theory concerning the time to reach the ground state poses more challenges. The experiments usually are based on transient variation of the ionization potential (IP) or H-fragments formation. In both, what is measured is not exactly the internal conversion to the ground state. In the first case, the IP strongly varies along the reaction coordinates, a feature that, although recognized, 19 has not been taken into account yet. As we can see in Figure 7, while in the Franck-Condon region (displacement zero), ionization from S 1 into D 0 requires 4.29 eV (289 nm), in the ring-distorted S 1 /S 0 intersection it requires 9.05 eV (137 nm) and in the H-elimination S 1 /S 0 intersection it requires 9.81 eV (126 nm). In Ref. 16 , up to 3×800 nm (4.65 eV) probe pulses were used to ionize pyrrole, while in Ref. 1 , the probe pulse was fixed at 266 nm (4.66 eV). In both cases, after the excitation, pyrrole should leave the ionization window much before reaching the intersection to the ground state, rendering time constants that should be shorter than the S 1 /S 0 internal conversion time. In fact, these experiments may be only probing the time for pyrrole to leave the Franck-Condon region.

In the case of H-fragment detection, the first problem is that such experimental setup neglects the fraction of the population deactivating through ring-distortion. Second, because these experiments are based on 1s→2s excitation of H atoms followed by its resonant ionization, 17,21 they may detect H + fragments arising from semi-dissociated pyrrole still in the excited state, leading to an underestimation of the lifetime. This happens because the atomic 2s state correlates with the * state of pyrrole, which becomes the S 1 after the S 1 -barrier along the H-elimination path (the barrier is indicated by a * in Figure 7; before the barrier, S 1 is 3s). Therefore, when the 1s→2s excitation takes place, pyrrole may still be in the S 1 state. In Ref. 17 , the time constant associated to H + detection after 200 nm (6.20 eV) excitation of pyrrole is found to be  f = 52±12 fs. How may we compare the simulations to this result? Sapunar and co-workers, 35 based on surface hopping dynamics also carried out with ADC(2), directly compared the S 1 /S 0 deactivation time (in their case, 150 fs) to  f from Ref. 17 . They attributed the difference between these two values to deficiencies of ADC (2) to describe the potential energy surface of pyrrole. However, since H fragments may arise from excited pyrrole, this comparison is not the most adequate.

What should be compared to the experiment is the time that trajectories following the Helimination path cross from 3s to * state. This crossing (indicated by * in Figure 7) happens at the NH distance of 1.29 Å. The computational quantity that holds the best equivalence to the quantity measured is the cumulative probability of the trajectories to have reached this threshold distance.

Fitting this probability with an exponential model similar to that used to fit the transient spectra in the experiments (see supporting information of Ref. 17 ) results in a H-fragment formation time of  * = 64 fs for Q sampling and 100 fs for T 300 sampling (Table 3). Compared to the  f = 52±12 fs experimental value, 17 these results mean, first, that there is nothing in principle wrong with ADC(2) simulations and, second, that Q sampling performs better than T 300 .

Choosing between quantum and thermal samplings

Along this paper, we have shown that the results for excited-state dynamics and the absorption spectrum simulations for pyrrole are very dependent on the method chosen for sampling the groundstate geometry distribution. Unsurprisingly, thermal samplings delivered narrower spectrum and slower dynamics than quantum sampling. Although we have studied a single molecule in this work, most probably these are general trends that should be observed for most of molecular systems investigated with similar methodologies.

We should admit that even before any simulation, we already favored the quantum sampling, as it should be the best representation of a quantum density projected on the classical phase space.

The comparison between the simulations and the experiments reinforced this view, as the Q sampling performed clearly better than the T sampling. Thus, our general prescription is that the Q sampling should be the preferred method to sample initial conditions for dynamics and spectrum.

We are aware, however, that for large molecular systems, the Q sampling may be computationally too expensive. As applied here (see Eq.( 2))-which could be considered the simplest approximation for a Q sampling-normal displacements are required. To get them may simply be prohibitive, for instance, for typical chromophore-protein systems. One possibility is to resort to a hybrid method like that discussed in Ref. 44 , where the Q sampling is applied only to the chromophore, and the remaining system is thermalized around each frozen chromophore geometry given by the Q ensemble. With such a method, only normal modes for the chromophore are needed.

A special situation may occur for systems with a large number of conformational or isomerization minima. In such a case, the ensemble and the derived quantities (spectrum and initial conditions) should be weighted by the corresponding Boltzmann factor of each minimum. Again, this may make the application of the Q sampling prohibitive because one needs to compute normal modes for each minimum. We have faced such a type of system before, studying the spectrum of urocanic acid. 70 This molecule has diverse tautomers and the final spectral shape was determined by few of them. To deal with a system with multiple minima, we may once more resort to a hybrid methodology, where a thermal sampling is initially used to locate the minima contributing with the largest Boltzmann factors, and then applying the Q to each of them.

Here, we have discussed a sub-picosecond process. However, for much longer processes (tens of picoseconds), the Q sampling may also pose other challenges. First, because of the faster dynamics than with the T sampling, the quality of the numerical integration of the Newton equations will degrade faster with the Q sampling, which may require shorter time steps to keep the same level of energy conservation than with T sampling. 40 Moreover, the initial distribution of energy (zero-point energy in each mode) will tend to leak from fast to slow modes, yielding a distribution similar to that discussed for T 912 sampling. Such situation may additionally require the usage of zero-point energy constraints. 71

Conclusions

Absorption spectrum and excited-state nonadiabatic dynamics simulations for pyrrole in the gas phase, initiated from thermal (T) and quantum (Q) samplings, show that:  The ground-state properties of the ensemble (energies and internal coordinates) are strongly dependent on the sampling method. By construction, Q sampling renders mean energies near the zero-point level, while T sampling renders much smaller mean energies, near a thermal distribution. The width of the energy distribution is also strikingly different, much broader with Q than with T sampling.

 The difference between the distribution of internal coordinates in the T and Q sampling depends on the normal frequencies contributing to the coordinate. The larger the frequencies are, the larger the differences.

 Overheating the T sampling until it reaches the zero-point level produces good ground-state energy distributions, but creates artifacts in the internal coordinate distributions.

 For all tested properties, T samplings prepared with either Born-Oppenheimer or molecular mechanics produce similar results.

 The absorption spectrum is strongly dependent on the sampling method. T sampling at 300 K produces too narrow spectra compared to the experiments, while the Q-sampled spectrum delivers satisfactory results. With an overheated T sampling, the spectrum is better, although still inferior to that obtained with Q sampling.

 Qualitatively, T and Q samplings lead to similar results in the dynamics. After excitation into S 3 -S 6 manifold (*), the S 2 and S 1 states are populated under 30 fs, and the ground state is reached between 90 and a 148 fs through two types of S 1 /S 0 intersections caused either by H elimination or by ring distortion.

 Quantitatively, dynamics stating from Q sampling is faster than starting from T by few tens of femtoseconds. When compared to the experimental results, Q sampling performs better in terms of time constants and distribution of reaction pathways. For instance, the time for H fragments formation is 64 fs with Q and 100 fs with T, while the experimental value 17 is 52±12 fs.

Based on these points, quantum samplings based on Wigner distributions are clearly superior to thermal samplings, and they should be preferred for creating initial conditions for spectrum and excited-state dynamics simulations. Cases where Q sampling is computationally too demanding may be treated with hybrid methods.

Although we do not discuss this aspect here, we would like to remark that the exceptionally narrow ensemble distribution arising from MM and BO trajectories may have impact on ground-state investigations as well. This may be especially relevant when the ensemble is used to evaluate the probabilities of reaching transition states. Clearly, geometries and velocities in thermal ensembles populate a too restricted region of the phase space, which may lead to an overestimation of the time to reach transition states nearby. The present simulations have been performed using ADC(2) method with the nuclear ensemble approach for spectrum simulations and surface hopping for dynamics. Compared to the experiments, ADC(2) provided excellent results for vertical excitations and very satisfactory results for spectrum simulations. It has also been able to describe the reaction path split of population during the dynamics between H elimination and ring distortion, a feature expected from photodissociation experiments, which previous simulations based on TDDFT and MRCI missed.

Finally, the present dynamics simulations for pyrrole reinforces the importance of further experiments on this molecule to take into account the substantial amount of internal conversion to the ground state occurring at ring-distorted pathways rather than at H-elimination pathways. The strong variations on the ionization potential along the main reaction coordinates in the excited states should also be considered.

  the T sampling each degree has vT kT   . Because at T = 300 K, vQ  (roughly, ~0.1 vQ  eV) is much larger than vT  (~0.03 eV, see Figure 1-top), this means that simulations based on one or other approach may render different results.

Figure 1 .

 1 Figure 1.Top: Illustration of the classical (T = 300 K) and quantum energies and amplitudes in a onedimensional harmonic potential. Bottom: Ratio between the classical and quantum amplitudes as a function of the wavenumber (Eq. (1)). The proportions in the upper figure correspond to 1600 cm -1 .

  Scheme 1. 1H-Pyrrole.

300T ) are shown in Figure 2 .

 2 The difference between Q and T samplings is striking. The mean kinetic and total energies in the ground state are much smaller in the T 300 and MM 300 T samplings than in the Q sampling. This happens because in the thermal sampling the mean kinetic energy is 3 /2 0.39 eV at N kT  , while it is / 2 1.18 eV ZP E  in the quantum sampling. Apart of a small shift between the distributions, sampling the ensemble either with a BO (T 300 ) or with a MM trajectory ( MM 300 T ) leads to equivalent results.

Figure 2 .

 2 Figure 2. Distribution of ground-state kinetic (top) and total (bottom) energies in the Q and T (BO: T = 300 and 912 K; MM: T = 300 K) samplings. E ZP = 2.36 eV is the harmonic zero-point energy computed at the same QM level.

Figure 3 .

 3 Figure 3. Distributions of NH distances, CN distances, and CNCC dihedral angles for the quantum and thermal (BO: T = 300 and 912 K; MM: T = 300 K) samplings. The ground-state harmonic oscillator wave function for the N-H stretching mode is indicated the upper graph as well (QHO). The triangles indicate equilibrium values. Wavenumbers of the main modes contributing to the variation of each coordinate are also shown. The areas of all graphs are normalized to one.

Figure 4 .

 4 Figure 4.Photoabsorption cross section computed with the quantum and thermal (T = 300 and 912 K) samplings compared to the experimental results.25

Figure 7 .

 7 Figure 7. Potential energy profiles of pyrrole along the H-elimination (positive coordinates) and ringdistortion (negative) paths. The ground and singlet excited states (up to S 7 ) are shown along with the cation (doublet) states D 0 and D 1 . The geometry at left corresponds to the ring-distorted S 1 /S 0 intersection; at the center (displacement zero) is the minimum of the ground state; at right is the H-elimination S 1 /S 0 intersection. The * indicates the top of the S 1 barrier.

  Origin of the state. b Estimated "experimental vertical excitation" (see text); its origin is at 5.698 eV and the band maximum is at 5.95 eV.

	1 1 A 2 1 1 B 1 2 1 B 1 2 1 A 2 1 1 B 2 2 1 A 1 3 1 B 1 3 1 A 2 4 1 B 1 4 1 A 2 2 1 B 2	5.27 [5.22] 5.94 [5.91] 6.04 [6.01] 6.11 [5.96] 6.25 [6.25] 6.33 [6.30] 6.58 [6.55] 6.60 [6.53] 6.82 [6.71] 6.85 [6.79] 7.02 [6.82]	0.000 0.006 0.034 0.000 0.143 0.001 0.012 0.000 0.024 0.000 0.002	5.10 5.85 5.99 5.86 5.98 6.37 6.47 6.43 6.72 6.50 6.63	→3s →3p yz →3s →3s →3p x →*/3p x  →3d yz →3d yz →3p yz →3d yz →3p x /*	5.861 a 6.12 b (0.12)

a

Table 2 the 6

 26 .1 ± 0.1 eV energy window. The number of initial conditions (IC) was taken approximately proportional to the nuclear ensemble (NE) up to a total of 100 initial conditions.

		Sampling	S 1	S 2	S 3	S 4	S 5	S 6	S 7	S 7 -S 25
	NE	Q	0	18	120	218	262	112	2	0
	NE NE IC	T 300 T 912 Q	0 0 0	13 19 0	246 108 19	286 192 30	408 208 36	126 118 15	0 1 0	0 0 0
	IC	T 300	0	0	23	27	38	12	0	0
	Figure 5. for the Q and T 300 samplings. Fraction of trajectories as a function of time in S 0 , S 1 +S 2 , and in the S 3 +…+S 7 groups of states

Acknowledgements

This work has been carried out thanks to the support of the A*MIDEX grant (n° ANR-11-IDEX-0001-02) funded by the French Government « Investissements d'Avenir » program. The authors thank Dr. Oliver Weingart, Dr. Michael Stenrup, Dr. Volker Haigis, and Prof. Rodolphe Vuilleumier for discussions.