article: The 5-year survival of non-small cell lung cancer patients can be as low as 1% in advanced stages. For patients with resectable disease, the successful choice of preoperative chemotherapy is critical to eliminate micrometastasis and improve operability. In silico experimentations can suggest the optimal treatment protocol for each patient based on their own multiscale data. A determinant for reliable predictions is the a priori estimation of the drugs’ cytotoxic efficacy on cancer cells for a given treatment. In the present work a mechanistic model of cancer response to treatment is applied for the estimation of a plausible value range of the cell killing efficacy of various cisplatin-based doublet regimens. Among others, the model incorporates the cancer related mechanism of uncontrolled proliferation, population heterogeneity, hypoxia and treatment resistance. The methodology is based on the provision of tumor volumetric data at two time points, before and after or during treatment. It takes into account the effect of tumor microenvironment and cell repopulation on treatment outcome. A thorough sensitivity analysis based on one-factor-at-a-time and latin hypercube sampling/partial rank correlation coefficient approaches has established the volume growth rate and the growth fraction at diagnosis as key features for more accurate estimates. The methodology is applied on the retrospective data of thirteen patients with non-small cell lung cancer who received cisplatin in combination with gemcitabine, vinorelbine or docetaxel in the neoadjuvant context. The selection of model input values has been guided by a comprehensive literature survey on cancer-specific proliferation kinetics. The latin hypercube sampling has been recruited to compensate for patient-specific uncertainties. Concluding, the present work provides a quantitative framework for the estimation of the in-vivo cell-killing ability of various chemotherapies. Correlation studies of such estimates with the molecular profile of patients could serve as a basis for reliable personalized predictions.
  Worldwide, lung cancer accounts for most cancer related deaths among both men and women [1]. Non-small cell lung cancer (NSCLC) represents the most common type [1]. The success of current treatment choices depends on the extent of the disease at diagnosis; however, overall prognosis remains poor. For locally advanced and metastatic NSCLC, accounting for more than a half of NSCLC incidence [2], the 5-year survival rate ranges between 14% and 1% [1]. The use of cisplatin in combination with another agent remains the standard of care in NSCLC [3]. For patients with resectable tumors, neoadjuvant chemotherapy can be proven particular beneficial in terms of operability, event-free survival, time to distant recurrence and overall survival [4]. However, if treatment fails, a considerable time will have passed during which the tumor may have advanced or even become inoperable [4]. Treatment choices have routinely been based on stage, tumor size, location, lymph node or distant metastasis and overall health status. The exploitation of the molecular profile of cancer cells as treatment selection criteria in NSCLC has only been limited to the consideration of EGFR or ALK mutations as therapeutic targets [5–6]. However, the consideration of the molecular landscape, not to mention the complete genome sequencing, of cancer cells to guide treatment choice is a promising new research area in the field of personalized medicine [7–8]. Mechanistic models that summarize our knowledge on cancer progression are potential candidates to bridge the gap between the gene and molecular world and the prediction of treatment failure or success. Simplified models summarizing cancer biology in few parameters (e. g. growth rate) cannot grasp the great inter-patient and intra-tumor heterogeneity of the disease. Personalized cancer treatment demands for more sophisticated approaches to better trace and understand the biology of cancer and overcome the barriers of current therapeutic practices. The success of such a model to predict the therapeutic outcome depends greatly on the accuracy of the patient-specific estimates of its parameters at the time of diagnosis. In the present work a methodology for the estimation of in vivo cell killing ability of chemotherapy based on data routinely measured in clinical practice has been developed. A discrete event-discrete state, clinically oriented simulation model of tumor response to chemotherapy is applied to translate tumor shrinkage after treatment to cytotoxic efficacy by taking into account the effect of tumor kinetics on treatment outcome. Measure of the antitumor efficacy is the ‘cell kill rate’ that expresses the fraction of treatment-induced cell losses, i. e. the fraction of cancer cells that are lethally hit, after each drug administration. The model addresses primary tumors during their clinical course of life, well beyond their initiation phase. It has been designed to incorporate patient-specific data such as imaging-based, histological, molecular and treatment data. In our modelling approach the high complexity of cancer is reduced to these essential mechanisms that prevail at cellular and super-cellular length scales. Hypoxia due to an insufficient tumor neovasculature, reversible dormancy, active proliferation and spontaneous, starvation-induced or treatment-induced cell loss are among the biological processes considered. Intra-tumor heterogeneity has been implemented via the cancer stem cell (CSC) hypothesis. Different resistance profiles of cancer cells have been taken into account. An essential feature of the model is its capability to simulate different proliferation patterns based on the assumed conditions of tumor microenvironment. Eventually tumor progression is regulated by the interplay of the above considered biological mechanisms. The core simulation model of tumor cell multiplication stems from previous work of the In Silico Oncology-In Silico Medicine Group, ICCS, National Technical University of Athens [9–12]. In the present work, several extensions and adaptations have been performed to account for the modelling of the combination chemotherapy and the mechanism of action of the drugs considered. Furthermore, by recruiting statistical sampling methods to overcome intrinsic uncertainties in model parameter values, a methodological framework has been proposed that enables the application of mechanistic simulation models on the retrospective analysis of accumulated clinical data in a systemic way. In order to facilitate the interpretation of histological subtypes to model parameters, a comprehensive literature survey on non-small cell lung cancer-specific proliferation kinetics has guided the consideration of biologically-reasonable and cancer-specific value ranges for critical proliferation features. Moreover, in terms of better understanding cancer, the effect of key biological processes that drive tumor progression on treatment response are studied and clinical implications are discussed. More specifically, a thorough sensitivity analysis has been performed to decipher which tumor features are determinant of the quantitative estimation of the drugs’ cell killing ability. The tumor features under examination include tumor growth rate, cell cycle time, growth fraction, cell loss and chemo-resistance. Means to improve the accuracy of our estimations by exploiting measurable tumor proliferation kinetics are discussed. Finally, a crucial advancement in relation to previous works is the use of real anonymized clinical data. The methodology is applied and early validated in the case of neoadjuvant chemotherapy treatment of NSCLC with various combinations of the antineoplastic agent cisplatin with gemcitabine, vinorelbine, and docetaxel. A real clinical dataset of 13 patients with primary NSCLC has been recruited for the purposes of this study. The anonymized data originate from the Institute of Pathology and the Clinical Cancer Registry, University Hospital of Saarland, Germany.
  Thirteen patients with newly diagnosed NSCLC are included in the study (Table 1). The patients were treated at the Institute of Pathology of the University Hospital of Saarland between 2006 and 2010. The patients were between 45 and 67 years of age. Nine of the patients were diagnosed with squamous cell carcinoma (SCC) and the rest four with adenocarcinoma (ADC). None of the patients had previous cancer occurrence. All patients had stage III or IV disease according to TNM Classification of Malignant Tumours, 7th ed. There were two patients with distant metastases at the time of diagnosis. All patients received primary systemic chemotherapy prior to surgery. All data exploited by the present study have been provided, following anonymization, through the security framework implemented within the Contra Cancrum European Commision-funded program [13] (Contra Cancrum: Clinically Oriented Translational Cancer Multilevel Modeling’, FP7-ICT-2007-2-223979, http: //contracancrum. eu/). Pseudonomized data from the patients and tumor tissue of the Contra Cancrum project were used according to the approval of the ethical committee Ethik-Kommission, Ärztekammer des Saarlandes, Faktoreistr. 4,66111 Saarbrücken, Germany—approval number: 104/10. (This Committee is responsible for any medical ethics decision in the state of Saarland including all affairs of the medical faculty of the University of Saarland.) As the study has to focus on retrospective data e. g. tissue data from stored tissue samples the ethical committee agreed that informed consent cannot be gained in retrospective cases. All patients were treated preoperatively with a cisplatin-based doublet regimen (Table 2): one patient received a combination therapy with cisplatin (75 mg/m2) plus docetaxel (75 mg/m2), three patients received cisplatin (80 mg/m2) plus gemcitabine (1250 mg/m2), and nine patients were given cisplatin (80 mg/m2) plus vinorelbine (30 mg/m2). The cisplatin/gemcitabine and cisplatin/vinorelbine doublet regimens were given as a three-week cycle, administered two or three times. On day 1 of each cycle the patients received a short IV infusion (10-30min) of gemcitabine (or vinorelbine), followed by a 1–2 h IV infusion of cisplatin, whereas on day 8 only gemcitabine (or vinorelbine) was administered. In the case of cisplatin/docetaxel doublet regimen, both chemotherapeutic agents (1h IV infusion of docetaxel followed by 1h IV infusion of cisplatin) were given on day 1 of a 21-day cycle, repeated three times. After systemic chemotherapy, the patients underwent surgical resection of the primary tumor and/or the regional lymph nodes. The CTs were performed by positioning the patients in supine position with both arms above the head. A spiral CT of the chest with intravenous contrast during mid–breath hold was obtained (two-channel Elscint Twin Flash CT or four-channel MX 8000 [Philips Medical Systems, Best, The Netherlands]) if there were no contraindications. A pitch value of 1. 2 was used. Imaging data were reconstructed at 3mm slice thickness and 1. 5mm reconstruction interval. All images were evaluated by a board-certified radiologist for the presence or absence of pathological lesions. In order to be able to compare model predictions qualitatively or quantitatively with patient 3D CT images, the images taken at different time instances must be properly aligned (registration) and the tumor shape must be determined (segmentation). The patient specific data that have been exploited by the model are the applied chemotherapeutic scheme (drugs, administration instants) and the 3D image of the tumor as segmented from CT imaging data (Table 2). The image processing of the initial DICOM files was implemented as described above. Interactive segmentation was performed at soft tissue window in the present study. The sets of imaging data refer to the primary tumor and were provided for two time instances before, and during or after the completion of the systemic treatment (prior to surgery). The reconstructed images contain information only about the external boundary of the tumor, whereas information concerning any distinct internal metabolic region is absent. Accordingly, the virtual tumor implemented is assumed homogeneous with a shape compliant to the segmented tumor shape. A detailed description of the model’s formulation can be found in [11–12]. The model is built on the concept of a discrete time and space stochastic cellular automaton. In particular, the anatomical region of interest is represented by a regular grid of voxels called “geometrical cells” (GCs). Each GC that belongs to the tumor region is occupied by a cluster of tumor cells that are distributed in various states, representing cell types and phases. The rules that govern the transition of tumor cells between these states, in each GC, are described by the cell kinetic model as presented below. A second set of rules regulate the movement of tumor cells between GCs, aiming at a spatial evolution conformal to the initial shape of the tumor [12]. Tumor propagation is modeled based on the ‘cancer stem cell’ hypothesis and is regulated by the balance between active cell cycle, quiescence, differentiation and death. The tumor population comprises the following five cancer cell categories: a. cancer stem cells, which possess an unlimited mitotic capacity, b. cells of limited mitotic potential (LIMP), c. terminally differentiated cells (DIFF) that have lost the ability to divide, d. cells that have died through apoptosis and e. cells that have died through necrosis. In addition, stem and LIMP cells can be either cycling, distributed in the four phases of the cell cycle (G1, S, G2, M), or quiescent (G0). The rules that describe the transition between the various categories/phases of the cancer cells are depicted in Fig 1. On the top of the developmental hierarchy lie the stem cells that have the ability of self-renewal and differentiation. Two types of stem cell division are allowed: symmetric that gives rise to two stem cells, with probability Psym, and asymmetric that gives rise to one stem and one LIMP cell, with probability 1-Psym. LIMP cells follow a type of aberrant differentiation pathway. After NLIMP divisions, LIMP cells are assumed to generate the population of the DIFF cells. Cycling cells (stem or LIMP) require a time equal to TC in order to progress throughout the active cell cycle. Cellular dormancy is considered to be due to both nutrient deprivation (hypoxia) and lack of growth-promoting stimuli. Proliferating cells found under either one of the aforementioned conditions are assumed to withdraw from the active cycle into a common G0 state upon completion of mitosis, with a mean probability of Psleep. Under conditions of insufficient nutrient supply and oxygenation, quiescent cells can survive for a limited period. Subsequently, they die through necrosis unless the local metabolic conditions stimulate the entry into the cell cycle. A mean residence time of cells in G0 phase equal to TG0 is considered. After its expiration cells reenter cell cycle with a mean probability PG0toG1. Cycling and quiescent cells may die through spontaneous apoptosis, with rate RA. Differentiated cells may die through apoptosis (rate RADiff) or necrosis (rate RNDiff). Apoptotic and necrotic cells are assumed to be present in the tumor bulk for a time length TA and TN respectively, before their final elimination. In our modelling approach, chemotherapy is assumed to affect only cancer cells with proliferative capacity, either cycling or quiescent, depending on the cell cycle specificity of the drug. The activation of apoptosis is regarded as the major mode of action of chemotherapeutic drugs against cancer cells at clinical relevant doses. Chemotherapy-induced apoptosis is implemented through the parameter ‘cell kill rate’ (CKR) that expresses the fraction of cancer cells that sustain lethal damage by the drug (s) and are destined to die. These cells enter a rudimentary cell cycle before the triggering of the apoptotic pathway. The exact time point within the cell cycle when lethally hit cells enter the apoptotic compartment depends on the mechanism of action of the specific drug. We assume that all lethal cell lesions induced by the drug take place instantaneously at the time of drug administration, and disregard any time delays due to the method of administration or the specific pharmacokinetics and pharmacodynamics of the drug. The consideration of such delays would postpone the permanent removal of lethally hit cells from several hours to a few days (S1 Text). Given that the time interval between the last drug administration and the final tumor volumetric measurement is comparatively long enough (Table 2), the bulk of lethally hit cells will have been removed and will not contribute to the final volume of the simulated tumor, in all case scenarios (S1 Text), despite the aforementioned delays. Likewise, since the dosage interval is comparatively long, no drug accumulation is anticipated between consecutive administrations of the drugs (S1 Text). Therefore, the consideration of such delays is not expected to alter the results of the present study with respect to the estimated CKR. Lethal lesions are assumed to occur in all cell cycle phases, since the drug absorption or the damage may take place in other phases too, besides the specific phase during which the final toxicity is manifested. Indeed, for all of the drugs considered in the present study, their ability to kill cells in all proliferative phases has been reported in literature (Table 3). As a first simplified approximation, for cell cycle specific drugs (gemcitabine, vinorelbine and docetaxel), the same value of cell kill rate is assumed for all proliferative phases, corresponding to a mean fraction of cells that will be lethally hit by the drug over the entire cell cycle. For cell cycle non-specific drugs (cisplatin) the same value of cell kill rate is assumed for cell cycle and G0 phases. When drugs are used in combination, their effects are considered additive, i. e. the total number of cells that are lethally hit in each proliferative phase is determined by the sum of the cell kill rates of the drugs. The additive action is applicable only for the proliferative phases (G1, S, G2, M), since cisplatin is the only drug that acts out of the active cell cycle, i. e. during G0 phase. Table 3 summarizes the action mechanisms and the modeling approach for the drugs considered in the present study. In order to study the effect of stem cell-resistance to chemotherapy, the model allows the independent handling of stem and LIMP cell categories. More specifically, the resistance to chemotherapy can be adjusted separately for stem or LIMP cell categories, through the cell kill factor (CKF) parameter. In the current work we have assumed that the CKF of LIMP cells is always equal to 1, i. e. the number of LIMP cells that are lethally hit by chemo is defined by CKR solely. Furthermore, the cell cycle and the cell cycle phase durations can be defined separately for the stem and the LIMP cells. Table 4 presents the input parameters of the simulation model and their range of values according to the conducted literature survey or based on logical assumptions supported by basic science or clinical experience. The effect of the proliferation kinetics of the tumor cells on the estimation of CKR has been explored. Measures of the kinetics of cell proliferation are the relevant parameter space of the model (Table 4) and the derived proliferation features of the simulated tumor, i. e. volume growth rate and tumor cell population constitution. The methodologies applied are based on the one-factor-at-a-time (OFAT) approach [20] and the Latin Hypercube Sampling/Partial Rank Correlation Coefficient (LHS/PRCC) sensitivity analysis [20,21]. The aim is to define a set of model parameter values that corresponds to a cancer representative proliferation pattern in terms of volume doubling time, Td, growth fraction (GF), apoptotic cell fraction (AF), expressed as the apoptotic cell population out of total cells, necrotic cell fraction (NF), expressed as the necrotic cell population out of total cells, and stem cell fraction (SF), expressed as the stem cell population out of living cells. We consider the fraction of newborn cells that enter the G0 phase, Psleep, the necrosis rate of differentiated cells, RNDiff, the duration of apoptosis, TA, the duration of necrosis, TN, and the number of mitotic divisions performed by LIMP cells before becoming terminally differentiated, NLIMP, as the dependent parameters of this multi-constrained problem. The independent variables comprise the rest of the model parameters that regulate tumor proliferation pattern, i. e. cell cycle time, TC, the duration of G0 phase, TG0, the apoptosis rate of stem and LIMP cells, RA, the apoptosis rate of differentiated cells, RADiff, and the symmetric division fraction, Psym, and the above mentioned tumor proliferation features: GF, Td, AF, NF and SF. Initially, biologically representative values, determined based on the cancer-specific cell kinetics and logical assumptions, are assigned to the independent input parameters TC, TG0, RA, RADiff and Psym. Following, the Psleep value is computed so as to achieve the given Td: Psleep=1−e (a+RA) Tc/ (1+Psym) 1−PG0toG1TG0/ (a+RA+1TG0) (5) (derived from Eq (7) in [12]) RNDiff, is calculated in order to achieve the given GF: RNDiff=1−Psym1A (1GF−1) −B−a−RADiff (6) (derived from Eq (47) in S2 Text) where A=a+RAe (a+RA) TC−1 (7) B=1+Psyma+RA+1TG0Psleep (8) By taking into account that the ratio of apoptotic, NA, to proliferating, Np, population is given by: NANP=AFGF1Living/Total=AFGF11−AF−NF (9) ΤA, can been calculated based on Eq (13) in S2 Text: TA= (RANANP+RANANP (1+Psym) PsleepAa+RA+1TG0+RADiffNANP (1−Psym) Aa+RADiff+RNDiff−a) −1 (10) Finally, ΤΝ has been determined for the given NF based on Eq (17) in S2 Text: NNNP=NFGF11−AF−NF (11) TN= ( (1−PG0toG1) TG0NNNP (1+Psym) PsleepAa+RA+1TG0+RNDiffNNNP (1−Psym) Aa+RADiff+RNDiff−a) −1 (12) The value of NLIMP can be estimated based on Eq (46) in S2 Text: ∑n=0NLIMP−12n (1+Psym) n+1=1 (1−Psym) NSNL (13) The ratio of stem to LIMP cells in Eq (13) can be derived based on SF, GF and ratio of dormant to proliferating cells given by Eq (5) in S2 Text: NSNL=SF (1+NG0NP) GF−SF (14) NSNL=SF (1+ (1+Psym) Psleepa+RA (a+RA+1TG0) 1 (e (a+RA) Tc−1) ) GF−SF (15) The above methodology is applied for each NSCLC subtype considered, i. e. ADC and SCC. The derived two value sets of model parameters consist the two sets of baseline values used in the OFAT sensitivity analysis. The Latin Hypercube Sampling method is used to generate a plausible collection of model parameter values that corresponds to virtual tumors having a common proliferation pattern in terms of Td, GF, AF and NF. We consider the fraction of newborn cells that enter the G0 phase, Psleep, the necrosis rate of differentiated cells, RNDiff, the duration of apoptosis, TA, and the duration of necrosis, TN, as the dependent parameters of this multi-constrained problem. The ‘lhsdesign’ function is run to produce N combinations of the independent model parameters: TC, TG0, RA, RADiff, PG0toG1, Psym, NLIMP, and CKF of stem cells. LHS output is modified in the case of parameters with value range other than [0,1] and parameters of integer type, as previously described. For each combination of parameter values, parameters Psleep, RNDiff, TA and TN, are derived from Eqs (5), (6), (10) and (12) as described previously, so as to achieve the given Td, GF, AF and NF. The above methodology is applied for each NSCLC subtype considered, i. e. ADC and SCC. The derived two LHS samples of parameter values are used for the estimation of the cell killing ability of the various cisplatin-based doublet regimens. The cell kill rate of cisplatin, CKRcis-DDP, is adapted to the observed, in the case of clinical tumors, or assumed, in the case of sensitivity analysis, tumor size reduction. The adaptation of CKRcis-DDP has been performed automatically using an optimization procedure (the built-in function ‘fzero’ of Matlab), that returns the value of the CKRcis-DDP for which the difference between the observed or given volume reduction and the simulation outcome is zero. During this optimization procedure the rest of the model parameters, including the CKR of the second drug, are kept constant.
  The model supports the division of tumor volume into regions of distinct metabolic activity (i. e. necrotic, quiescent and proliferative). However, due to the lack of relevant info in the segmented tumor images, only a single metabolic region is considered. This implies that all model parameters affecting tumor cell kinetics have the same value throughout a tumor’s volume. Furthermore, the kinetics of tumor cell proliferation (i. e. phase durations and transition rates or fractions) are assumed constant throughout simulation, reflecting the mean values over the simulation time window. Such an approximation is considered applicable for a relatively short time interval compared to the tumor’s lifetime. Finally, the model aims at simulating fully developed, clinical tumors, well beyond their initiation phase. For such a biological system it is rational to assume that a state of population equilibrium and of balanced growth has been achieved. The methodology developed for the initialization of tumor cell populations [11–12] ensures conditions of balanced growth at the beginning of the simulation. A tumor characterized by space and time invariant growth kinetics parameters and population equilibrium would demonstrate a grossly exponential growth pattern [12]. Even though it is generally accepted [23] that tumors grow in a Gompertzian manner during their lifetime, i. e. growth progressively slows down as the tumor enlarges, tumors that follow exponential law have been reported frequently in literature. Indeed, several lung cancer cases reported in literature exhibit fairly constant growth rates on the logarithmic scale [24–27] over prolonged periods of time, ranging from months to years. It has been argued that these tumors were probably still in the exponential phase of Gompertzian curve at the time of diagnosis and clinical surveillance [26]. Fig 2A and 2B show the simulated time course of a macroscopically homogeneous tumor under free growth or treatment conditions respectively. In the case of treatment a typical scheme for NSCLC has been assumed that is based on the combination of two chemotherapy drugs i. e. cisplatin and gemcitabine. The doublet regimen is administered three times as a three-weeks cycle. On the first day of treatment cycle the patient is given both gemcitabine and cisplatin. On the same day of the following week (day eight) only gemcitabine is administered. Qualitatively a fairly expected and reasonable tumor dynamics behavior can be easily noticed. An exponential pattern of growth in the absence of treatment and tumor regression followed by repopulation after each chemotherapeutic session in the case of treatment are successfully demonstrated. Due to the lack of clinical data related to the presence of internal regions of distinct metabolic activity (necrotic, quiescent, proliferating), the analysis has focused on macroscopically homogeneous tumors, in the sense that all model parameters affecting tumor cell kinetics are spatially invariant throughout tumor volume. Moreover, the tumor is assumed to have a spherical shape of 925 mm3 initial volume (approximately 12mm initial diameter). A typical scheme has been considered that is based on the combination of two chemotherapy drugs i. e. cisplatin and gemcitabine. The doublet regimen is administered three times as a three-week cycle. On the first day of treatment cycle the patient is given both gemcitabine and cisplatin. On the same day of the following week (day eight) only gemcitabine is administered. The virtual tumor is allowed to grow for one week prior to treatment. A 72% treatment-induced shrinkage (final tumor equal to 258 mm3) is assumed on simulation day 49. The LHS/PRCC method has been chosen to investigate the effect of the model input parameters as well as of selected macroscopical proliferation features of the tumor, i. e. volume growth rate and tumor cell population constitution, on the estimation precision of CKRsum, by taking into account the presence of correlations between the inputs. The histology, treatment and response data of 13 lung cancer clinical cases have been retrospectively analyzed in order to evaluate the in vivo efficacy of three different cisplatin-based doublet regimens: cisplatin plus gemcitabine, ciplatin plus docetaxel and ciplatin plus vinorelbine. The tumor proliferation features used for the OFAT sensitivity analysis have been also considered here. More specifically, for the clinical cases classified as ADC the assumptions are as follows: (a) Td = 225d, (b) GF = 0. 18, (c) AF = 0. 01, (d) NF = 0. 02 and (e) SF<0. 001. For the clinical cases classified as SCC we have considered: (a) Td = 109d, (b) GF = 0. 36, (c) AF = 0. 01, (d) NF = 0. 20 and (e) SF<0. 001. Furthermore, the ranges of the model parameters are the ones considered for the LHS/PRCC analysis, identical for both ADC and SCC subtypes, and are given in Table 4. LHS has been run to generate 2000 combinations of parameter values that fulfill the above requirements, as described in Material and Methods. After excluding combinations that result in biologically non-relevant tumors, i. e. negative cell class transition rates Psleep and RNDiff, or in tumors with non-relevant proliferation dynamics, i. e. Td beyond ±10% of the above assumptions and/or stem cell fractions higher than 1‰ (based on SoA results) and/or TA>24h and/or TN>200h and/or RNDiff>0. 02h-1, a subset of N = 146 and 175 combinations has been finally considered for the SCC and ADC cases respectively. For each clinical case and for each of the N sets, the CKRcis-DDP, and therefore CKRsum, that results in the observed volume reduction is estimated. The cell kill rate of the second drug, CKRB, has been assumed equal to 0. 2 with the exception of clinical case 9 for which a value equal to 0. 1 has been assumed. This selection of CKRB value has been made with the criterion to have a solution for CKRcis-DDP for all the N LHS combinations considered. Fig 6 shows the median (dot), 10th and 90th percentile (lower and upper boundary of solid line) and minimum and maximum (lower and upper boundary of dashed line) of the N CKRsum estimates for each clinical case. Even though the estimated values of CKRsum vary considerably for all clinical cases, the vast majority of the values (80%) is confined in a much smaller range, not exceeding, with the exception of one case, an average of 0. 1 units around the median value. We observe that the treatment efficacy is very low in 2 cases (15. 4% of the total cases) with a median CKRsum below 0. 3, low in 2 cases (15. 4% of the total cases) with a median CKRsum approximately 0. 38, moderate in 6 cases (46. 1% of the total cases) with a median CKRsum between 0. 47 and 0. 57, high in 2 cases (15. 4% of the total cases) with a median CKRsum approximately 0. 61 and very high in 1 case (7. 7% of the total cases) with a median CKRsum approximately 0. 71. The efficacy of the three cisplatin-doublet regimens seems comparable; however, the sample size is too small, especially for cisplatin/docetaxel and cisplatin/gemcitabine regimens, to draw statistically significant conclusions. In Fig 7 the scatterplots of the N CKRsum values for one SCC and one ADC clinical case are shown versus indicative model parameters and the proportion of dormant and terminally differentiated cells. We observe that for both ADC and SCC cases, the highest values of CKRsum are associated with short cell cycle durations (Fig 7A and 7J), shortly lived terminal differentiated cells (or equivalent removal rate of DIFF cells above 0. 01h-1) (Fig 7B and 7K) and resistant stem cells by a factor of 2 or above (or equivalent CKF of stem cells<0. 5) (Fig 7G and 7P). For the SCC case the highest values of CKRsum are found for low Psym (Fig 7C) and low Psleep (Fig 7D). Furthermore, tumors with higher fractions of dormant cells, at the expense of terminally differentiated cells (GF is constant), tend to require a higher CKRsum to achieve the same volume reduction (Fig 7H, 7I, 7Q and 7R). The three mutually perpendicular dimensions of the surgically resected tumor has been provided for the 12 out of the 13 patients (Table 7). An early validation of the proposed model/methodology is attempted here, by comparing the post-surgical excision measurements with the simulated tumor size at the time of surgery. It should be noted that the size of the resected tumor has not been considered for the estimation of CKRcis-DDP in the previous section; however, it is exploited here for validation purposes. In particular, for each clinical case and for each of the N parameter sets with corresponding CKRcis-DDP derived in the previous section, an extended simulation has been performed starting from the first CT acquisition until the date of surgery. The extra time interval simulated for each clinical case is reported in Table 7. For three clinical cases, additional drug administrations have taken place during this interval and have been modelled assuming the estimated values of CKRcis-DDP. Fig 8A depicts the boxplots of the N predicted tumor sizes for each clinical case in comparison with the real ones. Since it is customary in clinical practice to assess the size of a solid tumor by measuring its maximum dimension (s), tumor size is expressed by means of an equivalent diameter. The predicted tumor volume, V, is translated to the equivalent diameter, d, assuming that the tumor is approximated by a sphere of the same volume: d=2×3V4π3 (16) The equivalent diameter of the surgically resected tumor is calculated as before. The volume of the resected specimen is determined by its three mutually perpendicular dimensions, after assuming that its shape is approximated by a triaxial ellipsoid [41] (Table 7). In Fig 8B the predicted volume reduction is compared with the real one. Despite the large deviation observed between the minimum and the maximum values in both equivalent diameter and volume reduction graphs, the majority of the predictions are in fact confined around the median value. To that respect, the first, second (median) and third quartiles are considered for the validation of our predictions. A good agreement between the prediction and the actual data is observed for the six out of twelve cases. In particular, the equivalent diameter of the resected tumor and, hence, the observed volume reduction, falls between the first and third quartiles of the predicted ones in three cases (# 2,9, 13). Moreover, a deviation less than 10% (corresponding to less than 4mm difference) between the predicted median diameter and the one of the resected tumor is observed for three cases (# 1,6 and 7). However, due to the small size of the initial tumor, a somewhat larger deviation is observed in the case of volume reduction for cases # 1 and 6. It is noteworthy that the cases # 2 and 6 received additional therapy during the extra simulation period. A slightly higher deviation (12% to 20%) between the equivalent diameters of the median of the predicted values and the resected specimen is observed for the cases 5,8 and 11. An even lower deviation that does not exceed 14% exists for volume reduction in the same cases. It is worth mentioning that for these cases the surgery takes place shortly after the second CT acquisition (Table 7), corresponding to a limited (after the second CT acquisition) simulated period; therefore, it is possible that any uncertainties in estimating the true tumor volume by both the imaging techniques and the post-surgical excision measurements may be comparable to the change in tumor sizes and play a role in the observed deviations. For example, volume estimation errors of 10 to 20% are typically expected (more prominent for small lesions) [15], whereas based on our simulation results the volumetric change in a period of only a few days is less than 10%. The less agreement is observed for cases # 3,4 and 10. The model seems to underestimate the cell killing efficacy of treatment and predicts a tumor at surgery that deviates from the resected tumor around 1cm for # 4 and #10 and 2cm for #3. However, for #10 the median volume reduction differs only 15% from the real one. The tumor proliferation profile considered may be a possible source of the observed discrepancies. Our literature study (S3 Text) has revealed a considerable intra-patient variability in the macroscopic characteristics of NSCLC. It is therefore expected our assumptions to deviate from reality for some of the cases. For more accurate results, the clinical data should contain information about the proliferation profile per patient as well. By taking into account the existence of uncertainties in estimating the true tumor volume as described above, the limited number of data sets available for the estimation of CKR and the consideration of a single proliferation profile differentiated only on the grounds of the two histological subtypes of the available clinical cases, the gross agreement between the predicted tumor values and the prost-surgery measured ones indicate that the proposed methodology/ model has a clear clinical potential. Additional adaptation and validation work has been planned in order to further adapt and validate both the model and the methodology.
  Cancer is a highly complex phenomenon. Prediction of treatment failure and selection of the optimal therapeutic strategy can be greatly potentiated by the use of mechanistic models that summarize our knowledge on cancer progression and response to treatment. However, up to now, the adaptation, validation and establishment of mechanistic models in predictive medicine has been hindered by the limited availability of data directly linked to the input parameters. Needless to say, in several occasions, no established method exists for the direct measurement of model parameters. Determinant for the personalized prediction of cancer response is the resistance/sensitivity profile of tumor cells to the specific anti-tumor regimen [11,42,43]. In-vitro experiments for the evaluation of drug cytotoxicity on cancer cell lines have been proven only indicative, as they cannot reproduce the conditions of tumor microenvironment [44]. On the other hand, mechanistic models may provide means to estimate in vivo drug efficacy. In the present paper, a range of plausible values of the cell killing efficacy of chemotherapy has been estimated for NSCLC treated with various cisplatin-based doublet regimens by using a mechanistic model of cancer response to treatment. The cell killing efficacy of chemotherapy is expressed as the sum of cell kill rates, CKRsum, of the two drugs considered. Our modelling methodology takes into account cell re-population between drug administrations. Parameterization has been driven by principal mechanisms characterizing cancer biology, namely uncontrolled cell proliferation, reversible dormancy, clonal heterogeneity, hypoxia and attenuated apoptosis. Two types of tumor cell growth have been considered, corresponding to the ADC and SCC histological subtypes, respectively. A comprehensive literature survey has guided the selection of the most appropriate values of characteristic proliferation kinetics for the two histological subtypes. The more aggressive nature of SCC in respect to ADC has been reflected on the higher Td, the higher proliferation, i. e. GF and stem cell frequency, and the greater extend of necrosis. As a first approximation, inter-patient variation in tumor grade has not been taken into account, and a single differentiation profile has been assumed for each histological subtype. Since, in the patient cohort considered, most cancers are moderately to poorly differentiated, terminally differentiated cells have been assumed to comprise the minority of non-dividing cells, with a lifespan of only a few days. As a consequence, the majority of non-dividing cells are found at a reversible dormancy. It is noteworthy that long-lived differentiated cells, with a lifespan of a few months, similar to the normal bronchial epithelium, would most probably be associated with well-differentiated tumors. A number of sensitivity analyses have been performed in order to investigate the effect of model parameters and cancer proliferation kinetics on the estimation of CKR. OFAT analysis has been chosen as an essential first step in order to gain insight into the mechanisms, verify the biological relevance, validate the assumptions and reveal non-linear patterns or extreme behaviors of individual parameters. In particular, OFAT plots and local sensitivity studies for each model parameter and for each histological subtype have been presented. Subsequently, the OFAT analysis has been supplemented by a number of global sensitivity studies performed by applying the partial rank correlation coefficient (PRCC) method on a latin hypercube sample (LHS) either of the input parameters of the model or of characteristic cancer proliferation kinetics. For the global analyses no distinction between ADC or SCC subtypes has been considered. Although all input parameters are conceptually independent, the fact that the tumor as a whole should satisfy certain macroscopic characteristics (such as volume doubling time, typical fraction of stem cells etc.) dictates numerical correlations among the input parameters. The PRCC method, which by definition estimates parameters’ sensitivity, after controlling for the effects of the remaining parameters, is suitable for models with correlated inputs. A thorough investigation of input correlations is out of the scope of the present work. However, we intend to address it in subsequent publications. The local sensitivity analysis has demonstrated that from the regulating proliferation kinetics those with the highest influence on the estimation of drug efficacy are the fraction of stem cells that perform symmetric divisions and the fraction of newborn cells that withdraw from the active cycle to a dormant phase. Overall, the ADC representative case seems to be less sensitive to all model parameters. The global sensitivity analysis have revealed that the parameters related to the proliferation kinetics of stem and LIMP cells are the most important determinants of the accuracy of CKR estimates. However, these parameters seem to be reduced to primarily two variables, the Td and the GF at the time of diagnosis. In other words, what is of importance is not only the value of the parameters itself but also the value combination of the parameters that regulate Td and GF. Therefore uncertainties in the value of the model parameters could be compensated for by an estimation of Td and GF at diagnosis. More specifically, the global sensitivity analysis has demonstrated that the rate that a tumor grows in the absence of treatment, expressed by Td, has the strongest effect on the estimation of drug efficacy. This variable reflects the rate at which the tumor replenishes its treatment-induced losses. The estimation of Td would involve at least two volumetric measurements of the tumor prior to treatment. The use of CT imaging, and especially high-resolution CT, is necessary to derive more accurate estimations that take into consideration the usually asymmetrical growth rate of the tumor in the three dimensions [45]. The interval between CT measurements should be chosen with caution. Yankelevitz et al [46] has suggested that a minimum of 12% increase in the cross-sectional area between successive CT examinations is required for the reliable estimation of Td. For moderate to slow growing tumors an interval of 25 days would be required [45]. Shorter intervals would be needed for fast growing tumors. For example, based on the formula used in [45], a tumor with Td = 26 days would require a time interval of 6 days between successive CTs. Such waiting periods are common in the case of lung cancer. For half of our patient cohort, the time elapsed between diagnosis and treatment onset is on the order of 20 days and more, so the acquisition of two volumetric measurements prior to treatment might have been feasible. The GF at diagnosis seems to be the second factor of importance for the estimation of drug efficacy. GF estimations based on Ki-67 protein are routinely performed in clinical practice. In general, a deviation of Ki-67 expression levels between biopsies and surgical samples for the same patients has been observed, which has been attributed to differences in sampling, recruitment and Ki-67 staining procedures [47]. Additionally, since a tumor may include areas of low or high proliferation activity, a single small biopsy may not be accurate. However, a significant correlation for this marker between biopsies and resected tumors has been demonstrated [47,48]. Another method that has been used in literature for the estimation of GF is the proliferating cell nuclear antigen (PCNA). However, the latter systematically gives higher values than Ki-67 and is considered less reliable [49]. Based on the results of the sensitivity analyses and Fig 7, the duration of cell cycle, TC, seems to have a relatively significant effect on the estimation of treatment efficacy. Many techniques have been reported in literature, either in vivo or in vitro, for the measurement of cell cycle time, TC, and time-dependent indices in general; however, none has been established and routinely used in clinical practice, mainly due to their slow or laborious nature [23,50]. The best proposed in vivo technique is based on in vivo pulse labelling with a thymidine analogue, either bromodeoxyuridine (BrdUrd), or iododeoxyuridine (IdUrd), followed by flow cytometry analysis and application of relative movement (RM) method [23,50,51]. It requires just a single biopsy, taken a few hours, e. g. 6–8 hours [23], after BrdUrd (or IdUrd) administration. The method allows for the estimation of the duration of DNA synthesis, TS, and the labelling index, LI, which corresponds to the proportion of cells in S-phase. As a first approximation, TC can be calculated based on the equation TC = (TS/LI) *GF, on the condition that the GF is known. More elaborate relationships can be used that take into consideration the non-uniform distribution of an exponentially growing cell population through the cell cycle. In [52] a method to estimate TC and GF from flow cytometric analysis of a single tumor sample after BrdUrd labelling is proposed. Apoptosis is the main cell death mechanism triggered by chemotherapy. However, the global and local sensitivity analyses indicate that both the kinetics of apoptosis and the measure of apoptotic index at diagnosis have no effect on the estimation of CKR. The model’s robustness to the kinetics of apoptosis can be explained based on the rapid nature of the phenomenon. Since apoptosis is a process completed within a few hours, it is characterized by a time scale that it is more than one order of magnitude lower than the kinetics of the rest of the cell populations (days to weeks) and the simulation time window, i. e. the time interval between tumor volume measurements (month (s) ) or the time interval between the last drug administration and the last volume measurement (week). The model seems also rather robust to the kinetics of necrosis and the extent of necrosis at diagnosis. It is noted that the results of the local and global sensitivity analyses correspond to relatively fast necrosis kinetics, with TN not exceeding ten days, and a low to medium extent of necrosis. Furthermore, the importance of CSC resistance in short-term treatment response seems to be moderate due to the very low frequency of cancer stem cells. Finally, for the value range considered, both the local and the global sensitivity analyses demonstrate that the kinetics of terminally differentiated cells have a trivial effect on the CKR estimation. However, as indicated by Fig 3M, a significant sensitivity is expected in the case of long-lived tumor terminally differentiated cells with a lifespan of the order of months. Even though statin has been suggested as a marker of differentiated, non-proliferating cells [53], currently there is no established standard to distinguish terminally differentiated from quiescent cells and, hence, quantify them or measure their kinetics. Finally, this paper attempted an estimation of cell killing efficacy of cisplatin-based doublet regimen by adapting the model to 13 real clinical cases of NSCLC. The purpose of this process has been the determination of set (s) of values of the input parameters of the model that correspond to the clinical case under test (type of cancer, patient, treatment). CKR estimation can be viewed as a two-step adaptation process. The first step aims at adapting the growth kinetic parameters of the model (Table 4). As the modelling approach provides the ground to simulate tumors of diverse growth kinetic characteristics, by assigning the proper values to the model parameters, narrowing the vast number of possible value combinations, ideally to one, for a specific clinical case is a critical first step. In this context, the exploitation of measurable proliferation indices, such as Ki-67, or data that could allow the precise estimation of the tumor growth rate, such as at least two imaging data before therapy, is decisive. However, due to non-availability of patient-specific tumor cell kinetics data, a literature review survey has provided biologically reasonable values for critical proliferation features. The second step involves the adaptation of pharmacokinetic and pharmacodynamic parameters of the model in the case of chemotherapy. The driver in this process is the achievement of volumetric match between the clinical data, i. e. the observed shrinkage in volume, and the result of the simulation, after applying the treatment regimen that was administered in the clinical case under test. The pharmacokinetic and pharmacodynamic properties of the agents considered were summarized under the umbrella of one parameter, the apparent CKR. The purpose of the adaptation process was to evaluate a range of CKR values of the above parameter for the clinical case under test (treatment regimen, dose, patient) based on the results of the adjustment of the first step and the recorded change in tumor size due to the treatment. Since the regimen administered consists of two chemotherapeutic agents, it is not possible to accurately determine the CKR of each drug from the data provided, even in the ideal case of the availability of all required proliferation indices and tumor free growth kinetic features that would enable an excellent fitting of the model parameters to the clinical case examined. In the study presented here an ‘apparent’ combination of the CKR of the drugs involved, for the virtual tumor implementations considered, is determined, by assuming an arbitrary value for one of the drugs. In the clinical data set considered, the majority of the patients (6/13) appeared to have a moderate CKRsum, whereas 2/13 patients had a very low CKRsum and one had a very high CKRsum. Based on the results, the cell killing efficacy of treatment varies considerably among the patients, a fact that verifies the importance of a timely estimation of treatment response through in silico simulations. The accuracy of the conclusions on the treatment effect that was attempted in this paper was subject to the uncertainties with respect to the tumor kinetics. As the number of available data was the minimum required for the model adaptation, averages or ranges that were encountered in the literature were used. More specifically, the data were limited to the assessment of the tumor size in two time points (one before and the other after or during treatment) and to information about the treatment, whereas the data regarding the proliferation features of the tumor was absent. It is noteworthy that the provision of additional tumor size measurements, during the course of therapy, could allow for more accurate estimates. Furthermore, any errors in the estimation of the tumor size were not taken into account, as it is estimated that the uncertainty that they can introduce is much smaller than the lack of knowledge of the characteristics of free tumor growth. Despite the aforementioned sources of uncertainties, the short-term predictive potential of our model has been clearly demonstrated (Fig 8). The gross agreement between the tumor size predicted based on our CKR estimates and the real post-surgery measured one (Fig 8), substantiates the potential utility of the proposed cancer model/methodology for personalized predictions of treatment response at diagnosis. Studies between in vivo CKR estimates, as attempted in the present work, and the molecular profile of a patient, hold great potential for the identification of resistance and sensitivity profiles. For that purpose, statistical learning techniques, either classification or regression ones, can be recruited to distinguish among drug sensitivity phenotypes and to quantify tumor cell response to treatment, based on whole genome analysis data, gene expression profiles and transcriptomic or proteomic signatures. In vivo CKR estimates can be utilized together with the molecular data to train these machine learning models. Eventually, such studies could allow for the prediction of treatment outcome based on the molecular or genetic data of a newly diagnosed patient and will constitute a future step of our work. Our preliminary efforts in this direction [54,55], attempting to correlate prednisone CKR estimates with the pathway gene expression data in paediatric acute lymphoblastic leukaemia, have shown promising first results.
summary: Less than 14% of medically treated patients with locally advanced and metastatic non-small cell lung cancer are expected to be alive 5 years after diagnosis. Standard therapeutic strategies include the administration of two drugs in combination, aiming at shrinking the tumor before surgery and improving overall survival. Knowing the sensitivity profile of each patient to different treatment strategies at diagnosis may help choose the most appropriate ones. We develop a methodology for the quantitative estimation of the cytotoxic efficacy of cisplatin-based doublets on cancer cells by applying a simulation model of cancer progression and response. The model incorporates the proliferation cycle, quiescence, differentiation and loss of tumor cells. We evaluate the effect of in vivo microenvironment of real tumors, as expressed by measurable tumor proliferation kinetics, such as how fast the tumor grows, the percentage of cells that are actively dividing, the resistance of stem cells, etc. on treatment outcome so as to derive more accurate estimates. A literature survey guides the selection of values. The methodology is applied to a real clinical dataset of patients. Correlation studies between the derived cytotoxicities and the patients’ molecular profile could lead to predictions of treatment response at the time of diagnosis.
keywords: cell death
  medicine and health sciences
  cell cycle and cell division
  cancer treatment
  cell processes
  carcinomas
  cancers and neoplasms
  cell differentiation
  oncology
  developmental biology
  stem cells
  adenocarcinomas
  animal cells
  lung and intrathoracic tumors
  cell biology
  apoptosis
  biology and life sciences
  squamous cell carcinomas
  cellular types
  non-small cell lung cancer
title: In Silico Oncology: Quantification of the In Vivo Antitumor Efficacy of Cisplatin-Based Doublet Therapy in Non-Small Cell Lung Cancer (NSCLC) through a Multiscale Mechanistic Model