Methodology

Occupancy

Define the occupancy level occup, i.e. the fractional (enzymatic) effect or observed effect relative to maximal effect, for both compounds at given dose levels as $$ \textrm{occup}_{1}\left(d_{1}\right) = \frac{1}{1 + \left(\frac{\operatorname{EC50}_{1}}{d_{1}}\right)^{h_{1}}} $$ $$ \textrm{occup}_{2}\left(d_{2}\right) = \frac{1}{1 + \left(\frac{\operatorname{EC50}_{2}}{d_{2}}\right)^{h_{2}}} $$ Alternatively, the above equations can be rearranged to express dose in terms of occupancy so that $$ d_{1} = \operatorname{EC50}_{1} \left(\frac{1}{\operatorname{occup_{1}}} - 1 \right)^{-1/h_{1}} $$ $$ d_{2} = \operatorname{EC50}_{2} \left(\frac{1}{\operatorname{occup_{2}}} - 1 \right)^{-1/h_{2}} $$ Although the occupancy was considered here in the marginal case, it is equally well-defined when compounds are combined and is understood as the fraction of enzyme bound to any compound. It can thus be used to re-express classical Loewe additivity equations.

Classical Loewe model

In the classical Loewe model where both marginal models share upper (m) and lower (b) asymptotes, occupancy is defined as the solution to this additivity equation for each dose combination (d₁, d₂), namely $$\frac{d_1\left(\textrm{occup}^{-1} - 1\right)^{1/h_{1}}}{\textrm{EC50}_{1}}+ \frac{d_2\left(\textrm{occup}^{-1} - 1\right)^{1/h_{2}}}{\textrm{EC50}_{2}} = 1$$

Once occupancy is computed, in the classical Loewe model the predicted response at dose combination (d₁, d₂) can be calculated to be

Generalized Loewe model

Generalized Loewe model extends the classical Loewe model by allowing compounds to have different upper asymptotes so that when adjusted, the above predicted response is written instead as $$ y = b + \textrm{occup} \times \left[\frac{\left(m_{1} - b\right) d_{1}\left(\textrm{occup}^{-1} - 1\right)^{1/h_{1}}}{\textrm{EC50}_{1}} + \frac{\left(m_{2} - b\right) d_{2}\left(\textrm{occup}^{-1} - 1\right)^{1/h_{2}}}{\textrm{EC50}_{2}}\right]$$

In particular, if m₁ = m₂, then generalized Loewe is equivalent to the classical Loewe.

Highest Single Agent

A null model based on the Highest Single Agent (HSA) model does not attempt to model interaction effects at all and the predicted effect of a combination is either the minimum (if marginal curves are decreasing) or the maximum (if marginal curves are increasing) of both monotherapy curves.

Bliss independence model

Bliss independence implies that two agents do not cooperate, i.e. act independently of each other. Additionally, the assumption is that a decreasing monotherapy curves express the fractions of unaffected control populations, while increasing curves express the fractions of affected control populations.

Bliss independence model is formulated for the fractional responses f (“fraction affected”), where the predicted response f₁₂ at dose combination (d₁, d₂) is defined as:

f₁₂(d₁, d₂) = f₁(d₁) + f₂(d₂) − f₁(d₁)f₂(d₂),

with $$f_1(d_1) = \frac{y\left(d_1\right) - b}{m_1 - b} = \frac{1}{1 + \left(\frac{\operatorname{EC50}_{1}}{d_{1}}\right)^{|h_{1}|}}$$ $$f_2(d_2) = \frac{y\left(d_2\right) - b}{m_2 - b} = \frac{1}{1 + \left(\frac{\operatorname{EC50}_{2}}{d_{2}}\right)^{|h_{2}|}}$$ In the classical Bliss independence model, marginal models share baseline and maximum response.

To allow the compounds to have different maximal responses, the fractional responses are rescaled to the maximum range (i.e. absolute difference between baseline and maximal response). Then the predicted response is defined as:
$$ y = b + (m_{max}-b) \left[ \tilde{f_1}(d_1) + \tilde{f_2}(d_2) - \tilde{f_1}(d_1)\tilde{f_2}(d_2) \right], $$ where m_max is one of m₁ or m₂, for which the value of |m_i − b| is larger, and $$ \tilde{f_i} = f_i\frac{m_i-b}{m_{max}-b}~~\text{for}~~i = 1, 2.$$

This implementation of Bliss independence supports both compounds with decreasing and increasing monotherapy profiles. However using one compound with a decreasing profile and another with an increasing profile in combination is not supported.

Alternative Loewe Generalization

An alternative generalization of Loewe Additivity for the case of different asymptotes can be defined as a combination of Loewe and HSA approaches as follows:

