Multiobjective optimization identifies cancer-selective combination therapies

Otto I. Pulkkinen, Prson Gautam, Ville Mustonen, Tero Aittokallio

The authors have declared that no competing interests exist.

https://doi.org/10.1371/journal.pcbi.1008538, Volume: 16, Issue: 12, Pages: 1-17

Article Type: research-article Article History

Publisher: Public Library of Science

- Facebook
- Twitter
- Linkedin
- Whatsapp
Altmetric

Table of Contents

Introduction
Results
Discussion
Materials and methods
Supporting information

Abstract

Combinatorial therapies are required to treat patients with advanced cancers that have become resistant to monotherapies through rewiring of redundant pathways. Due to a massive number of potential drug combinations, there is a need for systematic approaches to identify safe and effective combinations for each patient, using cost-effective methods. Here, we developed an exact multiobjective optimization method for identifying pairwise or higher-order combinations that show maximal cancer-selectivity. The prioritization of patient-specific combinations is based on Pareto-optimization in the search space spanned by the therapeutic and nonselective effects of combinations. We demonstrate the performance of the method in the context of BRAF-V600E melanoma treatment, where the optimal solutions predicted a number of co-inhibition partners for vemurafenib, a selective BRAF-V600E inhibitor, approved for advanced melanoma. We experimentally validated many of the predictions in BRAF-V600E melanoma cell line, and the results suggest that one can improve selective inhibition of BRAF-V600E melanoma cells by combinatorial targeting of MAPK/ERK and other compensatory pathways using pairwise and third-order drug combinations. Our mechanism-agnostic optimization method is widely applicable to various cancer types, and it takes as input only measurements of a subset of pairwise drug combinations, without requiring target information or genomic profiles. Such data-driven approaches may become useful for functional precision oncology applications that go beyond the cancer genetic dependency paradigm to optimize cancer-selective combinatorial treatments.

Cancer is diagnosed in nearly 40% of people in the U.S at some point during their lifetimes. Despite decades of research to lower cancer incidence and mortality, cancer remains a leading cause of deaths worldwide. Therefore, new targeted therapies are required to further reduce the death rates and toxic effects of treatments. Here we developed a mathematical optimization framework for finding cancer-selective treatments that optimally use drugs and their combinations. The method uses multiobjective optimization to identify drug combinations that simultaneously show maximal therapeutic potential and minimal non-selectivity, to avoid severe side effects. Our systematic search approach is applicable to various cancer types and it enables optimization of combinations involving both targeted therapies as well as standard chemotherapies.

Pulkkinen,Gautam,Mustonen,Aittokallio,and Finley: Multiobjective optimization identifies cancer-selective combination therapies

Introduction

Combination therapies have become the standard of care to treat many complex diseases, including advanced cancers, which have developed resistance to monotherapies. Due to the heterogeneity of the disease, there is an increasing need to identify targeted combinatorial therapies for each cancer patient individually based on patient-derived cell models, such as ex-vivo cell cultures or patient-derived organoids and xenograft models. Systematic approaches to finding new combinatorial treatments are often based on high-throughput screening of pre-clinical models that emphasize combination efficacy and/or synergy as the key determinants of the drug combination performance [1–4]. A combination is called synergistic if it has a higher than expected effect on cancer cell inhibition, where the expected combination effect is calculated based on monotherapy responses and a reference model, such as the highest single agent, Loewe additivity or Bliss independence [5].

The mechanistic basis of therapeutically beneficial combinations remain still poorly understood, but recent findings in cancer cell line screens suggest that targets of highly synergistic drug combinations are likely to interact at the level of signaling pathways, either targeting the same or functionally opposite pathways, in such a way that the combined effect is greater than the sum of the individual effects [6]. This can occur for example by combinatorial inhibition of two orthogonal pathways or by targeting compensatory mechanisms via pathway cross-talks [7]. However, the same combination mechanism may affect not only the cancerous cells, but manifest itself also as unintended adverse effects in healthy cells [8]. In this case, the combination is nonselective against the particular cancer cells, and may result in a degree of toxic effects as well. Therefore, combination synergy and expected efficacy alone are not sufficient determinants when designing safe and effective treatment options.

Identifying cancer-selective drug combinations amounts to maximizing combination efficacy in the cancer cells, while minimizing potential toxic effects to healthy cells. Such drug combination design can be considered as a multiobjective optimization problem (MOP), in which the two goals are pursued simultaneously [9, 10]. Even if multiobjective optimization has been used as a tool in planning radiotherapies [11, 12], to our knowledge, the study of Matlock et al. [13] is the only attempt to apply the mathematical optimization framework to designing drug combinations. They proposed numerical algorithms that converge to most effective and least toxic combinations of drug targets based on so-called target inhibition maps [14]. Even though such target-centric approach provides clues about the targets to inhibit or avoid in a given cell context, it cannot optimize the drugs to be used in the combination, and is not applicable to non-targeted therapies, such as standard chemotherapies or emerging immunotherapies.

In this work, we develop a novel drug-centric mathematical optimization framework for maximizing drug combination efficacy and minimizing its nonselectivity against a given cancer type. Nonselectivity is defined through the average effect of a drug compound across a large number of cells from various cancer types, and it is used as a surrogate of toxic effects. This definition also makes it possible to apply the framework using drug sensitivity measurements in the cancer cells alone, without requiring measurements on healthy cells, nor information on targets of the drugs or genomic profiles of the cancer cells. To exemplify the framework, we use pairwise combination measurements of 104 drug compounds in 60 cancer cell lines from the NCI-ALMANAC resource [15], and subsequently experimentally validate the most interesting predictions in BRAF-V600E melanoma cell line. We propose a new quantity of selective synergy, Bliss multiexcess, a multiobjective generalization of the traditional Bliss excess measure.

