Computer Science & Systems Biology

Cancer chemoresistance (including adaptive resistance) has emerged as a barrier in developing successful chemotherapeutic strategies. We use Monte Carlo simulation based single cell analysis to provide insights into the regulatory mechanisms for generating chemoresistance under TRAIL (death ligand) induction. Based on stochastic computer simulations we elucidate systems biology of cancer cell apoptosis (at the level of single cells) and search for an optimal death ligand from a group of recently studied TRAIL affinity variants. In addition to assessing the population level behavior in cell death activation under induction of TRAIL/TRAIL-variant, Monte Carlo approach allows us to analyze cell-to-cell stochastic fluctuations in time-to-death that has implications for generating resistant cancer cells. We discuss application of Monte Carlo simulations in the context of developing more personalized approaches in treating various cancers. Initial findings indicate single cell in silico approaches can be utilized for disease subtype classification and in characterizing a given tumor, and, for finding an optimal strategy (such as network modules to target and ligands needed) in targeting a given tumor.


Introduction
Resistance to apoptotic cell death is a defining feature of malignant phenotype in cancer cells. Several oncogenes and tumor suppressor genes have been found and in many cases an increase in anti-apoptotic proteins has been implicated [1]. Even though cancer cells are known for their remarkable apoptosis resistance, over the years various chemotherapeutic agents (that includes genotoxic agents, oncogenic kinase inhibitors, direct apoptosis inducers) have been found to selectively eliminate cancer cells, raising the hope for curative therapies for cancer. In spite of such early promise it was becoming increasingly clear that cancer cells are capable of overcoming chemotherapeutic stress often through generating various chemoresistant phenotypes. However, the cellular and molecular mechanisms for selective apoptotic induction in cancer cells by chemotherapeutic agents, and, the mechanisms for generating chemoresistance, were not clearly elucidated. Lack of such mechanistic elucidations hindered progress in cancer biology in spite of prior bioinformatic analysis that revealed abnormalities in cancer genome and epigenome.
Study of systems level regulatory mechanisms of apoptotic cell death at the level of single cells, including that carried out in this lab, has revealed intricate mechanisms of cancer cell apoptosis that arise from a complex interplay between genetic and stochastic factors [2][3][4][5][6][7]. We have used Monte Carlo (stochastic) simulations to elucidate the type 1/type 2 phenotype choice in apoptotic cell death regulation [4,6]. Type 2 signaling regulation can lead to a slow activation phenotype with all-or-none activation and substantial cell-to-cell variability [8][9][10]. Large cell-to-cell stochastic fluctuations in slow type 2 apoptosis may allow survival of a fraction of cells leading to a bimodal behavior in cell death (fractional cell killing of kinetic origin) and eventual generation of apoptosis resistant phenotypes [2,11,12]. Thus, cell-tocell stochastic fluctuations in apoptosis seem to hold the key to many previously unanswered questions in cancer biology and cancer therapy. We have also studied apoptosis activation under combinatorial treatment options that may minimize the extent of fractional cell killing [13]. Recently, we have engaged in developing Monte Carlo models that started elucidating the mechanisms for TRAIL induced bimodal apoptotic activation (of purely stochastic kinetic origin [14]) leading to fractional cell killing and generation of resistant phenotypes. Affinity variant TRAIL ligands that are now available through protein engineering [15,16] are expected to modify the extent of fractional cell killing and chemoresistance generation in cancer cells but needs to be probed using single cell approaches.
Theoretical analyses, including that carried out in this lab, have provided crucial mechanistic insights into the origin of cell-to-cell variability in apoptotic cell death. Results obtained from our Monte Carlo simulation indicate that certain pro-and anti-apoptotic protein ratios (such as DR/DCR, Bax/Bcl2 and Smac/XIAP) keep the apoptotic threshold under tight regulation (in healthy cells) [2][3][4]7]. Below a threshold ratio (say of DR/DCR under TRAIL induction) cell death activation is inhibited (or remain low) whereas above the threshold cell death becomes increasingly robust. Clearly, cell-to-cell variability in cell death activation can originate from fluctuating levels of certain pro-and anti-apoptotic proteins (extrinsic noise through genetic/epigenetic variations and stochastic fluctuations in gene expression and protein synthesis). Unexpectedly, inherent stochastic fluctuations in apoptotic signaling reactions (intrinsic noise) may also generate significant cell-to-cell variability (with a slow activation phenotype), especially when cells operate close to transition thresholds of death activation (determined by ratios such as DR/DCR) [2,3,7,9]. Cancer cells seem to exploit both kinds of mechanisms, extrinsic and intrinsic noise, for generating a slow activation phenotype with large cell-to-cell variability to evade cell death. Hence, cells that are slow to activate the death pathway would get protected and may generate network), a molecule is sampled randomly and diffusion/reaction moves are carried out in proportion to probability constants in the governing master equations [18,19] (ensuring that the simulation captures the correct dynamics [20]). The current simulation model is a hybrid Monte Carlo that combines the following two approaches: 1) a probabilistic rate constant based (implicit free energy) kinetic Monte Carlo simulation for various reaction moves such as diffusion (of intracellular proteins), binding/unbinding and catalytic cleavage and 2) an explicit free energy based kinetic Monte Carlo that captures diffusion and clustering of receptor molecules utilizing free energy function for receptor-receptor interactions. In the implicit kinetic Monte Carlo simulation probabilistic rate constants are obtained from experimentally measured kinetic constants. Diffusion probabilities for different species (depending on intracellular/membrane-bound and free/complex/signalosome) were taken such that their diffusion rates mimic that obtained in biophysical experiments (such as for membrane receptors ~ 0.1 µm 2 /sec); this mapping allows us to set the time scale of MC simulation ∆T=10 -4 s -1 . The probability for dissociation/ catalytic reactions is then estimated by multiplying the corresponding experimental kinetic constants with ∆T (P off =k off ∆T). Association probability constants are obtained using MC simulations for molecular diffusion and binding/unbinding reactions (P on estimated from experimental k on values). For diffusion of membrane receptors (under the action of receptor-receptor attractive interactions), experimental kinetic constants (similar to k on /k off ) are not readily available and an explicit free energy based algorithm seems to be a natural choice. Energy values for receptor-receptor interaction was assumed to depend on ligand binding status in such a manner that free receptors do not form stable clusters (unless specified by previous experiments) but receptor clustering is inducible upon death ligand binding (previous work by us, Ref. [7,21] has details of energy parameter values). Prior Monte Carlo studies of Ising magnets and liquid-liquid phase separations [22] and biophysical measurements of lipid phase separations [23], have guided our estimation of energy parameter values namely free energy for receptor-receptor interactions (it is possible to vary free-energy values for various TRAIL variants). It is expected that ligand binding drives the cellular system away from its resting state ("live" state for a cell) and the subsequent dynamics is captured by Monte Carlo kinetics (with reaction moves satisfying detailed balance conditions). For probabilistic rate constant based implementation of Monte Carlo, detailed balance is implicitly satisfied through the ratio P on /P off . For explicit free energy based implementation of Monte Carlo, Metropolis method was used as the acceptance/rejection criteria (Metropolis sampling is based on detailed balance condition) [22] when receptor diffusion moves were carried out. It is expected that, in the absence of exact details for the membrane parameters (lipid composition, free energy of receptor-receptor interaction for ligand bound receptors), the hybrid Monte Carlo method used here should capture some of the dynamical mechanisms (including type 1/type 2 phenotype choice) for death and decoy receptor mediated activation of the membrane proximal death module (as DR/DCR ratio is varied). The time-scale of initiator and effector caspase activation that emerges from current simulations is consistent with that obtained in prior experimental studies for death ligand induced apoptotic activation. One may also use a max energy difference based criteria [24] that could capture the combined interaction effect of many receptors when energylowering receptor diffusion moves are carried out (resulting in distinct probability values for differential reduction in energy). Details of our simulation model and the reaction probabilities can be found in our earlier work [4,7].
Developing computational models for TRAIL induced activation more malignant phenotypes if not killed by repeated chemotherapy ("adaptive/acquired resistance"). Induction of rapid apoptotic death in cancer cells (such as by increasing TRAIL concentration), on the other hand, may increase healthy cell cytotoxicity. Extensive computer simulations complemented with theoretical analysis, carried out by us, have elucidated mechanisms for inducing cancer cell apoptosis in a manner that selectively activates cell death in cancer cells while keeping fractional cell killing and healthy cell cytotoxicity at a minimal level [3,4,7]. However, in our previous studies we did not consider generation of anti-apoptotic proteins during the course of apoptotic induction (such as Bcl2 protein synthesis under TRAIL induction), impact of which can be significant on fractional cell killing and generation of acquired resistance.
In this article, based on our ongoing effort in studying single cell biology of apoptosis induction by various chemotherapeutic agents, we report potential applications of stochastic simulations in developing effective cancer therapy. In a previous article, we studied the possibility of utilizing lower affinity TRAIL-like ligands in cancer chemotherapy [7]. Recently protein engineering has provided us with a specific set of affinity variant TRAIL molecules that have variable affinity for death and decoy receptors [15,16]. Here, based on recent experimental studies of apoptotic activation induced by affinity variant ligands [15,16], we carry out Monte Carlo simulation of apoptotic activation for two such altered ligands (DR5A and DR5B) in death receptor 5 (along with decoy receptor 2) expressing cancer cells [16]. In line with recent experimental studies, the extent of apoptotic activation (under induction of different ligands) in cancer cells was used as selective criteria for choosing the optimal ligand [15]. However, use of single cell approaches (Monte Carlo runs) also allows us to assess cell-to-cell variability in time-to-death and analyze the cumulative probability distribution in time-to-death (in deciding among the available ligands), as slow activation of apoptosis has been implicated in generation of adaptive resistance.

