Discovery of Drug Synergies in Gastric Cancer Cells Predicted by Logical Modeling

Discovery of efficient anti-cancer drug combinations is a major challenge, since experimental testing of all possible combinations is clearly impossible. Recent efforts to computationally predict drug combination responses retain this experimental search space, as model definitions typically rely on extensive drug perturbation data. We developed a dynamical model representing a cell fate decision network in the AGS gastric cancer cell line, relying on background knowledge extracted from literature and databases. We defined a set of logical equations recapitulating AGS data observed in cells in their baseline proliferative state. Using the modeling software GINsim, model reduction and simulation compression techniques were applied to cope with the vast state space of large logical models and enable simulations of pairwise applications of specific signaling inhibitory chemical substances. Our simulations predicted synergistic growth inhibitory action of five combinations from a total of 21 possible pairs. Four of the predicted synergies were confirmed in AGS cell growth real-time assays, including known effects of combined MEK-AKT or MEK-PI3K inhibitions, along with novel synergistic effects of combined TAK1-AKT or TAK1-PI3K inhibitions. Our strategy reduces the dependence on a priori drug perturbation experimentation for well-characterized signaling networks, by demonstrating that a model predictive of combinatorial drug effects can be inferred from background knowledge on unperturbed and proliferating cancer cells. Our modeling approach can thus contribute to preclinical discovery of efficient anticancer drug combinations, and thereby to development of strategies to tailor treatment to individual cancer patients.