Results

Joint modeling of therapeutic and nonselective effects

The reduction in the rate of growth of a cell line by a drug compound is typically measured as a growth fraction Q, the relative number of live cells with respect to untreated control after a certain time (e.g., inhibition % after 72h of treatment). The therapeutic effect E of a drug i at concentration c_i, or a combination of a pair of drugs i, j at concentrations c_i, c_j, on a target cell line l is here defined as the negative logarithm of growth fraction Q,

Hence, the higher the therapeutic effect, the less viable the cells are after the treatment. One advantage of the logarithmic formulation is that the effect is additive, in that the predicted effect of two drugs acting independently, i.e. without interaction in the sense of Bliss independence [], is given by the sum E_i + E_j [, ]. We assume that higher order combination measurements are not available, so all modeling needs to rely on the measured effects of single drugs and pairs of drugs [, ]. To this end, we define the therapeutic effect of a combination

c = {(c_{i})}_{i \in D}

from a set

D

of N possible drugs by the following pair interaction model (also known as regression model [, ] due to its original use in regression of drug interaction parameters []):

where

E_{i j}^{X S} (c_{i}, c_{j}; l) = E_{i j} (c_{i}, c_{j}; l) - E_{i} (c_{i}; l) - E_{j} (c_{j}; l)

is the effect excess of the pair over Bliss independence model, or Bliss excess for short. The first sum in () yields the Bliss model effect E_B(c; l) of the combination c, and the sum over pairs of drugs is its Bliss excess

E_{B}^{X S} (c; l)

Let m(c) = |{i: c_i > 0}| denote the number of components in a drug combination c. We then have three different types of equations for m(c) = 1, m(c) = 2, and m(c) ≥ 2, as shown in the following subsections.

1. Monotherapy, m(c) = 1

For a monotherapy, $c_{i_{0}} = c$ for a drug i₀, and c_i = 0 otherwise. The therapeutic effect in Eq (3) equals the effect of the single compound i₀ at concentration c, i.e.

The right-hand side equals the experimental value of E obtained from the measured growth inhibition if data for drug i₀ at concentration c exists. Otherwise, it can be obtained from experimental results by fitting a Hill curve [] to the other, measured concentrations (see Fig 1 and ).

Fig 1

Selectivity analysis of the 104 individual drug compounds in NCI ALMANAC.

(A) Data consists of drug sensitivity measurements of the NCI-60 cell lines at various concentrations c_i, c_j of drug pairs (here vemurafenib and gefitinib, respectively). Maximal single compound activity, Q_min, of a drug is obtained by fitting Hill curves to the dose-response data. (B) Hierarchical clustering of cell lines and drug compounds according to their Q_min. Highly effective but nonselective drugs are clustered together and show up as dark horizontal bands in the heatmap. Leukemias (red color columns) cluster strongly together as they are more vulnerable to a large variety of treatments than the other cancer types. Melanomas (violet columns) form one major cluster, and the rest of the melanoma cell lines are associated with other cancer types. In particular, all but one melanoma cell line stand out from most of the other cell lines because of their susceptibility to vemurafenib (light blue box). The drug activity profile of vemurafenib is clustered together most closely with estramustine, an estradiol derivative with a nitrogen mustard moiety, and the anthracycline neoplastics epirubicin and idarubicin. (C) Multiobjective optimization of single drug compounds (i.e. M = 1) acting on Malme-3M cell line. Vemurafenib is clearly the optimal single-compound treatment because it is both effective and selective. Bortezomib and valrubicin, near the rightmost tip of the convex hull, are the least optimal drugs. Mithramycin leads to highest therapeutic effect but it is also highly nonselective, far more than the common cytotoxic treatment cisplatin. EGFR inhibitor gefitinib alone is not a Pareto optimal solution. The multi-coloured dots show that a substantially better efficacy can be achieved by combining drugs. For the combination of vemurafenib and gefitinib, a higher therapeutic effect than predicted by Bliss independence model emerges, yet the boost in efficacy comes with an increase in nonselective effect as well.

2. Two-drug combination, m(c) = 2

For a combination consisting of a pair of drugs, we have $c_{i_{1}} = c_{1}$ and $c_{i_{2}} = c_{2}$ for some drugs i₁ and i₂, and zero otherwise. Eq (3) then reads

Emergence of a Bliss excess term quantifying drug-drug interaction is a notable difference to the corresponding formula () for a monotherapy. However, also in this case the right-hand side exactly matches the experimental value obtained from measured growth inhibition if the combinatorial growth inhibition of drug pair (i₁, i₂) was measured at (c₁, c₂). We note, however, that typically much fewer data points exist for the interaction part than for monotherapies, and the modeling of the combination effect is therefore more sensitive to measurement errors.

3. Drug triplet or higher order combination, m(c) ≥ 3

No large-scale datasets covering majority of all combinations with at least three drug components currently exist. Therefore, model (3) of the combination effect is based on measured pairwise effects, and the predicted effect is necessarily an approximation,

For example, in the case m(c) = 3,

The number of terms increases rapidly with an increasing m(c) because an interaction term with all other drugs is introduced upon addition of a new drug component to a combination. Accordingly, the effect becomes more sensitive to measurement errors.

As the level of abstraction in the modeling of the therapeutic effect changes with the number of drug components m(c), the three cases will also be discussed separately later in this article.

Nonselective effect and multiobjective optimization

We model the nonselectivity of a drug compound by the mean effect over the set of all cell lines $L$ in the data,

