CISNE: An accurate description of dose-effect and synergism in combination therapies

The precise determination of dose-effect curves and the combination effect of drugs is of crucial importance in the development of new therapies for the most dreadful diseases. We have found that the current implementations of the theory of Chou et al. are not accurate enough in some circumstances and might lead to erroneous predictions of synergistic or antagonistic behaviour. We have identified the source of inaccuracies and fixed it thereby improving the accuracy of those methods. Here we explain the main features of our approach and demonstrate its higher accuracy as compared to the standard methods. Therefore, this new implementation might have a huge impact in the reliability of future research on new Combination Therapies.

The precise determination of dose-effect curves and the combination effect of drugs is of crucial importance in the development of new therapies for the most dreadful diseases. We have found that the current implementations of the theory of Chou et al. are not accurate enough in some circumstances and might lead to erroneous predictions of synergistic or antagonistic behaviour. We have identified the source of inaccuracies and fixed it thereby improving the accuracy of those methods. Here we explain the main features of our approach and demonstrate its higher accuracy as compared to the standard methods. Therefore, this new implementation might have a huge impact in the reliability of future research on new Combination Therapies.
Combination Therapies have been a valuable approach for the treatment of diseases since the origins of medicine. Traditional Chinese medicine has been using mixtures of natural herbs for more than 2000 years 1 . Nowadays, the use of adequate drug combinations allows to increase beneficial effects using lower doses of each constituent and therefore to reduce their adverse effects, as well as to minimize the induction of drug resistance 2 . It can also provide selective synergism against a particular target. For these reasons, Combination Therapies are widely used for the treatment of the most dreadful diseases, such as cancer 3,4 or AIDS 5,6 .
Potential drug combinations are frequently tested first in cell cultures in preclinical environments, whereby the successful combinations are subsequently transferred to animal experimentation, clinical testing and clinical practice. However, the transition from a biologically controlled environment to the clinical practice often results in failures 7,8 . This lack of success is possibly caused by the difficulty to reproduce the overall physiological conditions in the laboratory. However, they would also occur if both drugs act on the same cellular sub-clone. We finally note that some of those failures might be originated by an insufficiently accurate statistical analysis of the cell cultures that may provide false positives or negatives in specific situations. The goal of the present article is to eliminate this last source of errors by implementing drastic improvements of the referred statistical analysis Computational methods can be applied to assist on the determination of the effect of Combination Therapies. First, they can be used in an effort to predict, previous to experimental verification, potential drug combinations [9][10][11] . Also, they can be used after experiments to analyse statistically the measured data in order to quantify the presence or absence of synergism 12,13 , that is the subject matter of the present work. Currently available methods 12,13 to quantify the synergistic effect of drug combinations from experimental data predict accurate dose-effect relationships in many cases ( Fig. 1(a1)), which lead to a correct prediction of the synergism results ( Fig. 1(a2)). However, we have found that they fail to reproduce the experimental data in quite a few others ( Fig. 1(b1)). This lack of precision may result in an inaccurate or even erroneous description of the synergistic properties of drug combinations ( Fig. 1(b2)). Therefore, there is an urgent need to develop a method that improves the accuracy of the dose-effect description, which then enables us to identify synergism effects in experiments with a minimal uncertainty.
Although drug synergism might seem an intuitive concept, its formal definition is not so simple and it has been a focus of controversy during the past century [14][15][16] . A combination of drugs is synergistic when the combined effect is larger than the additive effect of each individual drug. Similarly, the combination is antagonistic when its combined effect is smaller than the additive effect of each individual drug. Trouble arises however from the definition of the additive effect. Let us assume that two drugs A and B have an effect that we quantify with the numbers E A and E B (e.g.: the number of affected cells in a system). We define tentatively the additive effect of the two drugs as E A + E B . However, this definition would be true only if the dose-effect relationships of each drugs were linear, which is usually not the case. The modern definition of drug synergism [16][17][18][19][20][21] is based on the work by Chou et al., who define the median-effect equation (MEE) [22][23][24] as the universal equation for the non-linear dose-effect relationship. This definition is based on physicochemical processes, and its validity does not depend upon factors that are specific of each system such as the mechanism of the enzymatic reaction or the kinetic nature of the inhibitor. From that, they define the Combination Index (CI) 25,26 to quantify the synergistic effect of a drug combination. A given drug combination is said to be additive when its CI equals 1. The degree of synergism (or antagonism) is measured by the deviation of CI from 1, the smaller (the larger) the CI, the stronger the synergistic (antagonistic) effect of the combination is. Therefore, in order to analyse synergism in a drug combination, a full and accurate description of the dose-effect relationship of each drug separately, as well as of their combination at a fixed drug ratio, is needed 16 .
In this work, motivated by our need to obtain results as accurate as possible from synergism studies, we have analysed the numerical problems of the currently available methods such as calcusyn 12 and compusyn 13 . We present our new implementation of the theory of Chou et al., which solves these problems, thus increasing the precision and reliability of previous methods. We have made the new implementation freely available through the use of our new program CISNE (Code for the Identification of Synergism Numerically Efficient) 27 . Therefore, our new method keeps the generality of a theory that has been widely validated and applied to a variety of biological systems 28 , but, at the same time, it represents an improvement on the implementation of that theory. This is because it allows to determine accurately possible synergistic effects in biological systems within a wide range of characteristic parameters where the previous methods failed.