Methods
In the current Monte Carlo computational model, we plan to capture the stochastic dynamics of molecular reaction networks (inside a cell) perturbed by death ligands. Induction of death ligand (such as TRAIL) activates caspase 8 (initiator caspase) in the membrane proximal death module and thus acts as a trigger for apoptosis activation through type 1 and type 2 pathways [1] (such complex apoptotic networks are found in cells of higher organisms such as humans). Death ligand binding presumably induces some conformational changes in the death receptors leading to increased propensity to form oligomeric complexes that are capable of binding death adaptor proteins [17]. The death inducing signaling complex (DISC) thus formed could recruit pro-caspase 8 (along with anti-apoptotic cFLIP) and induce its activation depending on the extent of DISC formation and pro-caspase8/cFLIP ratio. In the type 1 pathway, caspase 8 directly activates effector caspases (caspase 3/7). In the type 2 pathway, a second initiator caspase 9 is needed for signal amplification, activation of which is controlled by pre-and post-mitochondrial regulatory mechanisms. Activation of caspase 3/7 (effector caspase) was taken as a downstream readout of apoptotic signaling activation.
We have developed kinetic Monte Carlo simulations (numerical simulation techniques) to solve underlying master equations (equations in probability space) that govern the stochastic dynamics of a biological system, such as intracellular reaction networks composed of different molecular species inside a biological cell. When simulating such intracellular reaction networks (such as the apoptotic signaling have been generated that could selectively induce apoptosis in various cancer cell lines [15] but one needs to find an optimal ligand for treating a given tumor. We simulated and compared the following ligands: natural TRAIL (high affinity for DR5 and DCR2), DR5B (TRAIL mutant variant with high affinity for DR5 but 400 fold lower affinity for DCR2) and DR5A (TRAIL mutant variant with high affinity for DR5 but 10 fold lower affinity for DCR2). In some types of cancer, death receptor 5 (DR5/TRAIL-R2) has been shown to be upregulated or its expression can be induced [16,29]. We chose the initial composition of death/decoy receptors to be DR5 ~ 50 molecules and DCR2 ~ 50 molecules (DCR2/DR5 ~ 1) on a 1.2 × 1.2 µm 2 surface area of cancer cells [7]. Death ligand DL=2 molecules was used to keep healthy cell cytotoxicity at a low level. In addition, XIAP and cFLIP levels are taken to be low (~ 10 nM), either due to low rate of expression or achieved by application of inhibitor molecules [13,30,31]. Monte Carlo simulation data was then analyzed to compare the three ligands in terms of their ability to induce death activation in cancer cells and to keep healthy cell cytotoxicity at a minimal level.
In Figure 1a, apoptotic cell death fraction and type 1 activation fraction are shown separately (for representative cancer cells for all three ligands). Even with significant expression level of decoy receptors (on cancer cells) selectivity for cancer cells was evident for all three ligands. Healthy cell apoptosis was found to be significantly lower than cell death in cancer cells (data not shown). Additional factors, such as lipid mediated effects and lower level of some of the apoptotic proteins may further inhibit the extent of healthy cell apoptosis. In addition, the type 1/type 2 phenotype choice was also found to be robust (type 1 dominant in cancer cells whereas type 2 dominant slow activation in healthy cells) across all ligands simulated. Type 2 activation was also found to be significant in cancer cells expressing high level of cFLIP and XIAP (when not inhibited by cFLIP/XIAP inhibitors as part of TRAIL combination therapy) [7,30]. Application of DR5B (selective ligand for DR5) could induce robust cancer cell death and showed a dominant type 1 activation phenotype.
In Figure 1b, we show the cumulative probability in time-to-death for the first 2/3 cells (fast responders) for all three TRAIL like ligands. The choice of optimal ligands (as well as the choice for optimal network modules to target) should address the issue of "adaptive resistance" arising from slow cell death and fractional cell killing. Under induction of TRAIL like ligands, slow cell death presumably allows stochastic generation of anti-apoptotic proteins (such as anti-apoptotic Bcl2 and/or XIAP through stochastic genetic regulations), which in turn may inhibit cell death resulting in a cellular bifurcation (survival/ death as shown in Figure 1). For natural TRAIL, high affinity binding to DCR2 leads to (in a fraction of cells): (i) lack of apoptosis inducible signalosome (DISC) formation or (ii) weak type 1 activation due to the combined effect of slow initiation in caspase 8 activation and inhibition of type 2 apoptosis through generation of anti-apoptotic Bcl2 molecules and survival of more than 10% of the cell population ( Figure 1). For DR5 selective ligands, XIAP and/or cFLIP inhibition allowed the type 1 activation to proceed even when ligand induced increase in Bcl2 levels prevent activation through the type 2 pathway [7,30]. Average and cellto-cell variability in time-to-death was found to be lowest for DR5B rendering it the ligand of choice. In some cases single cell analysis (such as cell-to-cell variability in time-to-death) and the extent of healthy cell cytotoxicity [7] may alter the choice for optimal ligands determined solely based on population average in cancer cell death.
For cancer cells equipped with low pro-caspase 8 level and/or high cFLIP level, but having significant expression levels of pro-and of the membrane proximal death signaling module is challenging due to intricacies of interactions among four different death and decoy receptors (DR4, DR5, DCR1, DCR2 [25] are TRAIL receptors), and, lack of sufficient experimental data on death/decoy receptor clustering. We have recently developed a Monte Carlo simulation model that could combine membrane proximal death/decoy receptor clustering and caspase 8 activation with intracellular apoptosis signaling [7]. In this work, we also allow stochastic generation of anti-apoptotic Bcl2 proteins under TRAIL induction (presumably through antiapoptotic pathways [26,27]), thus linking membrane proximal events, intracellular apoptotic activation and genetic regulatory mechanisms in one integrated mechanistic model. Synthesis step for anti-apoptotic Bcl2 proteins was carried out with constant probability in our current simulations. Additional complexities may arise because the kinetics of Bcl2 generation can be more intricate and/or it may vary depending on the ligand type.
Recently, protein engineering techniques have provided us with choice of ligands for targeting structurally/functionally similar surface receptors and other proteins in the apoptotic pathway. Affinity variant TRAIL-like ligands [15] and BH3 mimetics [28] are prominent examples of such altered ligands action of which can be simulated by varying the Pon/Poff constants for ligand binding/unbinding in Monte Carlo simulations. Thus it is possible to use in silico studies to determine an optimal ligand (or a combination of ligands), as well as the optimal network modules to target by that ligand. Here, we consider variable affinity binding to different death and decoy receptors by TRAIL like ligands [15] DR5A and DR5B. It is assumed that all three TRAIL molecules bind/unbind with the same probability to DR5 (P on =10 -2 and P off =10 -8 ; K A =10 9 M -1 ; approximation based on known nano-molar affinity of TRAIL for DR4/DR5 and SPR data in Ref. [15]).  [15]). We simulate a 1.2 µm × 1.2 µm × 1.2 µm cellular volume (60 × 60 × 60 lattice with lattice spacing ~ 20 nM) in such a manner that number of molecules included in simulation volume equals concentration (expressed in nanomolar). Total simulation time-scale was taken to be 10 8 MC steps (in one MC step each molecule is sampled once on average). Each run of our Monte Carlo simulation corresponds to apoptotic death activation in a single cell, thus, data obtained from many Monte Carlo runs can capture cell-to-cell stochastic variability (due to extrinsic and intrinsic noise) in signaling activation. Statistical analysis (average and standard deviation) was carried out using data obtained from many single cell runs (64 runs were used in current simulations; we also bootstrapped the data to carry out statistical variance estimation). Even though results shown here are based on a small simulation volume [7] (and the Monte Carlo simulation model has simplifying assumptions; e.g. TRAIL may impact the Rb pathway of cancer cells but is not considered here), it is expected that mechanistic insights provided by this study is robust.