For a large number of patient-derived samples from various tissue types, the average in comprises many somatically different cell lines. Hence, we can consider it as a proxy for the degree of nonselective toxic effects on an average (cancerous or non-cancerous) tissue. Indeed, it is known that widely cytotoxic (high mean effect over all the cell lines) chemotherapy drugs display a range of severe side effects due to their poor selectivity for cancerous tissue over normal tissue []. On the other hand, hormonal therapies and drugs targeting cellular signaling often have smaller mean effect, and their adverse effects are also typically more manageable. In the , we provide support to the assumption that the mean effect is a good proxy for nonselective toxic effects, and even clinical adverse effects, therefore serving as a gauge for combination safety.

As a drug combination should not benefit from additional drugs with hardly no effect, we add to each ${\bar{E}}_{i}$ a small regularizing parameter δ (we use the value 0.1) that counterbalances any arbitrarily small increase in therapeutic effect or decrease in the mean effect. We hence define the nonselective effect of a drug combination with m(c) components as

In the spirit of , the nonselective effect too can be split into nonselective effects

{\bar{E}}_{i} + δ

of isolated compounds, and mean excesses

{\bar{E}}_{i j}^{X S}

, which stem from the drug interactions.

As the nonselective effect acts as a surrogate for adverse effects and toxicity of a combination, it makes sense to ask, which drugs or their combinations have most therapeutic potential conditioning on a tolerance to nonselective effect. We divide the problem into parts of increasing complexity by restricting the number of drug components in combinations: We first study monotherapies only, then add in pairwise combinations, and finally allow for combinations of any order.

Let M denote the maximal number of drugs allowed in any drug combination c, in that m(c) ≤ M. To find optimal combinations, with at most M components, for a cancer cell line l that simultaneously maximize the therapeutic effect and minimize the nonselective effect is a multi-objective optimization problem (MOP) [9, 10]. Formulated using the so-called ϵ-constraint method [9, 26], the task is to find the set of concentration vectors

Here, Θ denotes the tolerance to nonselective effect, and a vector of drug concentrations that maximizes the therapeutic effect corresponds to each Θ. The set C of solutions of this continuous MOP traces, for each M, a Pareto front, a curve in the plane of nonselective and therapeutic effects, the points of which are called Pareto optimal. For each point on the Pareto front, no combination of drug concentrations with lower nonselective effect and higher therapeutic effect exists when M is fixed.

Dimensional reduction by the strongest effect of a combination

The fact that the concentrations c can be adjusted on a continuous scale means that the number of solutions at the Pareto front is infinite. The Pareto front is therefore a continuous curve in the plane of therapeutic and nonselective effects. Solving this continuum MOP is obviously a numerically demanding task, and a successful solution very much depends on availability of reliable combination measurements at a multitude of concentrations of all drug compounds, so that a faithful continuous approximation to the effect in the plane of concentrations can be found for each pair of drugs. To lower the risk of false predictions due to measurement error and/or insufficient data resolution, we inspect a discrete set of concentrations instead. The two-drug combination effects are typically quantified by n-by-n square matrices (e.g. n = 4, 8, see Fig 1) with n − 1 nonzero concentrations for both drugs in a combination. However, these concentrations are not necessarily the same for a single drug in each of those combination measurements, which complicates the analysis. Here we choose to reduce the problem further: We use binary variables N_i = 1 or 0, to determine whether a drug compound $i \in D$ is present in the combination or not, respectively. The effect of the whole drug combination is then predicted by the pair interaction model to be

where

is the strongest effect of drug i computed from a fitted Hill curve to growth fractions Q_i (see Fig 1 and ), and, with

(c_{i}^{*}, c_{j}^{*}) = {argmax}_{c_{i}, c_{j} > 0} E_{i j} (c_{i}, c_{j}; l)

is the Bliss excess at the maximal measured therapeutic effect of the combination i, j. This way of defining the parameters of the binary problem was found to constitute a robust model. In particular, as high-throughput drug screens often show variability within and between plates on which the cells are drugged, the strongest effect from a fitted Hill curve for each cell line yields a reliable single number for the therapeutic effect of a drug compound in monotherapy. Furthermore, the way the pairwise interaction constants

ε_{i j}^{X S}

are computed per matrix from measurements on the same plate, and not against the fitted Hill curves, respects plate-to-plate variance and significantly reduces the the risk of misestimating the drug interaction strengths.

In the binary model, the nonselective effect is again given by Eq (9) with the mean effect over all cell lines computed using the binary effect function. The notions of Bliss independence and Bliss excess of the whole combination carry over to the new definition in a natural way.

A consequence of dimensional reduction through Eqs (11–13) is that the concentration of a drug in combination with another drug gets fixed, and a different concentration may get chosen for the same drug in another pairwise combination. The combination triplet spanned by these two pairs would then paradoxically contain a single drug in two different concentrations, which may lead to overestimation of the combinatorial therapeutic effect. However, should this occur, the nonselective effect is likewise affected, and the MOP is stable against these changes.

We consider our solution to MOP a posteriori method [9], in which a human decision maker decides, which of the Pareto optimal points provided by the algorithm are clinically relevant and worth testing further in the laboratory. For example, the combinations taken into consideration should not be much more nonselective than broadly toxic standard chemotherapies. Even if this condition is controlled by the nonselective effect, the set of Pareto optimal or nearly optimal combinations may become large. Therefore ways of further delineating the benefits achieved by a particular combination in comparison to other Pareto optimal solutions are needed.

Therapeutic dominance

In multiobjective optimization, a vector is said to dominate another vector if it performs at least equally well as the other one in all features, and it outperforms the other vector at least in one feature [27]. Inspired by this concept, we define the therapeutic dominance of a Pareto optimal combination:

Let E(c*) and ${\bar{E}}_{δ} (c^{*})$ denote the therapeutic effect and the nonselective effect of a Pareto optimal combination c*, respectively, and let m(c*) denote its number of components. The therapeutic dominance of c* is defined as the difference of its therapeutic effect E(c*) and the highest therapeutic effect of all Pareto optimal combinations with at most m(c*) − 1 components and nonselective effect not exceeding ${\bar{E}}_{δ} (c^{*})$ . This directly quantifies the therapeutic advance gained by adding a drug to a Pareto optimal combination, or achieved by a monotherapy in comparison to no therapy at all.