The theory of synergism
In a dose-effect experiment several doses x of a drug are applied to a population (e.g.: a cell culture), and their effect is determined by measuring the fraction of the population affected by the drug after a certain period of time, f a . Current practice actually plots the fraction unaffected f u = 1 − f a , as a function of x, as we show in Fig. 1.
where D and m are parameters that depend on each given drug. D is the dose at which f a = 0.5, whereas m determines the shape of the dose-effect curve. Here m = 1, >1 or <1 indicate hyperbolic, sigmoidal or flat sigmoidal curves, respectively. Examples of these different curves are given in Fig. 2(a). Synergism is characterized by the combination index CI, that is always given as a function of f a . We discuss in some detail here the combination index for two drugs, although the formalism can be extended for N drugs. The combination index for drugs A and B combined in a ratio R A :R B is given by Figure 1. Example of a dose(in arbitrary units a.u.)-effect(unaffected fraction f u ) data (green diamonds) for which currently available methods (red lines) reproduces accurately the shape of the dose effect curve (panel (a1)) or completely fails to reproduce it (panel (b1)). This can lead to a correct determination of the synergism (panel (a2)) or to a completely wrong result (panel (b2)). The exact result (black dashed lines), obtained from exact data previous to introduce a experimental deviation, is also shown. Insets of panels (a1) and (b1) show the linear regression used to fit the data.
, are the doses of drug A or B that produce an effect f a alone and in the combination, respectively. Solving for x in Eq. (1) yields x = D(f a /f u ) 1/m . This expression for A, B and A + B can be inserted into Eq.
where m A,B,A+B , D A,B,A+B are the MEE parameters obtained for drugs A or B alone or in the combination A + B.
Notice that a not sufficiently accurate determination of any of these parameters may lead to erroneous estimates of CI.
We discuss now how popular implementations of Chou's theory 12,13 determine the important parameters m and D of a given drug and give some hints on the source of their inaccuracies. First, Eq. (1) is converted to the linear relationship.
where y′ = log(f a /f u ), x′ = log(x) and b = −mlog(D). Then the experimental data points are fitted to Eq. (4) using a linear least squares fit. This procedure delivers estimates for m and b. The remaining parameter is obtained from D = e −b/m . This approach is accurate whenever the experimental data points obey closely the theoretical dose-effect curves. In practice, the experimental dose-effect data obtained from biological samples always contain some non negligible dispersion whose impact need not be negligible.