Results
The goal of effective chemotherapy is to generate a rapid type 1 or type 2 (or even mixed) activation phenotype in cancer cells (to maximize cancer cell killing and minimize fractional cell killing mediated generation of chemoresistance) while reducing toxicity effects (minimize healthy cell toxicity) [7]. Several TRAIL like ligands anti-apoptotic Bcl2 family proteins and Apaf-1, one may use Bcl2 and XIAP inhibitors along with TRAIL to induce a rapid type 2 activation phenotype [7]. Thus, even when DCR/DR ratio is significant on a cancer cell (as indicated to be the case in certain cancer types [16,32,33]), it is possible to achieve selectivity in cancer cell targeting by use of affinity variant ligands and by correctly perturbing the apoptotic molecular network at the systems level (for generating an optimal signaling phenotype).

Discussion
Intense effort for more than a decade has led to generation of a large number of chemotherapeutic agents that directly or indirectly induces apoptotic elimination of cancer cells. Use of death ligands has generated considerable recent interest because it seems possible to induce apoptotic activation selectively for cancer cells, and, various signaling phenotype choices (type1/type2/type2-assisted-type1) are available to achieve rapid apoptotic activation (slow activation has been correlated with generation of chemoresistance). A key challenge is to decide an optimal ligand from several available death ligands (such as between TRAIL and FasL) but also from engineered ligands of similar kind that have recently been generated (such as affinity variants for TRAIL). It is reasonable to expect that the choice of an optimal death ligand (from TRAIL-like ligands) crucially depend on DR/DCR composition on the plasma membrane of cancer cells (also of potentially impacted healthy cells). Thus, depending on the tumor details (such as tumor subtype and/or tumor heterogeneity) and potential cytotoxicity for healthy cells one may choose the optimal ligand. In this work, we apply Monte Carlo simulations to elucidate mechanisms for generating a rapid type 1 activation phenotype in cancer cells having highly expressed DR5 and DCR2 under the induction of altered ligand DR5B (optimal choice of ligand). Hence, when DR5 expression is significant in a given tumor [16,34] (along with significant decoy receptor 2 expression), and, in healthy cells DR4 is the prevalent receptor (and/or the DR/DCR ratio is not high), then either of DR5A or DR5B might serve as the optimal ligand. For cholangiocarcinoma with highly expressed DR4, TRAIL variants having DR4 selectivity might serve as optimal choice for ligands [30]. Monte Carlo simulations can be carried out in a rapid manner, given proteomic data for membrane proximal death module, to assess the extent of selectivity in cancer cell targeting. One could simulate any combination of DR4/DR5/DCR1/DCR2 levels and estimate the extent of cell death (also long time survival/death bifurcation), probability distribution in time-to-death and type 1/type 2 phenotype choice at the level of single cells. It is reasonable to expect that such a single cell based in silico approach will allow us to derive an optimal strategy in cancer cell targeting.
Fluctuations in cell death outcome (fractional cell killing), originating from a complex interplay between stochastic intracellular/ genetic regulation and genetic/epigenetic variations (such as copy number variations), may underlie generation of resistant cancer cells ("adaptive resistance"). In many situations, slow activation of the death pathway will allow degradation of pro-apoptotic molecules that are being activated (under apoptosis induction) [12,35] and/or synthesis of anti-apoptotic proteins [5,25,36] resulting in survival of a fraction of cells (characterized by a death/survival bifurcation in long time cell fate choice) [2,12,37]. In the case of TRAIL (and its variants), ligand induction may induce expression of anti-apoptotic proteins such as Bcl2 and/or XIAP (in a stochastic manner) during the course of slow type 2 activation (~hours) and block apoptotic activation in a fraction of cells. Such bifurcation in cell fate at the level of single cells is of purely kinetic origin [18,19]. In this work, we showed that application of TRAIL variant DR5B along with a cFLIP and/or XIAP inhibitor results in robust type 1 apoptosis activation in cancer cells expressing DR5 and DCR2 (whereas for wild type TRAIL combined effect of slow caspase 8 generation and high level of Bcl2 induction stochastically inhibits apoptosis in a fraction of single cells). To minimize fractional cell killing and to overcome generation of resistant cancer cells one may need to apply such optimal ligands with combinatorial treatment options [13] and determine an optimal network module for inducing apoptosis.
Even though the study in this work has focused on TRAIL induced apoptosis, the approach taken here is more generally applicable. In certain cancers, one may want to initiate apoptosis by activating Bax either directly by inhibiting Bcl2-like proteins (using BH3 mimetics) or indirectly using tyrosine kinase inhibitors [3,13]. We expect that selectivity in cancer cell targeting is possible to achieve by use of affinity variant ligand/inhibitor molecules and by correctly perturbing the apoptotic molecular network at the systems level (for generating an optimal signaling phenotype). Choice of optimal ligands in cancer chemotherapy is intricately linked with elucidation of inherent vulnerability in cancer cells (such as due to higher DR/ DCR ratio, higher level of BH3 only proteins, altered lipid composition and receptor occupancy in raft domains [38] among others) and finding optimal network modules to target. Targeting different network modules in the apoptotic pathway may result in distinct signaling phenotypes. Theoretical approaches (such as stochastic/ probabilistic modeling) provide mechanistic insights into generation of different signaling phenotypes in death pathway activation (such as combinations of slow/fast with type1/type2) and fractional cell killing. In terms of dynamical systems, one aims to induce monostable death in cancer cells and monostable survival in healthy cells [35] (but might not be achievable due to the presence of extrinsic and intrinsic noise). It is expected that chemotherapeutic perturbation should be such that it would induce only slow activation (with large cell-to-cell variability) phenotype in healthy cells resulting in protection of most cells. Monte Carlo based in silico studies can rapidly simulate various scenarios of death pathway activation leading to identification of optimal network modules (to target by apoptotic agents) for a given cancer subtype. One may need to consider an expanded signaling network that will include the possibility of inducing apoptotic death in an indirect manner (such as through oncogenic kinase inhibition). In addition, one may utilize dynamic mechanisms (such as through inducing receptor-lipid clustering [39,40] in certain cancer cells) for achieving selectivity in targeting of cancer cells.
The scope of using affinity variant ligands (in cancer chemotherapy) is expected to be broader than that elucidated in this study. For example, targeting functionally similar proteins such as Bcl2/Mcl1/ BclxL by a single ligand inhibitor might pose challenges as the binding affinity may vary among members of the group. For example, targeting Bcl2 like anti-apoptotic proteins by small molecule inhibitors with high selectivity (such as by Bcl2-selective AbT-737) may present with resistance through over-expression of other members of functionally similar proteins; hence, cancer cells overexpressing Mcl1 may induce resistance to AbT-737 targeting [41]. In such a scenario one may have to choose a broadly neutralizing inhibitor or combine a set of different ligands. In the case of death receptor mediated apoptosis, TRAIL ligands have been generated that are (i) weaker agonist than natural TRAIL [26] or (ii) have variable affinity for different death and decoy receptors (DR4/DR5/DCR1/DCR2) [15,16]. Monte Carlo simulations can be carried out to determine optimal affinity values (for each of the death and decoy receptors) and concentration to be used for the TRAIL variant. One may use protein engineering and in silico modeling to determine the optimal ligands for a specific cancer subtype.
In addition to its crucial role in finding optimal strategy in cancer cell targeting, Monte Carlo based systems level elucidations may prove to be essential in linking large-scale data relevant for cancer detection, sub-type classification and therapy. Genomic instability is a characteristic feature of malignant cells and prior bioinformatic analyses have identified genetic changes associated with chemoresistant phenotypes. NGS methods now allow access to genomic/epigenomic/ transcriptomic data (for a given tumor) at the level of single cells [42]. On the other hand, clinical data has been in use for cancer subtype classification (such as growth factor and hormone receptor expression based classifications for certain cancers [43]) and therapeutic decisions. Systems level mechanistic studies, such as one that is carried out in this work, is expected to bridge the gap between large-scale genomic (Omics data in general) and patient data thereby providing rigorous methods for establishing cancer subtype classification. Monte Carlo based simulations (based on genomic/proteomic data for a given tumor) that integrate membrane proximal events, intracellular signaling and genomic regulatory mechanisms should provide insight into disease subtype classification (such as one based on DR/DCR composition of plasma membranes of cancer cells, e.g., DR4 -DR5 + DCR1 -DCR2 + , DR4 + DR5 + DCR1 -DCR2 -). However, the scope of such mechanistic elucidations in bridging large-scale data seems to be much broader than discussed here. For example, genomic data regarding gene amplification and mutations can be connected to autoantibody (bio-marker) data in patients [44] through mechanistic elucidations of humoral immune response.
Tumor heterogeneity has posed a considerable challenge in treating various types of cancer. Even though analyzing biopsy samples from different regions of a tumor can capture such heterogeneity to some extent [45], single cell genomic/proteomic analysis is able to reveal true complexities of a heterogeneous tumor [42]. Next challenge is to characterize a given tumor in terms of known subtypes (or even as a superposition of different subtypes). As discussed in the previous paragraph, integration of genomic data analysis (based on bioinformatics approaches) with systems biology based mechanistic modeling (such as Monte Carlo simulations) might be helpful in disease subtype classification. Additional complexities may arise from factors such as the tumor microenvironment in a given patient, analysis of which is also amenable to Monte Carlo based computer simulations. Once we have achieved a rigorous way of characterizing a given tumor (subtype or combination of subtypes based on omics data from primary cells), one may use a combination of different strategies to treat a given tumor in a patient (personalized approach). It is expected that elimination of cancer cells under apoptosis inducing therapy will help reduce the immunosuppressive tumor microenvironment and/ or molecular antigens released from cancer cells that had undergone apoptotic death will induce strong adaptive immune response against tumor cells [28]. High resolution in-vivo imaging experiments that have recently been developed [46] may provide insights into tumor microenvironment changes under apoptosis inducing targeted therapy. One may also consider dynamic changes in antibody production (to tumor antigens) [44] to assess the impact of cancer cell death on tumor microenvironment. In some cases, combination of chemotherapeutic agents with immunotherapy (including tumor vaccines based on individual genome sequencing) may lead to an optimal therapy for cancer [47,48].
We have addressed the issue of finding optimal network modules to target (by chemotherapy) as well as the optimal ligands needed for a given tumor (in a specific patient). One may also consider other factors such as the optimal dosing strategy when developing optimal treatment options. Finding an optimal strategy when several options are available would require rigorous analysis of cancer celldrug interaction and apoptosis induction at the level of single cells (including potential for generating resistance through fractional cell killing). Computer simulations can be carried out in a rapid manner for various chemotherapeutic strategies and also for various cellular parameters (representing different inherent/initial conditions of cancer cells). We are in the process of developing a scoring based approach to search for a few optimal strategies in targeting a given tumor [7]. We expect that Monte Carlo simulations and scoring based single cell data analysis, combined with bioinformatic analysis of genomic data, will help develop disease subtype classification schemes and effective personalized cancer therapy.