Selective synergy

A drug combination can be Pareto optimal for several reasons. First, the individual therapeutic effects of the drug components may be strong enough to bring the combination to the Pareto front without any synergistic interaction. In this case, the Bliss excess $E_{B}^{X S} (c; l)$ of the combination is small. Second, the combination may exhibit a strong effect due to synergistic interaction of the drugs. However, synergism occurs in a large number of cell lines, leading to a comparable but not excessive increase in nonselective effect. This scenario is characterized by a higher-than-average Bliss excess, and a mean Bliss excess ${\bar{E}}_{B}^{X S} (c)$ of the same order. Third, the combination may show high selective synergy, in that the drug interaction is manifested as an increase in therapeutic effect, and only little, or even as a negative change, in the mean effect over cell lines. Selective synergy is evidenced by high Bliss excess and low mean Bliss excess.

The degree of selective synergy of a combination can be boiled down to a single number, Bliss multiexcess

where ϕ weighs the importance of adverse effects. A motivation for definition (), and a systematic way of choosing the value of ϕ comes from the scalarized (i.e. projected to one dimension), convex optimization problem corresponding to our MOP: Bliss multiexcess () equals the interaction part of its objective function []

which is to be maximized. Thus the combinations that are scalar-optimal for synergistic reasons have high Bliss multiexcess. This argument can be extended to combinations outside the convex hull by appropriately choosing the weighing parameter ϕ. As the Pareto optimal combinations near the convex hull can be seen as optional solutions to the scalarized problem, we use a single value of ϕ for all Pareto optimal combinations between two convex solutions—the smallest ϕ that yields the higher effect combination as an optimal solution of the scalarized problem.

Application to melanoma

We applied the optimization framework to designing new combination therapies for malignant melanoma. The optimization is based on the sensitivity data of 104 drug compounds and their pairwise combinations in the NCI-60 cell lines available in the NCI ALMANAC resource [28]. In this data, the drug-induced growth inhibition is tested in four-by-four matrices, as depicted in Fig 1A. Hierarchical clustering of drug compounds and the cell lines based on monotherapy sensitivities reveals high variation across the cell lines (Fig 1B). In particular, vemurafenib shows a strong effect in most of the melanoma cell lines, whereas in other cancer types the effect is less potent. This is expected, as vemurafenib is a selective inhibitor of V600E-mutant BRAF, and this mutation is present in approximately 60% of melanomas. Vemurafenib has become one of the main treatments for late-stage melanoma since its discovery in 2008 and consequent FDA approval in 2011 [29–31]. Despite the initial efficacy of vemurafenib in many patients, the prognoses still remain poor as many patients develop resistance to vemurafenib, often within 6 to 8 months from therapy initiation, due to reactivation of the MAPK/ERK and other pathways. This has led to interest in combinatorial inhibition of the MAPK/ERK pathway components.

Fig 1C shows the solution of the MOP for all 104 drug compounds in the MALME-3M melanoma cell line with the BRAF-V600E mutation. The Pareto front traces out the optimal combinations in the plane of nonselective versus therapeutic effects. Vemurafenib stands out as a highly selective and effective treatment, as expected. Vitamin A derivative tretinoin has a lower nonselective effect than vemurafenib, but its therapeutic effect is much smaller. On the other hand, gefitinib, an EGFR inhibitor upstream of BRAF in the MAPK/ERK signaling pathway, has slightly smaller therapeutic effect than vemurafenib but it is more nonselective. Notably, the combination of vemurafenib and gefitinib (Fig 1C, top) shows an equally strong therapeutic effect as the most potent monotherapy, mithramycin (plicamycin), and the combination is synergistic, in that the effect is higher than that predicted by the Bliss independence model.

The solutions of the MOPs for combinations of up to two drugs (Fig 2A, S1 Table) and higher-order combinations (Fig 2B, Table 1 and S2 Table) confirm that vemurafenib-gefitinib combination is Pareto optimal regardless of the number of components allowed in combinations. Other Pareto optimal combinations involve vemurafenib with anastrozole, tretinoin, thalidomide and sirolimus as partners. Of these, the anti-inflammatory drug thalidomide and the aromatase inhibitor anastrozole have little effect when used alone, indicating a true synergistic interaction with vemurafenib. The same compounds appear in the Pareto optimal combinations of three and higher order in Fig 2B. With a few exceptions, Pareto optimal combinations consist mostly of monotherapies and duplets at smaller therapeutic and nonselective effects than vemurafenib-gefitinib. A cloud of triplets that all include vemurafenib and gefitinib as a pairwise subcombination dominates at higher effects in Fig 2B.

Fig 2

(A) Multiobjective optimization (MOP) solutions for monotherapies and two-drug combinations (M = 2). Monotherapies are represented by big monochrome dots, and multicoloured dots indicate combinations of selected compounds (see color legend). The small light blue dots represent all the other drugs alone and in two-drug combinations. Anastrozole alone has next to no effect on any cell line, yet in combination with vemurafenib and gefitinib, a selective effect emerges. The triplet combination of these three drugs is shown for comparison, and it yields a significant advance over all pairs of the drugs. Pareto fronts for M = 1, 2 are shown. (B) MOP for drug triplets and their subcombinations (M = 3). Due to high number of possible combinations, only monotherapies of drugs indicated in the legend and drug combinations involving vemurafenib, gefitinib and anastrozole are shown. The triplet of vemurafenib and gefitinib with sirolimus is on the convex hull, indicating a strong selective effect, yet it comes with a significant cost in nonselective effect. The Pareto front for M = 4 coincides with the one for triplets within this (and actually much wider) window of nonselective effects, indicating that none of the fourth-order combinations perform better than drug triplets and their subcombinations.