In a classical Loewe equation, predicted response y at a given dose combination (d₁, d₂) can be found by solving the equation:

$$ \frac{d_1}{D_1(y)} + \frac{d_2}{D_2(y)} = 1, $$

where $D_i(y) = \operatorname{EC50}_{i}\left(\frac{y-b}{m_i-y}\right)^{\frac{1}{|h_i|}}$, for i = 1, 2 is the dose of the i-th compound that gives response y. Note that here D_i is properly defined only if y is between b and m_i.

For the case of different asymptotes, say when y > m₁ (increasing curve) or y < m₁ (decreasing curve), we set D₁(y) = +∞, so that the y is determined from the equation d₂ = D₂(y), replicating what is done in the HSA approach.

Calculation procedure

In order to evaluate any of the null models described above, the fitSurface function will use the monotherapy parameter estimates from the previous step. The idea is if there are synergistic or antagonistic effects, then administration of both compounds will lead to important deviations from what combined monotherapy data would suggest according to the null model. Routines within fitSurface function do essentially the following.

1. Find occupancy for each combination of doses by solving the additivity equation of the classical Loewe model. This step does not require knowledge of the baseline or maximal response for either of the compounds. Occupancy solution is also reported in the HSA model case although occupancy plays no role in such a model.
2. Compute the predicted response based on the above described response equations and the previously computed occupancy rate for each dose combination.
3. If desired, the function will then calculate the selected statistic to evaluate the deviation of the predictions from the desired null model.

meanR

The meanR test will evaluate whether the null model globally fits well the observed data. It is derived using a lack-of-fit sum of squares technique. In particular, the test statistic is given by $$ \operatorname{TestStat} = \frac{R^{T}\left(D + C_{p}\right)^{-1}R}{n_{1}\sigma^{2}} $$ Assuming that residuals from the generalized Loewe model are normally distributed, it can be shown that this statistic follows an F_n1, df0 distribution under the null. If these assumptions are not satisfied, the null distribution can be approximated by bootstrapping.

maxR

The maxR test evaluates whether the null model locally fits the observed data. In particular, it provides a test score for each off-axis combination. Based on the sign of this score, it can be determined whether synergy or antagonism is more likely and a formal test can be constructed. Under the null hypothesis of no lack-of-fit and normally distributed effects, max |R^T(D + C_p)^−1/2|/σ ∼ max |Z₁, …, Z_k| where Z_j ∼ N(0, 1). More particularly, the test statistic for the k-th off-axis dose combination (d₁, d₂) is computed as TestStat (d₁, d₂) = [|R^T(D + C_p)^−1/2|/σ]_k where [⋅]_k indicates the k-th coordinate. This test statistic is then compared either to the null distribution based on normal approximation or a bootstrapped approximation.

Adapted meanR

The adapted meanR test uses two separate variance estimates for (a) the monotherapies (= σ₀²) and (b) the dose combinations (= Σ₁, a diagonal matrix). The notation for both unequal as model will be the same, but the estimation of Σ₁ will be different. The variance of the monotherapies σ₀² is estimated as σ² above by taking the MSE of the null model. The test statistic is given by:

$$ \operatorname{TestStat} = \frac{R^{T}\left(\Sigma_{1}D + \sigma^{2}_{0}C_{p}\right)^{-1}R}{n_{1}} $$

1. unequal method: The variance for the dose combinations is estimated by taking the variance in each dose combination and then taking the mean of all these variances, thus Σ₁ = σ₁²I_n1. The downside of this method is that the variance for all combinations is assumed to be equal. In reality the variance often depends on the mean effect. This is taken into account in the model method.

2. model method: In this method the diagonal elements of Σ₁ are no longer estimated as single number but rather as a vector of variances. Each off-axis point has now its own variance. A linear model is fitted on the original dataset, modeling the variance of each off-axis point as a function of its mean effect. The estimated model parameters are then used to predict the variance for the corresponding mean effect measured for that dose combination. These predicted variances are placed in the diagonal of Σ₁. Modelling the variance with a linear model may require a transformation, to achieve a better fit and to avoid negative variances being modelled. A log-transformation often makes a good impression.

Adapted maxR

The same approach is taken for the adapted maxR test statistic. Instead of using one estimated variance for both on- and off-axis points, two separate estimates are used. The estimates for Σ₁ are different depending on the method used (unequal or model). The methodology of estimating the variance is the same as was described in the “Adapted meanR” section above.

The maxR test becomes

max |R^T(Σ₁D + σ₀²C_p)^−1/2| ∼ max |Z₁, …, Z_k|

where Z_j ∼ N(0, 1). In particular, the test statistic for the k-th off-axis dose combination (d₁, d₂) is computed as TestStat (d₁, d₂) = [|R^T(Σ₁D + σ₀²C_p)^−1/2|]_k where [⋅]_k indicates the k-th coordinate.

