In order to perform simulations of cellular behavior upon drug perturbations, the prior knowledge network was converted into a logical model in the software suite GINsim, version 3.0.0b (Gonzalez et al., 2006). In the logical model, components of the regulatory network (genes, proteins, complexes, chemical components/therapeutic agents, and phenotypes) are represented by nodes (connection points). Activating or inhibiting interactions are represented by edges. Apart from ‘Oxaliplatin’, ‘Fluorouracil’ and ‘Overall_phenotype’, all nodes are connected to a minimum of two other nodes via directed edges (a minimum of one incoming and one outgoing edge). ‘Oxaliplatin’ and ‘Fluorouracil’ act as input nodes, i.e., nodes fixed at a specified activity level, without incoming interactions and regulatory rules. ‘Overall_phenotype’ works as an output node, providing rapid access to the simulated results. The activity state (active/inactive) of non-input nodes is precisely defined by their regulatory nodes and a set of logical operators (AND (&), OR (|), NOT (!)), assembled into a logical rule that describes the conditions for activity of the target node. The overall model fate is determined by the output node ‘Overall_phenotype’, which can take on one of three levels (1, 2, 3). Apart from the output node, all nodes are Boolean. A complete summary of nodes and associated logical rules of the models is presented in Supplementary Table S3 (Supplementary File S1: Results, figures, and tables). The initial model contained 81 nodes and 142 edges.

In the initial unadjusted model, the activity of all non-input nodes except DSB, representing double stranded DNA breaks, and Overall_phenotype followed a standard default (equation 1), Target node = (Activator 1 | Activator 2 | etc.) &! (Inhibitor 1 | Inhibitor 2 | etc.)! (1) i.e., for a node to be active, at least one of its activating regulators and none of its inhibiting regulators had to be active. We constructed the model so that the absence of perturbations, as well as perturbation of single nodes, resulted in an output corresponding to a proliferating cell phenotype (‘Overall_phenotype’ at its maximum level, 3). To meet the assumption of no overall effect of single node perturbations, subsequent model adjustments were guided by reference data (drug combination effects) generated from low to moderate drug concentrations, for which no significant effect on apoptosis was observed upon single-drug administration.

Following initial analysis, the model was adjusted in order to improve its ability to correctly classify synergistic and non-synergistic drug combinations observed in HCT-116 cells in vitro (Folkesson et al., 2020). The model was adjusted by stepwise implementation of one or several of the following changes.

- Changing of logical rule(s) for existing nodes.

Each step was counted to keep track of the number of changes. Since the addition of a new node to a self-sustaining regulatory network involves addition of a minimum of two edges as well as modifications of at least one logical formula, we count these related changes together as one model adjustment. Likewise, since adding an edge between two existing nodes requires updating of the logical formula of the target node, these changes are counted as one model adjustment. In general, when adjusting a logical formula because of edge addition, the formula was first defined by the standardized form above (equation 1) (Mendoza and Xenarios, 2006). Any additional changes to the formula were counted as separate adjustments. Changes were made in the order in which they are listed in Supplementary Table S4 (Supplementary File S1: Results, figures, and tables). Overall, the order of model adjustments was motivated by the synergy scoring order of the combinations of interest (Supplementary Table S5, Supplementary File S1: Results, figures, and tables). The adjusted model consisted of 85 nodes and 158 edges. A complete summary of nodes and associated logical rules of the models (initial/adjusted) is presented in Supplementary Table S3 (Supplementary File S1: Results, figures, and tables). GINsim model definition files are available in the Figshare repository (see Data availability).