Table 1

Selected Pareto optimal combinations for MALME-3M (See for the complete list).

Monotherapies

1. tretinoin

2. vemurafenib

Two-drug combinations:

3. anastrozole + vemurafenib

4. thalidomide + vemurafenib

5. tretinoin + vemurafenib

7. sirolimus + vemurafenib

9. gefitinib + vemurafenib

Drug triplets:

6.anastrozole + thalidomide + vemurafenib

8. anastrozole + sirolimus + vemurafenib

10 anastrozole + gefitinib + vemurafenib

11 gefitinib + tretinoin + vemurafenib,

12. gefitinib + vemurafenib + vismodegib

13. gefitinib + sirolimus + vemurafenib

The numbering of combinations is the same as in Figs 2B and 3.

Therapeutic dominance and the measure of selective synergy were used to further investigate the selected combinations (Fig 3). Vemurafenib monotherapy, along with vemurafenib-gefitinib and vemurafenib-gefitinib-sirolimus combinations are clearly the most dominant therapies, since they yield the most therapeutic advance in comparison to a lower order Pareto optimal combinations. However, such advance can be associated with a significant increase in nonselective effect in the case of combinations. For instance, the duplet of vemurafenib and gefitinib is synergistic, but not selectively synergistic. The vemurafenib duplets with anastrozole, thalidomide, tretinoin and sirolimus, and the triplets of vemurafenib, gefitinib and anastrozole, or sirolimus are truly selectively synergistic combinations with significant therapeutic dominance over their lower order combinations.

Fig 3

Therapeutic dominance (bars, left y-axis) and selective synergy measured by the Bliss multiexcess (red curve, right y-axis) of Pareto optimal combinations as a function of nonselective effect.

The numbered combinations correspond to Pareto optimal combinations highlighted in Fig 2B, from tretinoin (1) monotherapy to vemurarfenib+gefitinib+sirolimus (13) triplet. Vemurafenib monotherapy (2) is dominant at low mean effect. Vemurafenib + gefitinib (9) duplet has high a dominance but very low selective synergy, indicating that the combination is optimal because of its strong therapeutic effect, not due to selective drug interaction. The combination vemurafenib+gefitinib+sirolimus (13) has the highest selective synergy among shown combinations. The combinations vemurafenib+anastrozole (3) and vemurafenib+gefitinib+anastrozole (10) exhibit both considerable dominance and local maxima in selective synergy.

In-house validations

To further test the most interesting predictions from the MOP solution, we experimentally validated in-house two Pareto optimal drug triplets and their subcombinations in the BRAF-V600E melanoma cell line MALME-3M. This validation experiment was performed to test the accuracy of the predicted therapeutic effect in other experimental assay (FIMM) than that used for model predictions (NCI-60). The in-house drug combination assay was performed in 8x8 matrix format, similar to previous studies (25) (see Materials and methods for details and S1 Fig for an example).

Table 2 shows the maximal effect over the measured concentrations in both assays. In accord to predictions, the triplet combinations tested in-house (E_FIMM) show stronger effect than any of the duplets, and the observed therapeutic effects of triplets are very close to the predicted values (E_NCI). The vemurafenib monotherapy and its combinations with gefitinib had stronger effects in FIMM assay than those reported in NCI ALMANAC. This makes the Bliss excesses smaller in the in-house assays. However, as the concentration of the third component was fixed in the in-house assays, we lack the full three-dimensional optimization over concentrations, and hence the most synergistic vectors may have been missed.

Table 2

NCI ALMANAC and FIMM drug assays on MALME-3M.

Treatment	E_NCI		E_FIMM
vemurafenib	1.60	-	2.10(10)	-
gefitinib	1.13	-	1.39(6)	-
anastrozole	0.02	-	0.02(1)	-
vismodegib	0.04	-	0.03(1)	-
vem+gef	2.99	0.26	3.57	0.08
vem+ana	1.91	0.29	2.33	0.21
vem+vis	1.53	-0.11	2.41	0.28
gef+ana	1.57	0.42	1.73	0.32
gef+vis	2.13	0.96	2.15	0.73
vem+gef+ana(1 μM)	3.73	0.98	3.74	0.23
vem+gef+vis(1 μM)	3.89	1.12	3.73	0.21

E_NCI: Therapeutic effect predicted by the binary pair interaction model based on NCI ALMANAC data.

E_FIMM: Therapeutic effect measured at FIMM. The figure in parenthesis is the error of the last digit from fitting the Hill curve.

$E_{B}^{X S}$ : Bliss excesses computed either from the first or the third column.

As the pairwise matrices in our in-house experiments have higher resolution (8x8) than in the NCI ALMANAC assays (4x4), we were able to observe interesting nonlinearities in the growth inhibition patterns. For example, the combination of vemurafenib and gefitinib shows antagonistic pattern at high concentrations of vemurafenib and low concentrations of gefitinib, while it results in a synergistic effect at high concentrations of both components (see S2 Fig). Further, the experiments revealed a highly synergistic combination between gefitinib and vismodegib, a SMO inhibitor of the hedgehog pathway, and this combination scores high also in our MOP solution because of this synergy along the high additive effect and nonselective synergy of vemurafenib and gefitinib.

Discussion

