4 Synergy and sensitivity analysis
4.0 Baseline Correction
Typically, the models for estimating synergy or sensitivity of drug combinations expect that the drug’s effect to be expressed as continuous values ranged between 0 and 1 (or between 0% and 100%). However, in practice the readouts from the drug combination screening could be negative. In this cases, the technical error might be introduced in the dose-response data. The base line correction function is designed to adjust the “baseline” of the dose-response data closer to 0. Following is the main steps for the process of “baseline correction”:
- Extract all the single drug treatment (monotherapy) dose-response data from the combination.
- Fit the four-parameter log logistic model for each monotherapy with the data table from step 1.
- Pick the minimum fitted response value from all the models fitted in step 2 as “baseline”.
- Adjusted response value = observed response value - ((100 – observed response value) / 100 * baseline).
By setting the value of the parameter correct_baseline while calling function CalculateSynergy or CalculateSensitivity, user could chose:
- “non” - do not correct base line;
- “part” - correct base line but only adjust negative response values in matrix;
- “all” - correct base line with adjusting all values in matrix.
4.1 Drug synergy scoring
The synergistic effect can be determined as the excess of observed effect over expected effect calculated by a reference models (synergy scoring models). All of the models make different assumptions regarding the expected effect. Currently, 4 reference models are available in SynergyFinder.
- Highest Single Agent (HSA) [2] states that the expected combination effect equals to the higher effect of individual drugs:
\[y_{HSA} = max(y_1, y_2)\] where \(y_1\) and \(y_2\) are the monotherapy effect of combined drug 1 and drug 2
- Bliss model [3] assumes a stochastic process in which two drugs exert their effects independently, and the expected combination effect can be calculated based on the probability of independent events as:
\[ y_{Bliss} = y_1 + y_2 - y_1 \cdot y_2 \] where \(y_1\) and \(y_2\) are the monotherapy effect of combined drug 1 and drug 2
- Loewe additivity model [4] defines the expected effect \(y_{Loewe}\) as if a drug was combined with itself, i.e. \(y_{Loewe} = y_1(x_1 + x_2) = y_2(x_1 + x_2)\). Unlike the HSA and the Bliss independence models, which give a point estimate using different assumptions, the Loewe additivity model considers the dose-response curves of individual drugs. The expected effect \(y_{Loewe}\) must satisfy:
\[ \frac {x_1}{\chi_{Loewe}^1} + \frac{x_2}{\chi_{Loewe}^2} = 1 \]
where \(x_{1,2}\) are drug doses and \(\chi_{Loewe}^1,\ \chi_{Loewe}^2\) are the doses of drug 1 and 2 alone that produce \(y_{Loewe}\). Using 4-parameter log-logistic (4PL) curves to describe dose-response curves the following parametric form of previous equation is derived:
\[ \frac {x_1}{m_1(\frac{y_{Loewe}-E_{min}^1}{E_{max}^1 - y_{Loewe}})^{\frac{1}{\lambda_1}}} + \frac{x_2}{m_2(\frac{y_{Loewe}-E_{min}^2}{E_{max}^2 - y_{Loewe}})^{\frac{1}{\lambda_2}}} = 1 \] where \(E_{min}, E_{max}\in[0,1]\) are minimal and maximal effects of the drug, \(m_{1,2}\) are the doses of drugs that produce the midpoint effect of \(E_{min} + E_{max}\), also known as relative \(EC_{50}\) or \(IC_{50}\), and \(\lambda_{1,2}(\lambda>0)\) are the shape parameters indicating the sigmoidicity or slope of dose-response curves. A numerical nonlinear solver can be then used to determine \(y_{Loewe}\) for (\(x_1\), \(x_2\)).
- Zero Interaction Potency (ZIP) [5] calculates the expected effect of two drugs under the assumption that they do not potentiate each other, i.e. both the assumptions of the Loewe model and the Bliss model are met:
\[ y_{ZIP} = \frac{(\frac{x_1}{m_1}) ^ {\lambda_1}}{(1 + \frac{x_1}{m_1})^{\lambda_1}}+\frac{(\frac{x_2}{m_2})^{\lambda_2}}{(1 + \frac{x_2}{m_2})^{\lambda_2}} - \frac{(\frac{x_1}{m_1})^{\lambda_1}}{(1 + \frac{x_1}{m_1})^{\lambda_1}} \cdot \frac{(\frac{x_2}{m_2})^{\lambda_2}}{(1 + \frac{x_2}{m_2})^{\lambda_2}} \] where \(x_{1, 2}, m_{1, 2}, and\ \lambda_{1, 2}\) are defined in the Loewe model.
For the input data without replicates, the statistical analysis for synergy scores of the whole dose-response matrix is estimated under the null hypothesis of non-interaction (zero synergy score). Suppose that \(S_1, S_2, ...S_n\) are the synergy scores for all the concentration combinations (excluding monotherapies) for the dose-response matrix with mean of \(\bar s\). The p-value is:
\[p = exp(-0.717z - 0.416z^2),\ where\ z = abs(\bar s')/\sqrt{\frac{1}{n - 1}\sum_{i = 1}^n(s_i' - \bar s')^2}\]
The function CalculateSynergy is designed to calculate the synergy scores. Important parameters are:
data: The R object generated by functionReshapeData;method: The reference model for synergy score calculation. ParameterEminandEmaxare used to set the corresponding parameters in the 4 parameter log-logistic model which will be fitted to calculate “ZIP” and “Loewe” synergy scores;correct_baseline: Whether to adjust base line of response matrix.
res <- CalculateSynergy(
data = res,
method = c("ZIP", "HSA", "Bliss", "Loewe"),
Emin = NA,
Emax = NA,
correct_baseline = "non")The output adds one data frame “synergy_scores” to the input R object. It contains:
- “block_id” - The identifier for drug combination blocks;
- “concX” - The concentration for combined drugs;
- “ZIP_ref”, “HSA_ref”, “Bliss_ref”, and “Loewe_ref” - The reference additive effects calculated from corresponding models;
- “ZIP_synergy”, “HSA_synergy”, “Bliss_synergy” and "Loewe_synergy - The synergy scores calculated from corresponding models;
- “ZIP_fit” - The fitted %inhibition for combinations which is used to calculate the ZIP synergy score.
The mean of synergy scores for the combination matrix (excluding the monotherapy observations) and the p-values for them are added to the “drug_pairs” table in the input R object.
res$drug_pairs
#> # A tibble: 2 × 15
#> block_id drug1 drug2 conc_…¹ conc_…² input…³ repli…⁴ ZIP_s…⁵ HSA_s…⁶ Bliss…⁷
#> <int> <chr> <chr> <chr> <chr> <chr> <lgl> <chr> <chr> <chr>
#> 1 1 ispine… ibru… nM nM viabil… FALSE 3.51e-… 5.31e-… 1.14e-…
#> 2 2 canert… ibru… nM nM viabil… FALSE 9.21e-… 2.72e-… 6.99e-…
#> # … with 5 more variables: Loewe_synergy_p_value <chr>, ZIP_synergy <dbl>,
#> # HSA_synergy <dbl>, Bliss_synergy <dbl>, Loewe_synergy <dbl>, and
#> # abbreviated variable names ¹conc_unit1, ²conc_unit2, ³input_type,
#> # ⁴replicate, ⁵ZIP_synergy_p_value, ⁶HSA_synergy_p_value,
#> # ⁷Bliss_synergy_p_value
str(res$synergy_scores)
#> tibble [72 × 12] (S3: tbl_df/tbl/data.frame)
#> $ block_id : int [1:72] 1 1 1 1 1 1 1 1 1 1 ...
#> $ conc1 : num [1:72] 2500 2500 2500 2500 2500 2500 625 625 625 625 ...
#> $ conc2 : num [1:72] 50 12.5 3.125 0.781 0.195 ...
#> $ ZIP_fit : num [1:72] 93.3 92.3 82.7 75.2 57.3 ...
#> $ ZIP_ref : num [1:72] 79.4 79.4 67.9 37.2 37.2 ...
#> $ ZIP_synergy : num [1:72] 13.9 12.9 14.8 38.1 20.1 ...
#> $ HSA_ref : num [1:72] 71.3 54.2 54.2 54.2 54.2 ...
#> $ HSA_synergy : num [1:72] 20.89 38.96 30.7 21.29 7.75 ...
#> $ Bliss_ref : num [1:72] 86.9 78.6 73.1 45.9 52.5 ...
#> $ Bliss_synergy: num [1:72] 5.34 14.54 11.85 29.6 9.47 ...
#> $ Loewe_ref : num [1:72] 62.3 62.3 62.3 62.3 -14.9 ...
#> $ Loewe_synergy: num [1:72] 29.9 30.8 22.6 13.2 76.9 ...4.2 Sensitivity scoring
SynergyFinder mainly calculates 3 sensitive scores: relative IC50 and relative inhibition (RI) for single drug treatment, and combination sensitivity score (CSS) for drug combinations.
- Relative IC50: SynergyFinder extracts the single drug treatment response data from the combination matrix to fit the 4-parameter log-logistic curve:
\[y = y_{min} + \frac{y_{max} - y_{min}}{1 + 10^{\lambda(log_{10}IC_{50} - x')}}\] where \(y_{min}, y_{max}\) are minimal and maximal inhibition and \(x' = log_{10}x\)
- RI: RI is the abbreviation of “relative inhibition”. It is the proportion of the area under the dose-response curve (fitted by four parameter log logistic model), to the maximal area that one drug can achieve at the same dose range (Figure 1). For example, RI of 40 suggests that the drug is able to achieve 40% of the maximal of inhibition (e.g. a positive control where the cell is 100% inhibited in each dose that in the range). RI is comparable between different scenarios, even for the same drug-cell pair but tested in different concentration ranges.
Figure 1. Concept of RI
- CSS: CSS is the abbreviation of “Combination Sensitivity Score”. It is derived using relative IC50 values of compounds and the area under their dose-response curves. CSS of a drug combination is calculated such that each of the compounds (background drug) is used at a fixed concentration (its relative IC50) and another is at varying concentrations (foreground drug) resulting in two CSS values, which are then averaged. Each drug’s dose-response is modeled using 4-parameter log-logistic curve. The area under the log-scaled dose-response curve (AUC) is then determined according to:
\[AUC = \int_{c_1}^{c_2}y_{min} + \frac{y_{max}-y_{min}}{1 + 10^{\lambda(log_{10}IC_{50}-x')}}dx'\] where \([c_1, c_2]\) is the concentration range of the foreground drug tested.[6]
The function CalculateSensitivity is designed to calculate the synergy scores. Important parameters are:
data: The R object generated by functionReshapeData;correct_baseline: Whether to adjust base line of response matrix
Following sensitivity scores are added to the “drug_pairs” table in the input R object:
- “ic50_x”: The relative IC50 for drug x;
- “ri_x”: The RI for drug x;
- “cssx_ic50y”: The css for drug x while drug y is at its relative IC50;
- “css”: The CSS for the combination.
sensitive_columns <- c(
"block_id", "drug1", "drug2",
"ic50_1", "ic50_2",
"ri_1", "ri_2",
"css1_ic502", "css2_ic501", "css")
res$drug_pairs[, sensitive_columns]
#> # A tibble: 2 × 10
#> block_id drug1 drug2 ic50_1 ic50_2 ri_1 ri_2 css1_i…¹ css2_…² css
#> <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 ispinesib ibrutinib 2500 2.74 60.0 27.0 85.9 82.8 84.3
#> 2 2 canertinib ibrutinib 973. 1.44 -48.0 45.2 -35.0 -6.99 -21.0
#> # … with abbreviated variable names ¹css1_ic502, ²css2_ic501