An overall model output was defined based on upstream signaling component activities, and ranged from 1 to 3, where lower numbers correspond to growth-inhibited phenotypes and higher numbers correspond to growing phenotypes. To determine the level of the model output (1, 2, or 3), as well as the activity of non-input nodes at baseline and upon perturbation, stable states analysis was performed using the Compute stable states algorithm implemented in GINsim. By making use of so-called decision diagrams for representation of the Boolean functions of the model, the Compute stable states algorithm enables computation of stable state(s) from an initial given condition without the necessity to construct a state transition graph describing the dynamical behavior of the model (Naldi et al., 2007). This allows assessment of stable states also when working with larger networks where the construction and simulation of state transition graphs is not possible due to the size of the resulting graph. Due to the size of the network, all simulations were run using synchronous updating. In cases of multiple stable states, the Assess attractor reachability tool was used to estimate the probability of reaching each of the suggested stable states. Here, the evaluation of reachability is based on stochastic algorithms for evaluation of probabilistic outcome(s) resulting from repeated analysis of randomly sampled inputs (Naldi et al., 2018). Multiple such algorithms are implemented in GINsim and in this work, attractor reachability simulations were performed using the historically commonly featured MonteCarlo algorithm. The simulations were performed with 1,000 runs and a max depth of 1,000 as running parameters. Also here, the simulations were run using synchronous updating of the network. The stable state associated with the highest probability was selected as the main stable state. The level of the model output given by the stable state analysis was used for evaluation of combination effects in silico. As the baseline (unperturbed) model was associated with the highest output level (3) and since single-node perturbation did not have a reducing effect on the overall output level, perturbation of combinations of nodes which resulted in an output level ≤ 2 were, in the case of pairwise combinations, considered synergistic. Mathematically, this quantification of combination effects in simulations is analogous to synergy quantification using the Highest Single Agent (HSA) reference model, one of several mathematical reference models often employed for cell laboratory experiments (Lehár et al., 2007). Analyzed model perturbations and their in vitro counterparts (Folkesson et al., 2020) are listed in Supplementary Table S6 (Supplementary File S1: Results, figures, and tables). Nodes targeted by model perturbations adhere exclusively to the list of primary targets reported for each of the drugs that were tested in vitro (Supplementary Table S1; Supplementary File S1: Results, figures, and tables).

For in vitro validation of predictions made by the computational model, HCT-116 cells (CVCL_0291), obtained from the US National Cancer Institute (NCI), were used. Cells were cultivated in 1X RPMI-1640 medium (Thermo Fisher Scientific) supplemented with 10% fetal bovine serum (FBS, Sigma Aldrich), 2 mM L-Glutamine (Sigma Aldrich) and 100 U/mL Penicillin–Streptomycin (Thermo Fisher Scientific) and maintained at 37°C with 5% CO₂ and 80% relative humidity. Passaging was performed according to in-house protocols as described in (Folkesson et al., 2020). Cells were used at passage numbers <20. Drugs used in the validation screen were olaparib (OLA), oxaliplatin (OXA), PI-103 (PI) (Selleckchem), PD0325901 (PD), 5-fluorouracil (5FU) (Sigma Aldrich) and 5Z-7-Oxozeaenol (5Z) (Enzo Life Sciences). The CellTiter-Glo 2.0 Viability Assay (Promega) was used for viability assessment of cells in the screen. A complete list of reagents and other material can be found in Supplementary Table SM1 (Supplementary File S2: Materials and Methods).

For high-throughput validation screening of third-order combinations PI + PD + OXA, PI + PD+5FU, OLA + PI + PD, 5FU + OLA + PD, 5FU + OLA + PI and 5FU + OLA+5Z, HCT-116 cells were robotically seeded (Biomek FX^P Laboratory Automation Workstation; Beckman Coulter) with 35 µL complete growth medium in 384-well black tissue culture-treated plates (Corning) at a density of 1,200 cells/well. Following seeding, plates were shaken (1,600 rpm, 30s) to ensure uniform distribution of cells. The cells were then incubated for 24 h (37 °C with 5% CO2 and 80% relative humidity), whereafter drugs were added. Drugs (single, pairwise, and triple combinations), as well as positive (digitonin, staurosporine) and negative controls (DMSO) were added using a Tecan D300e system. The DMSO concentration in the culture medium never exceeded 0.5%. All treatment conditions are listed in Supplementary Table SM2 (Supplementary File S2: Materials and Methods). Details regarding dose selection can be found in Selection of doses for third-order validation screening (Supplementary File S1: Results, figures, and tables). Each condition (drug/combination, concentration) was tested in four technical replicates. Drug addition was followed by shaking of the plates (1,000 rpm, 30 s) and incubation at 37°C with 5% CO₂ and 80% relative humidity. Following 48 h of treatment, cell viability was assessed using CellTiter-Glo 2.0, as described in (Folkesson et al., 2020). Cells were also imaged (2 sites per well) at 0, 24 and 48 h post drug addition. Imaging at 0 and 48 h was performed prior to drug addition and viability measurement, respectively. A SpectraMax i3x plate reader was used for all readouts. Two biological replicates of the high-throughput validation screen were performed.