We introduced a multiobjective optimization framework for designing effective and safe drug combinations by simultaneous optimization of therapeutic and nonselective effects. The MOP approach requires only the measurements of single drugs and a subset of their pairwise combinations, and it enables optimization of also higher-order combinatorial treatments of three or more drugs that are often required for treating advanced cancers. The nonselectivity was modeled here by the mean effect over all the cell lines in the NCI ALMANAC data set, yet any other control data can be used for the purpose of estimating toxic side effects. The exact MOP solutions from Pareto-optimization provide the decision maker with a targeted set of patient-specific combinations, illustrated and quantified by their therapeutic and nonselective effects. We further defined Bliss multiexcess, a multiobjective generalization of the traditional Bliss excess, as a new measure of selective synergy, and demonstrated how it enables one to solve two critical challenges of combinatorial therapeutics: (i) combinatorial explosion of the number of potential higher-order combinations to consider for a given patient, and (ii) optimization of the balance between therapeutic and nonselective effects to rule out broadly toxic combinations.

Our results confirmed the previously suggested combinatorial inhibition of the MAPK/ERK pathway as an effective and safe treatment for malignant melanoma with BRAF V600E mutation. However, we note that the approved combination of vemurafenib with cobimetinib [28, 29], or with any MEK inhibitor, was not tested in the NCI ALMANAC study, but our mechanistic-agnostic model pairs the EGFR inhibitor gefitinib with vemurafenib instead. Several higher order combinations were also found to be Pareto optimal. Of these, an aromatase inhibitor anastrozole is typically used to treat estrogen-responsive breast cancer, and it is therefore quite surprising to see it having selectivity for melanoma. On the other hand, already the unsupervised hierarchical clustering of drugs without any interaction information paired vemurafenib with estramustine, an estrogen receptor binder. There is increasing evidence that hormone levels play a role in the growth and progression of melanoma, although the underlying mechanisms are still unclear [30]. These examples demonstrate how our MOP solutions identify among a massive number of potential combinations those with optimal efficacy and safety profiles that have higher likelihood of success in further pre-clinical and clinical testing.

To select the Pareto optimal multi-drug combinations without relying on approximative modeling of combination efficacies, the current MOP framework requires as input the efficacies of a majority of combinations both in the target cell line and in a healthy control, or in a sufficiently large number of other cell lines, in order to compute the nonselective effects. More specifically, for finding the Pareto optimal drug combinations of at most M components from a set of N drug compounds, an order of N^M of measurements are required for each cell line used in the modeling. Such large-scale drug screens currently exist for monotherapies and two-drug combinations (M = 1, 2), and to apply the MOP framework to higher-order combinations, it must be augmented with a model that predicts the higher-order combination efficacies. In the present work, we chose the pair interaction model that uses only the measured efficacies of monotherapies and two-drug combinations for this purpose, but there are also other alternatives, which may be based on pair combination efficacies [21, 22] or use additional chemical and biological information such as quantitative structure-activity relationship analysis of the compounds, gene expression data, or signaling pathway and network annotations [8, 32].

The prediction of the combination effects was based on Hill equations fitted to monotherapy dose-response data, and pairwise Bliss excesses, to increase the model robustness against measurement errors. Since the maximal effect Q_min of a combination may be sensitive to measurement errors, alternative measures of combination efficacy, such as the mean effect over the concentration pairs, could be used instead. However, since the most promising Pareto optimal combinations require anyway further experimental validation over a finer grid of concentrations, any potential erroneous solutions based on outlier measurements can be then identified in this validation phase. We further note that, unless the concentration vectors of single drugs were fixed in all the dose-response measurements, even averaging of the pairwise effects over concentrations will not completely resolve the apparent paradox, where a drug may be present at two different concentrations in one combination. Full multiobjective optimization in the space of continuous concentrations will be needed to remove the superfluous degree of freedom, and developing for it a theoretical and numerical methodology that meets the challenges posed by both low data resolution and possible measurement errors requires further studies.

The multiobjective optimization framework presented in this article was designed with high-throughput experimentation in mind. Datasets such as the NCI ALMANAC, consisting of nearly 3 million data points, allow this kind of analysis. Even if datasets of this size are still rare, and can be produced only with larger resources, the multiobjective optimization method is also applicable to smaller-scale combinatorial datasets that can be produced in many laboratories with the help of robotic instruments, in particular if combined with pre-existing data from other studies. There area also public resources, such as PharmacoDB [33] and SYNERGxDB [34], that integrate monotherapy and drug combination data from published studies in several cancer types, and can therefore be used in the MOP approach, either as target cell lines or as background data for toxicity modelling. Dose-reponses in normal cells can be also used to asses the nonselective toxicity of new drugs to be incorporated in the model and in the optimization, in those applications where such control cell lines are available.

The caveat of all the models that are based on cancer cell lines is that they are monoclonal and miss many factors important for in vivo drug responses, including tumor microenvironment. This is in contrast to real tumors, which comprise heterogeneous collections of cells with subclonal structure. Therefore, also the response to a therapy in vivo is a mixture of responses from the subclonal populations, bringing the total therapeutic effect closer to the nonselective effect. The tumors grow within special microenvironments that permit carcinogenesis and allow tumor maintenance. As a consequence, most effective combination therapies, in vivo, are probably multi-targeted and attack the cancer from multiple sides by slowing down or stopping tumor cell proliferation, by thwarting drug resistance, and by blocking the required support coming from the tumor microenvironment [35]. Nonselectivity may thus be advantageous in some cases. However, the complementary objectives of drug efficacy and safety are still valid. The multiobjective method can also be extended to drug profiling data of primary patient samples (e.g. blood cancers) or complex 3D models like organoids, co-cultures or even tissue slices, which contain clonal heterogeneity and tumor microenvironmental structure. In particular, when considering heterogenous cell population within tumors and tumor microenvironment one would need to design targeted drug combinations that co-inhibit only the malignant cells and avoid severe inhibition of healthy cells. Extending the MOP formulation to sub-clonal drug response data is one important future aim, once such data will become available in larger-scale.