Why current implementations fail
We illustrate the impact of actual experimental data dispersion via an example. Imagine an experiment in a culture where the measured fraction unaffected of the population f u leads to the any of the curves plotted in Fig. 2(a).
Here we give the same standard deviation σ to each experimental data point. We plot now the same data points in a logarithmic scale in Fig. 2(b), so that the axes represent x′ and y′. The figure demonstrates that a constant σ for f u leads to a non constant standard deviation σ′ for y′. Importantly, σ′ increases hugely at the sides of the curve specially if m is sufficiently large. This pathology is made more evident in Fig. 2(c), where we show how y′ and σ′ depend on f a . The figure shows that σ′ increases and eventually diverges as we move away from f a = 0.5. This means that the linear least squares fit procedure gives artificially more weight to those experimental points having f a closer to 0 or 1. Therefore, the results are artificially biased. We will show below that this artifact can harm substantially the quality of the results even leading to false positives or false negatives verdicts of synergistic behaviour.
Our new implementation in CISNE avoids the above pathology from the outset. We take advantage of modern curve fitting algorithms, that obtain m and D directly from a nonlinear least squares fitting 29 of equation (1). This approach presents two main advantages: 1. All the data points share the same weight, that is independent of the value of f a thus avoiding any bias. 2. Experimental data points having f a ≥ 1 or f a ≤ 0 happen experimentally due to experimental fluctuations and carry information that should not be ignored. The linear least squares fitting does ignore them and misses relevant information while the nonlinear least squares fitting does not present this limitation, and therefore no data information is lost. We show now two examples where the current implementation fails while CISNE reproduces the correct behaviour. In the first example, two drugs A 1 and B 1 show purely additive behaviour (CI = 1), however the current implementation predicts strong synergism (  CI 1). We generate dose-effect data points numerically for the two drugs and for their combination A 1 + B 1 in a 1:1 ratio, that we show in the upper left panel in Fig. 3. Notice that we have introduced random deviations in the data points following a standard deviation σ = 0.05 to try and mimic experimental data. The generated dose-effect data is given in Table 1 of the Additional Information. The upper left panel of Fig. 3 shows the results of the dose-effect fitting curve with the standard (red lines) and the CISNE (blue lines) methods, as well as the exact curve before introducing the experimental deviations (black dashed lines). The panel shows that the current implementation reproduces well some of the dose-effect relationships but it fails for others, while CISNE fits accurately the experimental data points. More important however are the predictions for the combination index obtained from the fitted curves, that are shown in the upper right panel of the figure, together with their 95% confidence interval. Here CISNE reproduces quite correctly the additive behaviour of the combination while the standard method predicts a strong synergistic effect leading to a false positive.
In the second example, two drugs A 2 and B 2 show a strongly synergistic behaviour (CI = 0.6). We generate dose-effect data points for the two drugs and for their combination A 2 + B 2 in a 1:1 ratio, where we include random deviations as before. The dose-effect curves are plotted in the lower left panel of Fig. 3 and the corresponding data are given in Table 2 of the Additional Information. The panel also shows the curves fitted with the standard method as well as with CISNE. As in the former example, we find that CISNE fits the experimental data points much more accurately. We use these curves to determine the combination index. We plot our results in the lower right panel of Fig. 3 together with their 95% confidence interval. We find that CISNE reproduces accurately the synergistic behaviour while the standard method predicts additivity or even antagonism leading to a false negative.
Notice that the wrong results with the standard method are a consequence of the inaccurate description of some of the dose-effect data by the fitted curves, especially in the central range of f u , which is not visible in the  linear fittings shown in the insets of Fig. 3 (where r > 0.96 in all cases indicating an accurate linear fitting). In the first example the fitted curve underestimates the effect of drug A 1 (the fitted curve is shifted to the right respect to the dose-effect data points), but it overestimates the effect of the drug combination A 1 + B 1 (the fitted curve is shifted to the left respect to the dose-effect data points). This combined effect leads to a large overestimation of the effect of the combination of drugs A 1 + B 1 respect to their individual effect, and in this case to a false positive. Similarly, in the second example the effect of drug B 2 is overestimated and the effect of the drug combination A 2 + B 2 is underestimated, leading to a large underestimation of the effect of the combination of drugs A 2 + B 2 respect to their individual effect, and, therefore, to a false negative in this case.  Table 2. Generated dose x and fraction affected f a data of a combination of drugs A 2 + B 2 with synergistic behavior.