Section I -Signaling network building and manual curation 1.1 Identification of pathways to include in the manually curated model
The AGS gastric adenocarcinoma cell line was chosen because AGS cells harbor mutations in numerous genes encoding proteins that are central in signaling pathways known to be affected in gastric adenocarcinomas: MAPK-, PI3K-, Wnt/β-catenin and NFκB-pathways [55,56,[69][70][71]. These four signaling pathways were chosen as a focus for our signal transduction model. Mutations reported in AGS cells can be found in the databases listed in Table 1.

Manual curation of the signal transduction model 1.2.1 Inclusion of signaling components
In order to construct a model depicting current knowledge on the MAPK, PI3K, Wnt/β-catenin and NF-κB pathways we used publicly available models in databases as guidance for inclusion of signaling components and interactions. Databases consulted are shown seen in Table 2.

CellDesigner model
The collected pathway information was encoded in a CellDesigner SBML model [84,85] where edges are annotated with literature references from which the interactions are substantiated, see Figure 1. All proteins are annotated with the following information: Uniprot ID, gene name, presence at transcript level in AGS, and presence of mutations in AGS cells as reported in the databases listed in Table 1. Yellow boxes denote signal transduction components, and red boxes denote proteins whose gene harbor known mutations in AGS cells. Components identified by black circles mark the components targeted with specific inhibitors. The two diamond-shaped boxes "Prosurvival" and "Antisurvival" correspond to the two main readouts of the model. The curated model and its annotations are available in SBML format (S1 Dataset).

Exclusion of signaling components
The curated model was trimmed so that only nodes that are regulated by other nodes also present in the model were retained. Thus all model nodes are part of feedback loops, except for the two output nodes and some of their close upstream components (BAD, BAX, BCL2, Caspase3/7, Caspase8, Caspase9, CFLAR, CytochromeC and ITCH).

Final SBML signaling model
The final SBML signaling network consisted of 75 signaling components and two phenotypic output nodes, and 149 interactions, and the full SBML model with annotations can be found in S1 Dataset.

Default logical functions
We formulated a set of Boolean equations to formalize the evolution rules of each component. Boolean equations were formulated by the following general formula: A signaling component is active if any of its regulatory activators were active, while at the same time none of its regulatory inhibitors were active. For activators B, C and D and inhibitors E, F and G of the node A, this translates to: A = (B or C or D) and not (E or F or G) Note that several nodes only have activating regulators, and several nodes only have inhibitory regulators.
As described in the main text a few exceptions to this rule were made. The logical formulae of all nodes can be found in S5 Table. 1.3.2 Exceptions to the default logical functions

Adjustments to comply with knowledge on protein activity regulation
The general Boolean equation stated above does not always comply with the biological understanding of protein activity regulation. Two Boolean rules were changed because of mechanistic knowledge on the protein activity regulation: 1) The equation describing the activation of the β-catenin destruction complex, which activates β-TrCP, resulting in degradation of β-catenin, and 2) the equation describing β-catenin activity.
For the β-catenin destruction complex to be active all subunits must be available, and so this multi-protein complex was said to depend on all three of Axin, GSK3 and CK1 (note that APC was excluded from the model as it wasn't regulated by any component present in the model, and also note that Dvl was modelled as an intermediary between the Frizzled receptor and GSK3 [86], and not as part of the destruction complex). The logical function is: betaTrCP = Axin and GSK3 and CK1 instead of betaTrCP = Axin or GSK3 or CK1 The activity of β-catenin should further be governed by β-TrCP in a way so that the absence of β-TrCP causes β-catenin to be active. In the model IKKA is also positively regulating β-catenin activity, and the activity of β-catenin can, if IKKA is active, be sustained even when β-TrCP is active, since IKKA is phosphorylating and protecting β-catenin from destruction [64]. Thus, the equation describing activity of β-catenin was changed from the general formula: betacatenin = IKKA or not betaTrCP instead of betacatenin = IKKA and not betaTrCP    Table. Literature reports of protein activity levels were found for 46 of 75 model signal transduction components.

Adjustments to comply with reported observations on protein activity in AGS cells
We reasoned that whenever the information from several scientific publications consistently implied that a specific protein should be active or inactive in proliferating AGS cells (for example AKT was consistently found to be active), the observed stable state attractor should reflect this (the component AKT value should equal 1). For 21 of the 46 components several, independent and consistent reports were found, and these activity state values were used for guiding and validating the Boolean equations. Of the 25 component observations that were less well substantiated, ten were found to be in conflict with the predicted state in the attractor. We consider compliance with the 21 well documented components to be of highest importance and therefore concluded that the stable state of the model adequately matches experimental observations of protein activities in actively growing AGS cells.
In searching PubMed for relevant articles, we initially searched for all AGS cell papers reporting Western blots published after the year 2000. This search resulted in 376 papers (the literature survey was performed in June-September 2013). For model proteins with none or only few observations in these initially retrieved papers, a second round of PubMed searches was done retrieving all papers that could be found with the keywords "AGS [missing protein]". This search was not restricted by the time of publication. In addition, Google was used to search for papers with keywords "AGS [missing protein]", identifying a few articles that were not returned by PubMed's internal search engine (even though the paper was indexed by PubMed).
For all observations their suitability was assessed according to the quality of the reported Western blots (subjective perception), specifications of whether cells were in complete medium or starved, whether phosphorylation events or only total protein levels were reported, and if the observation agreed with observations in other papers. An assessment of whether a protein was active or not was then done and activity levels were dichotomized into the two states ON or OFF. See S2 Table.

Literature review of AGS signal transduction perturbation
As an additional test of the validity of the model, the asymptotic behavior of the perturbed model was verified by comparing it with the behavior observed in experimental perturbations. We retrieved information from 56 published experiments where specific signaling proteins had been activated or inhibited (S6 Table), and recorded the reported effects of these perturbations on signaling proteins. The experimentally observed effects were then compared with the effects found for the corresponding nodes in model simulated perturbations. The model accurately reproduced 21 published observations. None of the 56 observations were contradicted, meaning that the model never proposed the result opposite to findings in published reports. For 35 observations, the reported change in activity was not captured by the logical model. We consider this to be related to the likelihood that not all activity changes observable on a continuous scale by methods like Western blot will translate into a qualitative signaling change.
In order to capture signal transduction dynamics, the literature was surveyed for reports on AGS cell perturbation experiments where the target was a node in our curated model, and the observed perturbation effect pertained to another node in our model. The papers retrieved from the following PubMed-searches, "AGS siRNA" (163 papers), "AGS [relevant specific inhibitor]", and "AGS inhibitor", were assessed according to the criteria specified above and resulted in 56 observations spanning 20 unique publications that were of relevance to the model (see appendix C).

Model validation
The predicted stable state of the model fit data from literature, as seen in Table 4. This list of 21 components is a summary of the observations shown in S2 Table, with only the observations of highest confidence included. exceptions for a few of the less confident literature observations, but as the quality and/or consistency across independent reports for these observations was substantially lower than the observations in Table 4 these discrepancies were not analyzed further.

Section IV -Model reduction and dual interaction from p38alpha to Antisurvival 1.7 Model reduction
Using the built-in model reduction tool of GINsim a logical model can be reduced to a model consisting of fewer nodes, while stable states of the original logical model are preserved. During model reduction, all logical equations are replaced by logical parameters, which amounts to having the logical rules in a lengthy normal disjunctive form. These are provided in S7 Table for the reduced logical model depicted in Figure 3 in the manuscript, along with the reduced model in S3 Dataset. S7 Table was not used for drug synergy analysis, but can provide details to the interactions in model. Compressed expressions of the logical parameters are provided in S8 Table. For each interaction S7 Table specifies the context determining if the interaction is positive, neutral or negative. The context is a discrete vector containing the activity values of all nodes contributing to the nature of the interaction. We use it here to describe the context-dependent nature of the dual interaction from p38alpha to Antisurvival.

Dual interaction from p38alpha to Antisurvival
After model reduction the interaction from node 'p38alpha' to 'Antisurvival' is shown as 'dual' (Figure 3, Manuscript and S7 Table), where p38alpha can increase the value of 'Antisurvival' in some contexts, and decrease it in other contexts. See S7 Table for a full list of contexts that affects the interaction from p38alpha to Antisurvival. From S7 Table it can be seen that two scenarios can hypothetically happen: (1) p38alpha inhibition can decrease Antisurvival if AKT is inactive (2) p38alpha inhibition can increase Antisurvival if AKT is active Scenario (1) is never observed in our simulations. For (1) to happen one of the drugs in a drug combination would have to cause p38alpha to be activated as p38alpha is inactive in the unperturbed case, and AKT to become inactivated (AKT is active in the unperturbed case). Then if p38alpha inhibition is the second drug in the drug pair we could observe behavior (1). However, among the drug inhibitions we simulated no single drug has the ability to both activate p38alpha and inactivate AKT, and therefore (1) is never observed.

Section V -Simulations of FOXO knockout
From model simulations it was hypothesized that FOXO is important to the predicted and validated synergies. Model simulations where FOXO was forced to be inactive (i.e. knock-out simulations) showed that the synergies of TAK1 and PI3K, and AKT and MEK, and MEK and PI3K were all dependent on FOXO, while the synergy of AKT and MEK was only partially dependent on FOXO (the synergy was weakened but not abolished when FOXO knock-out was simulated). The results from simulations can be seen in Table 5.  -AKT  3  1  2  --AKT  OFF  3  1  2  -p38alpha  3  0  3  --p38alpha  OFF  3  0  3  -TAK1  3  0  3  --TAK1  OFF  3  0  3  -Unperturbed -3  0  3 --In addition to the importance of FOXO activation, we found AGS cell growth to be highly dependent on β-catenin, since administration of the β-catenin inhibitor PKF118-310 even at low concentrations profoundly impeded cell growth in comparison with its effect reported for other cell lines [87][88][89]. Our AGS logical model simulations correlate well with these experimental observations since the growth output exhibited high dependence on the state of the node representing β-catenin, yet model simulations predicted that this trait was not dependent on β-catenin induced activation of FOXO. We thus searched for inhibitor combinations targeting β-catenin and any other node accounted for by the model, to suggest new potent drug combinations. While we found that none of the pairwise inhibitions of components included in the model were able to inactivate β-catenin and simultaneously activate FOXO, model simulations revealed that the strong anti-proliferative state induced by β-catenin inhibition could be further strengthened by simulating activation of FOXO (output values: Prosurvival = 0, Antisurvival = 2, contrasted with Prosurvival = 0, Antisurvival = 1 for β-catenin inhibition alone), emphasizing the central position of FOXO regulation in AGS growth..

Section VI -Cell growth experiments 1.9 Chemical inhibitors
Seven chemical inhibitors were chosen for cell growth experiments. Criteria that guided the selection of inhibitors include: • Chemical properties similar to clinically approved drugs (Lipinski's rule of five) • Well characterized in in-vitro kinase screens, in order to know of any off-target effects, and preferentially choose inhibitors with few off-target effects, following the recommendation by the MRC lab whenever possible (http://www.ppu.mrc.ac.uk/).The inhibitors used for perturbation experiments are shown in Table 6.

Determination of GI50 concentrations
For each drug a dose-response profile was initially determined, in order to find a dose that inhibited growth by 50%, and which was the dose used in combination experiments. If there was no knowledge on the effective concentration range to be used a 10x dilution series was performed initially, before a 2x dilution series could be done. If the GI50 or a proxy could be found in public databases (Genomics of Drug Sensitivity in Cancer) or deduced from literature (either for AGS cells or other cells) the dilution series would start with a 5x or 2x dilution series.

.2 β-catenin inhibitor
The β-catenin inhibitor PKF118-310 was highly potent in AGS cells, and proved to be very difficult to work with in experiments, as there was either an all-or-nothing response for a given concentration. 150 nM was used as an estimate of the GI50, even though this would typically give a full effect. (At 100 nM there would typically be seen no effect). In other cells the GI50 value has been found in the range 0.31-3.48 µM [87][88][89], indicating that AGS cells are indeed very sensitive to this inhibitor.

MEK inhibitor
The MEK inhibitor PD0325901 was found to have a GI50 of roughly 35 nM. This fits data from the Genomics of Drug Sensitivity in Cancer database, where the GI50 for AGS cells has been found to be 0.013-0.032 µM.
1.9.1.2.4 TAK1 inhibitor The TAK1 inhibitor (5Z)-7-oxozeaenol was found to have a GI50 of 0.5 µM in AGS cells. In literature the GI50 for other cell lines has been reported in the range 0.02-60 µM [91,92]. In KRAS-dependent colon cancer cell lines concentrations in the range 0.625-1.25 µM of (5Z)-7-oxozeaenol has been found to promote apoptosis [92], and in another work AGS cells have been described as KRAS dependent [65], and the GI50 found in our experiments are thus in agreement with the few reports of the effects of (5Z)-7oxozeaenol on proliferation in cell line experiments.

AKT inhibitor
The AKT inhibitor AKTi-1,2 (AKT inhibitor VIII) was found to have a GI50 in AGS cells of 10 µM. This fits data from the Genomics of Drug Sensitivity in Cancer database, where the GI50 for AGS cells is estimated to be 7.6 µM.
1.9.1.2.6 GSK3 and p38alpha inhibitors For the two inhibitors CT99021 and BIRB0796 calculated GI50 values are included for reference. This GI50 value is however probably not biologically meaningful as the concentration of each inhibitors is very high compared to concentrations that are typically used to inhibit their intended target. This also coincides with data from the Genomics of Drug Sensitivity in Cancer database, where the GI50 for these inhibitors is estimated by extrapolation to be > 200 µM for AGS cells.
1.9.1.2.7 GI50 vs time The GI50 varied with time, and so an initial experiment was performed to find a time point at which the GI50 could be estimated. The GI50 stabilized for all inhibitors and combinations of inhibitors within 48 hours after addition of inhibitors, and so this was chosen as the time point for further analysis (both single drug dose response curves, and for drug combination experiments). All growth curves were still recorded in real time, and the 48h time point was used for calculating GI50s and combinatorial indexes (see below). Figure 16 and Figure 17 show two examples of GI50 estimates versus time.  . Initially, a doseresponse experiment of each inhibitor was performed to find a concentration of each inhibitor that inhibited growth of cells by half compared to the control (no inhibitor). This dose of the drug is the GI50 (growth inhibition 50%), and for combination experiments a dose close to but less than the GI50 was sought. Since the exact effect of all drugs will vary slightly between biological replicates, experiments were set up so that both single drugs and their combinations would be tested within one experiment, i.e. within the same biological replicate.

Experimental design
All single drugs were tested in three doses within one experiment, in a 2-fold dilution series from the GI50 dose. Similarly the combinations were tested in three doses within one experiment. The highest dose was obtained by mixing the two inhibitors with equal volumes from the GI50 concentration solutions, and thus the concentration of each inhibitor in the combination was half of the GI50. Then a 2-fold dilution series was made so that the doses of each inhibitor were a quarter and an eight of their GI50 doses, which would give three doses where synergy could be assessed according to the Loewe definition. Note that the highest dose of inhibitors in a combination were created from the same aliquots as the single inhibitors, thus ensuring that the concentration in a combination was halved, even if the exact concentration of the single inhibitor deviated slightly from the GI50 concentration.

Combination indexes
A combination index (CI) was also calculated for all combinatorial experiments based on the analysis of Chou and Talalay [29], where the CompuSyn software was used for calculating the CI values [68]. The combinatorial index is based on a mechanistic approach to synergy assessment. A median-effect plot is made from specific dose-response data, where the potency and shape of the dose-response is determined, and thus the entire dose-effect curve is predicted. This analysis also takes into account the steepness of the dose-response curves when the effect of combining halves is analyzed (http://www.combosyn.com/) [68,93]. For analyzing the combination index the effect on growth 48 hours after adding inhibitors were used. The GI50 was found to vary with time, but had reached a stable level at this time, as described in the paragraph "GI50 vs time". Note that a CI value can only be reliably computed if a dose-response relationship to an inhibitor can be found, which made CI value calculations impossible for drugs that target proteins that are inactive in proliferating AGS cells (CT99021 targeting GSK3 and BIRB0796 targeting p38).  Figure 18 to Figure 21 show the four synergies found. For all figures the combination of two drugs has half the concentration of each single drug in the combination. If the effect of a combination is more profound than the best-performing single inhibitor the effect is synergistic. The plots show the reported Cell Index, which is proportional to the number of cells, on the y-axis, and the time in hours on the x-axis.

Appendix C
Literature references for the AGS signal transduction perturbation experiments (for full reference details, see Supplementary references next page): [66, 90,159,162,164,166,172,176,178,183,185,189,194,195,197,205,219,224,225].