In summary, we have introduced a novel conceptual approach for optimizing patient-specific combinatorial treatments, and demonstrated its potential both in published large-scale drug combination data as well as using in-house validation experiments. The approach is fast, does not require any genomic profiling of the patient cells, is applicable to various cancer types, and enables optimization of combinations involving both targeted therapies as well as standard chemotherapies. We hope that the systematic optimization framework will prove useful in designing new effective and safe cancer therapies to further reduce the death rates and the toxic effects of treatments.

Materials and methods

Cell line for in-house experiments

Melanoma cell line Malme-3M (ATCC HTB-64) was purchased from ATCC. It was maintained at 37°C with 5% CO2 in a humidified incubator in IMDM (Cat#12440046) supplemented with 20% fetal bovine serum (Cat#10270-106) and 1% Penicillin-Streptomycin (Cat#15140-122). All the reagents were purchased from ThermoFisher Scientific for the in-house assays.

Drug combination assays

The drug combination screening approach described previously [36] was adopted on MALME-3M cell line. Seven different concentrations in log3-fold dilution ranging from 10 nM to 10000 nM of vemurafenib and gefitinib, vemurafenib and anastrozole, and gefitinib and anastrozole were combined with each other in 8x8 matrix formats in duplicate. For the triple combination, 500 nM of vemurafenib or 500 nM gefitinib or 1000 nM of anastrozole was added to the alternate drug combination matrix (in one copy). Exactly as described above, another set of combination was prepared for vemurafenib, gefitinib and vismodegib. For triple combination, 1000 nM of vismodegib was used whereas the concentrations of vemurafenib and gefitinib were the same as mentioned above (Supplementary file drug_combo_design.xlsx). The compounds were plated to black clear bottom 384-well plates (Corning #3764) using an Echo 550 Liquid Handler (Labcyte). 100 μM benzethonium chloride (BzCl2) and 0.1% dimethyl sulfoxide (DMSO) were used as positive and negative controls respectively. All subsequent liquid handling was performed using MultiFlo FX multi-mode dispenser (BioTek). The pre-dispensed compounds were dissolved in 5 μl of culture media containing cytotoxicity measurement reagents CellTox Green (Promega) and left in plate shaker at RT for 30 min. Twenty microliter cell suspension (75000cells/ ml) was dispensed in the drugged plates. After 72 h incubation, cytotoxicity (fluorescence, 485/520 nm excitation/emission filters) was recorded using PheraStar plate reader (BMG Labtech). Subsequently, 25 μl per well of CellTiter-Glo (Promega) reagent was added, and after 10 min of incubation at room temperature, luminescence (cell viability) was measured using PheraStar plate reader (BMG Labtech).

Acknowledgements

The authors thank Kaisa Miettinen (University of Jyväskylä) for comments concerning multiobjective optimization, Alexander Ianevski and Anil Kumar (FIMM) for sharing their Excel spreadsheet of the NCI ALMANAC data, and Laura Turunen (FIMM HTB unit) for the drugs plates in the in-house experiments.

References

LHe, EKulesskiy, JSaarela, LTurunen, KWennerberg, TAittokallio, et al Methods for High-throughput Drug Combination Screening and Synergy Scoring. Methods in molecular biology (Clifton, NJ). 2018;1711:351–398. 10.1007/978-1-4939-7493-1_17

AMalyutina, MMMajumder, WWang, APessia, CAHeckman, JTang. Drug combination sensitivity scoring facilitates the discovery of synergistic and efficacious drug combinations in cancer. PLOS Computational Biology. 2019;15(5):1–19. 10.1371/journal.pcbi.1006752

BAl-Lazikani, UBanerji, PWorkman. Combinatorial drug therapy for cancer in the post-genomic era. Nature Biotechnology. 2012;30(7):679–692. 10.1038/nbt.2284

AHCVlot, NAniceto, MPMenden, GUlrich-Merzenich, ABender. Applying synergy metrics to combination screening data: agreements, disagreements and pitfalls. Drug Discovery Today. 2019;24(12):2286–2298. 10.1016/j.drudis.2019.09.002

JTang, KWennerberg, TAittokallio. What is synergy? The Saariselkä agreement revisited. Frontiers in Pharmacology. 2015;6:181 10.3389/fphar.2015.00181

MPMenden, DWang, MJMason, BSzalai, KCBulusu, YGuan, et al Community assessment to advance computational prediction of cancer drug combinations in a pharmacogenomic screen. Nature Communications. 2019;10(1):2674 10.1038/s41467-019-09799-2

SJaeger, AIgea, RArroyo, VAlcalde, BCanovas, MOrozco, et al Quantification of Pathway Cross-talk Reveals Novel Synergistic Drug Combinations for Breast Cancer. Cancer Research. 2017;77(2):459–469. 10.1158/0008-5472.CAN-16-0097

KCBulusu, RGuha, DJMason, RPILewis, EMuratov, YKMotamedi, et al Modelling of compound combination effects and applications to efficacy and toxicity: state-of-the-art, challenges and perspectives. Drug Discovery Today. 2016;21(2):225–238. 10.1016/j.drudis.2015.09.003

KMiettinen. Nonlinear Multiobjective Optimization. vol. 12 of International Series in Operations Research & Management Science. Boston, USA: Kluwer Academic Publishers; 1999.

JBranke, KDeb, KMiettinen, RSłowiński. Multiobjective Optimization: Interactive and Evolutionary Approaches. Berlin, Heidelberg: Springer-Verlag; 2008.

HWHamacher, KHKüfer. Inverse radiation therapy planning—a multiple objective optimization approach. Discrete Applied Mathematics. 2002;118(1):145–161. 10.1016/S0166-218X(01)00261-X

DLCraft, TSHong, HAShih, TRBortfeld. Improved Planning Time and Plan Quality Through Multicriteria Optimization for Intensity-Modulated Radiotherapy. International Journal of Radiation Oncology*Biology*Physics. 2012;82(1):e83–e90. 10.1016/j.ijrobp.2010.12.007

KMatlock, NBerlow, CKeller, RPal. Combination therapy design for maximizing sensitivity and minimizing toxicity. BMC Bioinformatics. 2017;18(Suppl 4):116 10.1186/s12859-017-1523-1

NBerlow, LEDavis, ELCantor, BSéguin, CKeller, RPal. A new approach for prediction of tumor sensitivity to targeted drugs based on functional data. BMC Bioinformatics. 2013;14(1):239 10.1186/1471-2105-14-239

SLHolbeck, RCamalier, JACrowell, JPGovindharajulu, MHollingshead, LWAnderson, et al The National Cancer Institute ALMANAC: A Comprehensive Screening Resource for the Detection of Anticancer Drug Pairs with Enhanced Therapeutic Activity. Cancer Res. 2017;77(13):3564–3576. 10.1158/0008-5472.CAN-17-0489

CIBliss. The toxicity of poisons applied jointly. Annals of Applied Biology. 1939;26(3):585–615. 10.1111/j.1744-7348.1939.tb06990.x

DYBaeder, HozéYu G N, JRolff, RRRegoes. Antimicrobial combinations: Bliss independence and Loewe additivity derived from mechanistic multi-hit models. Philosophical Transactions of the Royal Society B: Biological Sciences. 2016;371(1695):20150294 10.1098/rstb.2015.0294

ATendler, AZimmer, AMayo, UAlon. Noise-precision tradeoff in predicting combinations of mutations and drugs. PLOS Computational Biology. 2019;15(5):1–17. 10.1371/journal.pcbi.1006956

KWood, SNishida, EDSontag, PCluzel. Mechanism-independent method for predicting response to multidrug combinations in bacteria. Proceedings of the National Academy of Sciences. 2012;109(30):12254–12259. 10.1073/pnas.1201281109

KBWood. Pairwise interactions and the battle against combinatorics in multidrug therapies. Proceedings of the National Academy of Sciences. 2016;113(37):10231–10233. 10.1073/pnas.1612365113

AZimmer, IKatzir, EDekel, AEMayo, UAlon. Prediction of multidimensional drug dose responses based on measurements of drug pairs. Proceedings of the National Academy of Sciences. 2016;113(37):10442–10447. 10.1073/pnas.1606301113

AZimmer, ATendler, IKatzir, AMayo, UAlon. Prediction of drug cocktail effects when the number of measurements is limited. PLOS Biology. 2017;15(10):1–16. 10.1371/journal.pbio.2002518

AWeiss, RHBerndsen, XDing, CMHo, PJDyson, Hvan den Bergh, et al A streamlined search technology for identification of synergistic drug combinations. Scientific Reports. 2015;5(1):14508 10.1038/srep14508

AVHill. The possible effects of the aggregation of the molecules of hæmoglobin on its dissociation curves. The Journal of Physiology. 1910;40:i–vii.

ROun, YEMoussa, NJWheate. The side effects of platinum-based chemotherapy drugs: a review for chemists. Dalton Trans. 2018;47:6645–6653. 10.1039/C8DT00838H

YYHaimes, LSLasdon, DAWismer. On a Bicriterion Formulation of the Problems of Integrated System Identification and System Optimization. IEEE Transactions on Systems, Man, and Cybernetics. 1971;SMC-1(3):296–297. 10.1109/TSMC.1971.4308298

MTMEmmerich, AHDeutz. A tutorial on multiobjective optimization: fundamentals and evolutionary methods. Natural Computing. 2018;17(3):585–609. 10.1007/s11047-018-9685-y

AVan Hoeck, NHTjoonk, Rvan Boxtel, ECuppen. Portrait of a cancer: mutational signature analyses for cancer diagnostics. BMC Cancer. 2019;19(1):457 10.1186/s12885-019-5677-2

JTsai, JTLee, WWang, JZhang, HCho, SMamo, et al Discovery of a selective inhibitor of oncogenic B-Raf kinase with potent antimelanoma activity. Proceedings of the National Academy of Sciences. 2008;105(8):3041–3046. 10.1073/pnas.0711741105

RFisher, JLarkin. Vemurafenib: a new treatment for BRAF-V600 mutated advanced melanoma. Cancer Manag Res. 2012;4:243–252. 10.2147/CMAR.S25284

AKim, MSCohen. The discovery of vemurafenib for the treatment of BRAF-mutated metastatic melanoma. Expert Opin Drug Discov. 2016;11(9):907–916. 10.1080/17460441.2016.1201057

FCheng, IAKovács, ALBarabási. Network-based prediction of drug combinations. Nature Communications. 2019;10(1):1197 10.1038/s41467-019-09186-x

PSmirnov, VKofia, AMaru, MFreeman, CHo, NEl-Hachem, et al PharmacoDB: an integrative database for mining in vitro anticancer drug screening studies. Nucleic Acids Research. 2017;46(D1):D994–D1002. 10.1093/nar/gkx911

HSeo, DTkachuk, CHo, AMammoliti, ARezaie, SAMadani Tonekaboni, et al SYNERGxDB: an integrative pharmacogenomic portal to identify synergistic drug combinations for precision oncology. Nucleic Acids Research. 2020;48(W1):W494–W501. 10.1093/nar/gkaa421

APetrelli, SGiordano. From Single- to Multi-Target Drugs in Cancer Therapy: When Aspecificity Becomes an Advantage. Current medicinal chemistry. 2008;15:422–32. 10.2174/092986708783503212

PGautam, LKarhinen, ASzwajda, SKJha, BYadav, TAittokallio, et al Identification of selective cytotoxic and synthetic lethal drug responses in triple negative breast cancer cells. Molecular Cancer. 2016;15(1):34 10.1186/s12943-016-0517-3