More computational experiments
In order to evaluate the precision of the new methodology implemented in CISNE compared to the standard procedure more rigorously, we have analysed the MEE parameters and CI obtained from experimental data generated computationally by means of Monte Carlo simulations 30,31 . In these simulations we consider the data points shown in Fig. 2(a) Fig. 4 we show the obtained values of D and its standard deviation σ(D). Notice that σ(D) is a measurement of the confidence interval of D from each of the experimental data sets separately, and it is not related to the standard deviation of the gaussian-like representations of D in the left panels of Fig. 4, which contain the information from all the experiments together. As expected, the obtained results for D follow a normal distribution centered at its expected value E[D] = 10. For the flat sigmoidal curve (m = 0.5), the results of the standard procedure and of the new CISNE method are very similar. This makes sense, as the selected data shown in Fig. 2(a) presents values of f a which do not come very close to 0 or 1, and therefore almost all our data points stay in the central region of Fig. 2(c), where the change in the standard deviation of the points is not very important. However, the benefit of the new approach becomes apparent for hyperbolic (m = 1) and sigmoidal (m = 2) curves. This is accompanied by smaller values of σ(D). Therefore, CISNE allows to obtain values of D much closer to the real value, and with better confidence intervals, especially for m ≥ 1. Figure 5 shows the obtained results for m and its standard deviation σ(m). The difference between the standard and the CISNE approach becomes even more clear in this case. CISNE presents gaussian distributions centered at the corresponding expected values E[m] = 0.5, 1 and 2, as expected. On the other hand, the standard procedure leads to a distribution with a maximum at m slightly lower than 1 for the hyperbolic curve and much lower than 2 for the sigmoidal curve. Moreover, the standard deviations σ(m) obtained with the standard procedure are lower than the ones obtained with the CISNE method (notice that in both cases σ(m) increases with E[m] = 0.5, but its relative value σ(m)/E[m] = 0.5 remains almost constant). Therefore, meanwhile CISNE reproduces accurately the results of m, the standard procedure can not only lead to wrong results, but could also lead to wrong confidence intervals which do not correspond to reality, leading to a false impression of accuracy in the results.
In order to analyse the accuracy of both methods in the determination of the CI, we perform blocks of 3 simulations with the same methodology as before, and we assign the first one to drug A, the second to drug B and the third to the drug combination A + B in a 1:1 ratio. As the 3 simulations have the same expected values of m and D, the obtained effect should be just additive and CI(f a ) = 1. However, deviations from this value are expected from the effect of the experimental fluctuations. We repeat these blocks of simulations 100000 times for E[m] = 0.5, 1, 2 and the results of CI(f a = 0.5) and its standard deviation σ(CI(f a = 0.5)) are summarized in Fig. 6.   The results obtained with the standard and the new CISNE methods are similar for E[m] = 0.5, but the advantage of the new method is clear for E[m] = 1 and E[m] = 2. We show in Table 3 the percentage of the simulations that are inside an interval where we can consider synergistic (f a < 0.8), additive (0.8 < f a < 1.2), and antagonistic (1.2 < f a ) effect. CISNE allows to strongly reduce the number of false positives of synergistic and antagonistic effect, especially for E[m] = 2, where the standard method leads to more than one third of simulations indicating non-additive effect, meanwhile this value is reduced to less than 3% of the experiments with the CISNE method. Notice that, as we are using a theoretical additive combination, wrong results appear as false synergistic or antagonistic results. However, the increase in the precision of the obtained CI with CISNE is not limited to reduce false positives, as shown in previous sections. If we consider a synergistic combination, the standard method can lead to wrong results of additive (or, even worse, antagonistic!) effects, and CISNE allows to minimize the probability of obtaining these wrong results as well.

Conclusions
Combinations Therapies that present a synergistic effect is a valuable resource for the treatment of the most dreadful diseases. We have demonstrated in this article however that the accuracy whereby the available schemes determine whether a drug combination is synergistic or not can easily lead to false positives and false negatives. We have also identified the main source of inaccuracy and fixed it. Our new implementation allows us to determine dose-effect relationships more accurately, that conveys a more precise determination of the CI. In particular, our Monte Carlo simulations show that the number of false positives that can be obtained from the mathematical analysis can be reduced by a factor of 10. We find that the defects of the standard method become specially damaging whenever the experimental data deliver values of the fraction of the population affected f a larger than 0.9 or smaller than 0.1. Although our new method presents advantages for any range of experimental data, we stress that the use of our new method becomes especially relevant in these cases.