For validation screening of third-order drug combinations OXA+5FU + PD and OXA+5FU + PI, and fourth-order combination OXA+5FU + PI + PD, HCT-116 cells were seeded with 90 µL complete growth medium in 96-well black tissue culture-treated plates (Corning) at a density of 4,800 cells/well. Seeding was followed by shaking of the plates (750 rpm, 30 s). The cells were then incubated for 24 h (37°C with 5% CO₂ and 80% relative humidity), whereafter drugs (single, pairwise, triple and four-way combinations) and a negative control (DMSO) were added manually. All treatment conditions are listed in Supplementary Table SM3 (Supplementary File S2: Materials and Methods). Each condition (drug/combination, concentration) was tested in two technical replicates. Drug addition was followed by shaking of the plates (750 rpm, 60 s) and incubation at 37°C with 5% CO₂ and 80% relative humidity. Following 48 h of treatment, cell viability was assessed using CellTiter-Glo 2.0. A POLARstar Omega plate reader was used for the readout. Three biological replicates were performed.

Confluency (high-throughput screen) was estimated by re-analyzing brightfield images using the SoftMax Pro 6.5.1 software. For each well, the percentage of covered area was estimated using settings specified in Supplementary Table SM4 (Supplementary File S2: Materials and Methods). For both confluency and viability (CellTiter-Glo 2.0), treatment effects were normalized to the plate-internal vehicle control (DMSO) and reported as average ± standard deviation. Pearson’s correlation coefficient has been used to quantify the association between variables (high-throughput screen) and biological replicates (high- and low-throughput screen).

The highest single agent reference model (Lehár et al., 2007) was used for evaluation of combination effects in the validation screens. HSA scores were computed by integrating the HSA excess value per data point to generate an overall HSA score for each combination across the tested doses. The latter was done to account for the fact that while effect size and synergy resulting from perturbations in vitro to some extent was dose-dependent, the corresponding perturbations in silico were logical to their nature. An integrated HSA excess score of 0 was used as cut-off, and combinations resulting in an integrated score < 0 and > 0 were deemed synergistic and antagonistic, respectively. Hence, in our in vitro experiments, synergy was called when the overall effect of the combination AB (for pairwise combinations), ABC (for triple combinations) and ABCD (for the four-way combination) was larger than the effect of the strongest underlying effector component (powerset), i.e., A and B for AB; A, B, C, AB, AC and BC for ABC; A, B, C, D, AB, AC, AD, BC, BD, CD, ABC, ABD, ACD, BCD for ABCD. Viability data (48 h) were used as the main data for assessment of combination effects. Confluency data were used for supplementary purposes, with combination effects evaluated at 24 h and 48 h. For each readout, combination effects were evaluated based on the average response of the biological replicates. In addition, for all tested higher-order combinations the statistical difference in relative viability between the highest-order combination versus the most effect underlying component (per dose) was calculated using a two-sided Student’s t-test. For each combination (third- or fourth-order) this was done per dose step and involved all technical replicates from all biological replicates.

For the high-throughput screen, intra-experiment reproducibility was assessed for both viability and confluency data by computing the Pearson correlation between the two biological replicates. The correlation coefficients were 0.99 and 0.96 for viability and confluency, respectively (Supplementary Figure S4A, Supplementary File S1: Results, figures, and tables). Also, technical variability was assessed by first computing the coefficient of variation (CV) per treatment condition (4 replicates), followed by averaging of all CV values across biological replicate and readout (Supplementary Table S8 Supplementary File S1: Results, figures, and tables). The technical variability was considered low for both readouts with a CV of 6%–9%. In addition, inter-experiment reproducibility was assessed by comparing data points (viability data only) common for the reference screen (Folkesson et al., 2020) and the high-throughput validation screen. The correlation coefficient for responses tested in both setups was 0.96 (Supplementary Figure S4B; Supplementary File S1: Results, figures, and tables). Response assessed by viability was strongly correlated (0.96) with the confluency response (Supplementary Figure S5, Supplementary File S1: Results, figures, and tables), which was in concordance with results from the previous screen (Folkesson et al., 2020). For the low-throughput screen, intra-experiment reproducibility was assessed for viability data by computing the Pearson correlation between all pairs of biological replicates. The correlation coefficients ranged from 0.84 to 0.96 (Supplementary Figure S4C, Supplementary File S1: Results, figures, and tables).