Bootstrapping under unequal variances

In case of the unequal variance assumption, the bootstrap proceeds as before, with the off-axis residuals being pooled and resampled. With the model assumption, the resampling is more complicated, as the residuals are no longer exchangeable. One option is to rescale the observed residuals according to the mean-variance model (i.e. dividing them by their standard deviations), resample from this pool of standardized residuals, and then scale back to the true variance (by multiplying by the standard deviation). Yet this approach has proven to be unstable as it leads to extreme observations. An alternative (the default) is to generate zero-mean normal data with the modelled variances (see the rescaleResids argument in fitSurface()).

Advantages of unequal and model methods compared to assumption of equal variances

The assumption of equal variances between monotherapies and off-axis dose-combinations fails to control the type I error rate around pre-specified level, when the variance of off-axis points increases (natural variance or outliers). This results in false positive synergy calls when in reality there were none. Both the unequal and the model methods control the type I error rate far better, with slightly better results obtained by the model method.

Furthermore, the sensitivity and specificity of the maxR test statistics are higher with the methods assuming variance heterogeneity compared to the methods where equal variances are assumed.

Confidence interval

Let e_i denote the true effect size on off axis point i, and call E_i its estimate based on the data. Relying on the asymptotic normality of the estimator, an approximate (asymptotic) confidence interval would be formed as the set:

with ŝ_i the estimated standard error of E_i, and z_α/2 the 1 − α/2 quantile of the standard normal distribution.

We know, however, from the meanR and maxR tests that the asymptotic distributions provide poor approximations. Therefore, we use the bootstrap here too to build the confidence intervals.

For every bootstrap instance, bootstrap observations are sampled for on- as well as off axis points. For the on-axis points, a parametric bootstrap based on the estimated monothereapy curves is used, as for the calculation of the meanR and maxR statistics. Based on these on-axis bootstrap samples, new monotherapy curves are fitted with resulting residual variances bootstrap variance MSE₀^b, and corresponding response surfaces with expected outcomes R̂_i^b are derived. For the off-axis points, with n_i replicates at a given point, A_ij = R_ij − R̄_i are the pointwise residuals j = 1, ..., n_i. Here R_ij is the observed outcome and R̄_i, the average outcome at point i, serves as an unbiased estimator of the true response. Note that these residuals are different from the residuals E_i = R̄_i − R̂_i used to construct the test statistics; for A_ij the departure with respect to the mean outcome at that off-axis point is used. These residuals are resampled with replacement from the observed residuals, possibly using rescaling as explained below. The resampled residuals are then added to the estimated effect sizes to obtain bootstrapped observations R_ij^b = R̄_i + A_ij^b, with b = 1, ..., B denoting the bootstrap instance. This leads to the bootstrap effect sizes E_i^b = R̄_i^b − R̂_i^b and test statistics $\left\vert \frac{E^b_i -E_i}{\hat{s}^b_i} \right\vert$, using the standard deviations ŝ_i^b from the bootstrap.

Over all bootstrap instances, we then find the distribution of

Call t_α the threshold such that

Finally, we find for every off-axis point the confidence interval as the collection:

Standard error

As an estimate of the standard error we use

with F^b equal to MSE₀^bD, MSE₁^bD or S^b, depending on which variance model is used. This standard error resembles the variance estimators of the meanR and maxR statistics discussed there.

By using only the diagonal elements of F^b, we ignore the covariance between test statistics. We will use the studentised-range concept (which is also central to Tukey’s method for multiple testing) for controlling the family-wise error rate (FWER). For the covariance matrix C_p we use the one estimated for the observed data. This matrix turns out to be quite stable over the bootstrap runs; moreover the calculation of the matrix for each bootstrap sample would imply a time-consuming nested bootstrap procedure. On the other hand, MSE₀^b is re-estimated with each bootstrap sample. Also F^b is re-estimated based on the bootstrapped data.

Single effect measure

Researchers may also want to have a single measure for the strength of the synergy for a single experiment, e.g. in view of ranking compound combinations according to their synergistic effect. For this we calculate the average of the pointwise off-axis effect sizes. This is estimated as $\bar{E} = \frac{1}{\sum_{i= 1}^{n_1}n_i }\sum_{i=1}^{n_1}\sum_{j= 1}^{n_i}(R_{ij}-\hat{R}_i)$. It may be considered as a measure of the “volume” between the expected and observed response surfaces, similar to the “integrated synergy” of . Note that this effect size may cause synergistic and antagonistic points to cancel out against one another, but we believe this scenario is unlikely. As before we would like a measure of imprecision for this effect size.