To adopt a regulatory network-based approach for the study of cellular effects of second and higher-order drug combinations we began with constructing a prior knowledge network encompassing all proteins and associated signaling pathways targeted by drugs in our previously published high-throughput combination screen (Folkesson et al., 2020) (Figures 1A, C). The resulting regulatory network comprised pathways mediating mitogenic signaling (MAPK, PI3K and TAK1/NF-kB pathways), pathways signaling the effect of DNA damage and repair induced by the chemotherapeutic agents oxaliplatin and 5FU (Bakkenist et al., 2018; Pilié et al., 2019), cell cycle components, and pathway-crosstalk. For more details on network construction, see Details on network construction, Supplementary File S1: Results, figures, and tables.

To enable study of network behavior upon perturbations, we formulated the regulatory network as a logical model using GINsim (Figure 2A), as described in Materials and Methods. In the model, two input nodes represented the effect of the chemotherapeutic agents 5FU and oxaliplatin, whereas all other drug targets were represented by ‘internal’ nodes. For assessment of model phenotype at baseline and upon perturbation, three phenotypic pre-output nodes (representing G1/S transition, growth arrest, apoptosis) were included in the model. These nodes summarized into the output node ‘Overall_phenotype’, capable of taking on values 1, 2, and 3, where 3 corresponds to a proliferating cell. While a fourth output level, Overall_phenotype = 0, was defined as an output level resulting from inactivity of all pre-output nodes, this state was not achievable, since either “G1/S transition” or “Growth arrest” was defined to always be active.

To study the inherent signaling features of the initial model, prior to directed adjustments, model stable states were computed for several different conditions (Figure 2B). Here, ‘Baseline’ corresponds to what we would define as an unperturbed cell in vitro (i.e., no drugs added, and no ongoing DNA damage). The stable state associated with the ‘Baseline’ condition displayed activity of growth-promoting pathways and anti-apoptotic proteins, whereas pathways involved in DNA damage response and repair were all inactive (Figure 2B). The output value associated with the baseline stable state was the highest possible, which is expected as the features mentioned have been described for viable and growing cancer cells in vitro (Evan and Vousden, 2001; Fattahi et al., 2020; Guo et al., 2020; Mukhopadhyay and Lee, 2020). Using the initial computational model, we studied the in silico effects of all single drugs tested in our high-throughput screen (Folkesson et al., 2020). This was done by computing stable states in the presence of in silico perturbations that either inactivated nodes (MEK, PI3K, CDK4_6_CycD, TAK1, PARP) or fixed them in the active state (Fluorouracil, Oxaliplatin) (Supplementary Table S6, Supplementary File S1: Results, figures, and tables). Each such perturbation was associated with a single stable state (Figure 2B). A common observation with all perturbations was that although we observed an effect on the activity of components (nodes) within the pathway that was primarily targeted by each specific drug, the activity of nodes associated with other pathways was rarely affected. Additionally, none of the single drugs showed an inhibiting effect on the global model output (i.e., Overall_phenotype = 3). These aspects of our model are necessary prerequisites to enable detecting responses to combinatorial perturbations that exceed the response of single perturbations. Overall, the results from the analysis of the initial model indicated that the model correctly captured biologically plausible signaling features present in unperturbed cancer cells in vitro. Furthermore, upon single node inhibition or induction of DNA damage, the model accurately engaged response and repair mechanisms reported in the literature, still while allowing additional effects of combination treatment to be observed. Altogether, these results indicate that the model sufficiently represents the cell fate decision mechanism to study synergistic drug combinations affecting DNA integrity as well as cellular growth.

As we seek to use the computational model for prediction of cellular response to higher-order drug combinations, i.e., combinations comprising more than two drugs, we first improved the model’s ability to predict the response to pairwise drug combinations. This was done by implementation of adjustments in model topology and logical rules as described in Materials and Methods. Drug response data for HCT-116 cells from our previous combination screen (Folkesson et al., 2020) were used as reference data for model adjustments, and the effects of model adjustments on the output level were closely monitored by continually comparing stable state analysis to reference data. The reason for focusing on HCT-116 cells when performing model adjustment was that the strongest as well as the largest number of synergies were observed for this cell line in the reference screen. Despite this, we wanted to investigate if it was possible to detect an increased effect by the addition of a third and possibly fourth agent. More details concerning specific model adjustments are summarized in Supplementary Table S4 (Supplementary File S1: Results, figures, and tables).

We found that the model’s predictive capacity improved by each stepwise adjustment(s) guided by predictive performance for HCT-116 cells (Figure 3). Of note, the inclusion of the node MLH1 as a knockout (KO) played an important role for the predictive power of the adjusted model. Here, the node knockout represents the HCT-116 cells’ homozygous loss of MLH1 (Liu et al., 1995). Following implementation of this and other (in total 13) adjustments (Figure 3), the model correctly captured 12 out of 13 synergies, and 7 out of 8 non-synergies observed in HCT-116 cells in vitro (Folkesson et al., 2020). Of note, due to the inherent nature of the model construction, where single-node perturbations were set to have to have no effect on the overall phenotype, none of the antagonistic pairwise drug combinations detected in our previous screen (many of them involving the chemotherapeutic agent oxaliplatin) were predicted as antagonistic by the model. Instead, all pairwise perturbations which did not show an effect on overall phenotype were classified as “non-synergies”. However, as our primary aim was to use the model for prediction of synergistic higher-order combinations and since the model was now configured to closely mimic the synergy pattern of pairwise combinations in HCT-116 cells (Figure 3), we hypothesized that the model would also be effective for prediction of synergistic effects of higher-order drug combinations in the same cell line (HCT-116).

Next, we used our adjusted computational model to study the effect size of pairwise combinations, to conclude whether synergistic higher-order combinations (>2 drugs) would at all be detectable by the model. To study the effect size of pairwise perturbations evaluated during model adjustment (Figure 3), associated stable states were computed for all combinations, followed by recording of the Overall_phenotype level and activities of model nodes. Maximum effect was defined as Overall_phenotype = 1, and hence an apoptotic phenotype, whereas Overall_phenotype = 2, characterized by cell cycle arrest in G1/S phase, was considered a medium effect. Perturbations with no effect on overall phenotype were associated with the highest possible value for Overall_phenotype, 3. Also, to gain further insights into possible mechanisms underlying observed synergistic effects of pairwise combinations, underlying single node perturbations were subjected to the same analysis (stable state analysis and recording of node activity/level), whereafter the differential activity of model nodes could be computed for all pairwise combinations as compared to the effect of underlying single node perturbations (Figure 4). We found that none of the 12 correctly predicted synergies entirely associated with the lowest output level that corresponds to an apoptotic phenotype, Overall_phenotype = 1 (Figure 4). Instead, we found that combined perturbations, when synergistic, often resulted in arrest of the cell cycle in G1/S phase (Figure 4), altogether resulting in the Overall_phenotype = 2. When studying model attractor reachability, one combination (PI + PD) was found capable of inducing apoptosis in a fraction of the simulations (Figure 4), but since the probability for such an outcome (assessed using the Assess attractor reachability algorithm in GINsim) was found to be less than 0.3, growth arrest (Overall_phenotype = 2) was selected as the main output also for this combination. The prediction of a cytostatic rather than cytotoxic synergistic effect of evaluated pairwise drug combinations was supported by data generated in the reference screen, which demonstrated an increase in confluency over time for all tested drug combinations (Supplementary Figure S6, Supplementary File S1: Results, figures, and tables). Here, a decrease in confluency relative to the starting point (drug addition) would have been sign of cytotoxicity. However, as no marker of apoptosis was included for most of the combinations evaluated in the reference screen, an additional apoptotic effect in vitro cannot be ruled out. Computationally, as none of the pairwise perturbations tested by the model did on their own induce a maximum effect (apoptosis/Overall_phenotype = 1), such an effect induced by e.g., perturbation of additional nodes would theoretically still be detectable. We therefore consider the model in its current adjusted state to be well suited for exploring combination effects of higher-order drug combinations in HCT-116 cells.

The model was next used for prediction of synergistic third and fourth-order combinations in HCT-116 cells. Combinations subjected to stable state analysis and all associated stable states are listed in Supplementary Table S9 and Supplementary Table S10 (Supplementary File S1: Results, figures, and tables). For third and fourth-order combinations, synergy was called when the effect of the highest-order combination was larger than the effect of the most effective lower-order treatment subset.

From 35 possible third-order drug combinations tested, three (PI + PD+5FU, PI + PD + OXA and PI + PD + OLA) were predicted as synergistic by the model (Figure 5A, Supplementary Table S9, Supplementary File S1: Results, figures, and tables). In all three cases, the probability of the stable state for the third-order drug combination was 1. Of note, all three synergistic third-order combinations involved combined inhibition of PI3K and MEK (PI + PD). Two of them also involved chemotherapeutic DNA damage. All three combinations were associated with Overall_phenotype = 1, and an apoptotic phenotype. According to the model, apoptosis was enabled by the activity of the pro-apoptotic factor BAX, which was synergistically activated upon third-order perturbation, but not by any of the simulated single or pairwise drug responses (Figure 5B). The model suggested that the presence of functional BAX was crucial for synergy of these third-order combinations. The latter was confirmed by performing simulations where the activity of BAX was turned off (knocked out), whereupon none of the three combinations showed synergy (results not shown). The suggested mechanism is supported by literature stating that higher-order drug combinations might be needed in order to effectively kill CRC cells and that direct targeting of processes such as apoptosis and cell cycle might be crucial for enhanced cytotoxicity (Horn et al., 2016). According to further stable state analysis of all possible fourth-order perturbations (35), none of these combinations contributed with any additional effect compared to third or second order combinations (Supplementary Table S10, Supplementary File S1: Results, figures, and tables).

Altogether, results from the tests indicate that the adjusted model can be used for prediction of especially third-order drug combinations with possible synergistic effects in HCT-116 cells.

To evaluate the predictive performance of the adjusted model, the effect of higher-order drug combinations in HCT-116 cells was assessed by a follow-up validation screen. An overview of the general screening procedure is shown in Figure 6A. During screening, cells were treated with all third-order drug combinations predicted to be synergistic by the model (PI + PD + OXA, PI + PD+5FU, OLA + PI + PD), two randomly selected third-order combinations with predicted non-synergy (OXA+5FU + PD, OXA+5FU + PI), one of the predicted fourth-order non-synergies (OXA+5FU + PI + PD), as well as all underlying lower-order treatments. Combination effects of second- and higher-order drug combinations were evaluated based on viability data by calculating the HSA excess per dose step (see illustration of the meaning of HSA excess in Figure 6A), followed by integration of HSA excess values across the entire dose range per combination, as described in Materials and Methods. We applied a cut-off at HSA = 0 for classification of synergy for higher-order drug combinations; hence, drug combinations with integrated HSA excess scores < 0 and > 0 were considered synergistic and antagonistic, respectively. While HSA as a reference model for synergy scoring often is considered less stringent in its classification of synergy compared to several other reference models (such as the Bliss independence and Loewe additivity models), and therefore should be used with extra care, its mathematical analogy to how we in this study defined synergy in silico made it suitable for quantification of synergy also in vitro. A complete list of all treatment conditions evaluated in the validation screen can be found in Supplementary Table SM2 and Supplementary Table SM3 (Supplementary File S2: Materials and Methods).

All three third-order combinations predicted as synergies by the computational model were confirmed statistically significant synergistic (p < 0.05) at multiple doses in HCT-116 cells in vitro in the high-throughput validation screen (Figure 6B; Supplementary Figure S7 and Supplementary Table S11, Supplementary File S1: Results, figures, and tables). We also found that the two combinations OXA+5FU + PI and OXA+5FU + PD, which were both predicted as non-synergies by the model, were non-synergistic in HCT-116 cells in vitro (Figure 6B; Supplementary Figure S8 and Supplementary Table S12, Supplementary File S1: Results, figures, and tables), indicating that the adjusted model has the power to correctly predict synergistic, as well as non-synergistic properties of multiple third-order combinations in HCT-116 cells.

To hypothesize on further model refinements related to modelling of DNA repair, we also screened three additional third-order drug combinations; 5FU + OLA + PD, 5FU + OLA + PI and 5FU + OLA+5Z. The computational model predicted none of these combinations to be synergistic, as 5FU combined with either of PD, PI and 5Z alone resulted in non-repairable DNA damage (Figure 4). Hence, the model did not indicate additional repair-compromising effects of PARP inhibition (OLA), although this could be expected from a biological point of view (see Selection of additional third-order drug combinations for screening, Supplementary File S1: Results, figures, and tables). These combinations were found to be synergistic when screening HCT-116 cells (Supplementary Figure S7), in line with biological reasoning, and indeed indicating that certain predictions most likely would benefit from a multileveled representation of non-repairable DNA damage.

Based on model predictions, none of the potential fourth-order combinations (35) were expected to act synergistically in HCT-116 cells (Supplementary Table S10, Supplementary File S1: Results, figures, and tables) and due to the experimental complexity of fourth-order combinations (14 effector components to be individually tested per fourth-order component), we chose to limit the experimental follow-up to only one such combination, OXA+5FU + PI + PD. The selected fourth-order combination displayed a negligible negative HSA excess value in HCT-116 cells (Figure 6C).

To summarize, our combined computational-experimental framework predicted and validated three novel third-order drug combinations demonstrating synergistic effects in HCT-116 cells.