M o n t e Carlo comparison of rival experimental designs for t w o - a g e n t c o m b i n e d action studies

comparison of rival experimental for two-agent combined OBJECTIVE: The combined action of two or more chemotherapeutic agents and/01· biological agents can be quantitatively described wilh empirical multidimensional concenlration-effecl response surface models. This intuitive statistical approach provides a framework for suggesting experimental designs for in vitro. in vivo and possibly clinical experiments of agent combinations. Five rival 32-poinl cxperimenlal designs for in vitro continuous re ponse two-agent combined action studies were compared using Monte Carlo simulation. D ESIGN: The designs were: factorial; central composite: one-ray in duplicate; four-ray; and D-optimal. SETTING: Datasets were simulated by generating ideal data with the authors' flagship lwo-agcnl combined action model. which includes six parameters: the control sunrivaL Econ=lOO (where Econ is lhc full range of response that can be aJfecled by the drug); median effective concentrations. /C50, 1=10. /Cso.2= 1 for drug l and drug 2, 1·cspectively; slope parameters. m1 =- 1. m2=-2 for drug l and drug 2. respectively; and the interaction parameter, a=l or a=5. For each design, for each of four types of error (absolute. relative wiU1 1% coefficient of variation [cv]. relative with 10% cv. and relative with 10% cv plus a noise conslanl of 1% of Econ). for each of two values of the true a (1. 5). 500 Mon le Carlo datasets were generated. and then flt via weighted nonlinear regression wilh lhe flagship model. MAIN REsULTS: For the a parameter. for relative error-containing datasets. the D-oplimal designs had lhe smallest variances. CONCLUSION: The counterintuitive D-optimal designs may be useful for studies in which the experimental units arc relatively precious. and frugal designs are essential. In addition. il may be fruitfu l lo add lhe D-oplimal design points lo standard experimental designs. ct·association.


THREE-DIMENSIONAL COMBINED ACTION MODELS
Our f1agship combined action model is given as Equation 1. Equation 1 was derived with an adaptation of an approach suggested by Berenbaum (10). with the assumption of Equation 2. the Hill model (11). as the appropriate model for each agent alone.
In Equations 1 and 2. E is the measured effect (response). D is the concentration (or dosage) of drug. Econ is the full range of response that can be affected by lhe drug, /C50 is lhe med ian effective dosage (or concentration) of agent. and m is a slope parameter (additional subscripts. 1 and 2 . refer to agents 1 and 2, respectively). When m has a negative sign, the curve falls with increasing agent. concentration; when m is positive. the curve rises with increasing agent concentration. The interaction parameter is a. When a is positive. Loewe synergism is indicated, when a is negative. Loewe antagonism is indicated , and when a is 0 , Loewe additivity is indicated. (This terminology is the central element of a recent consensus on combined action nomenclature and concepts [12]. The use of 'Loewe· as an adjective stresses Lhe connection to the ·no interaction· reference model inherent in the classi-52A cal isobologram approach lo interaction assessment pioneered by Loewe [13].) The magnitude of a indicates the intensity of U1e interaction. Thus. although Equation I is not the model for Loewe synergism (or Loewe antagonism). it is a model for Loewe synergism (or Loewe antagonism).
The derivation of Equation 1, the flagship equation for two-drug interactions, is provided in detail in Greco el al (6) . Equation 1 allows the slopes of Lhe concentration-effect curves for the two drugs lo be unequal. Il is this key feature that distinguishes Equation I from many other response surface models used by others to describe drug interactions (eg, 14). Since Equation I is in unclosed form, a one-dimensional bisection rool finder (eg, 15) is used to calculate E for sinrnlations. Equation 1 was not derived from biological theory. rather it is an empirical equation which often matches the shape of real data (eg, 2,6-8).

STATISTICAL EXPERIMENTAL DESIGN
The main decisions that must be made regarding experimental design are: first. where to choose Lhe concentrations; second, numbers of replicates: and Lhird, numbers of experiments. These seemingly simple questions have spawned many full careers for statisticians, who have delved deeply into Lhem to reveal their inherent complexity. The adoption of a response surface paradigm for the assessment of combined action of agents facilitates the understanding and use of formal statistical experimental design. First. Lhe experimenter must decide whether he/she is in an exploratory or a confirmatory mode. Screening expe1iments (exploratory mode) should first include, for each agent individually, agent concentrations which span lhe anticipated response region. Logarithmic spacing of lhe concentrations over a thousand -fold to a million -fold range is probably necessary. depending upon U1e previous knowledge of the researcher about Lhe concentrationeffect behaviour of Lhe compound. After the individual agent concentration-effect. curves are well characterized, a combination experiment should be conducted    that repeats single agent data points, and that includes a set of combination points. A full factorial (checkerboard) design (Table 1). a single ray (fixed ratio) design (Table 3) or a multiple ray design (Table 4). all with logarithmically spaced concentrations, may be appropriate. If a complex three-dimensional concentrationeffect surface is anticipated. then the entire interesting region of agent 1 and agent 2 concentrations should be sampled. either with a checkerboard or multiple ray design. However, if a well beh aved three-dimensional concentration-effect surface is anticipated, and the specific combination being studied is only one of many candidates being screened. then a single ray may be sufficient. Composite designs consisting of a checkerboard and some rays might also be used. Of course. if the in tended data analysis approach is firmly tied to a particular design, then that design will have to be used. After comp leting the analysis of the fi rst mixture experiment in exploratory mode, the researcher may want to switch to confirmato ry mode . The repeat of the combination experiment may use the same design as in the exploratory experiment, but probably the knowledge gained from the first run will help to refine the design for the second run. If a complex three-dimensional concentration-effect surface was found in the exploratory experiment, then agent concentrations in the interesting regions of the surface should be accented in the confirmatory experiment. Increasing the n u mbers of replicates probably also will be necessary. If a simple three-dimensional concentration-effect surface was fo u nd in the exploratory experiment. ie. on e with pure Loewe synergism or Loewe antagonism,        then a design lhat facilitates the estimation of parameters with the smallest variance may be appropriate, A D-oplimal design (eg, 4) may be indicated with many replicates. Inlereslingly, the number of design points in a D-opUmal design (Tables 5-10) is equal lo the number of estimable parameters. For example, if one assumes that Equation l, which contains six parameters, will adequately describe lhe three-dimensional combined action concentration-effect curve, then a D-optimal design will include only six design points, with or wilhou l replicates. The D-optimal designs may, al first. seem to be very strange and potentially noninfom1ative. This type of frugal experimental design may have great potential for animal and human experiments, in which the experimental units are very dear. However, because D-optimal designs are very counterintuitive. a Monte Carlo computer simulation study was designed to compare several rival experimental designs, including D-optimal designs. for combined action studies.

METHODS
Monte Carlo simulations: For lhe Monte Carlo simulation, Equation 1 was assumed lo represent the combined action of two agents, with realistic parameter values: Econ= 100: JC50. L = 10; m, =-1; JC50.2= 1; m2=-2; and a= 1 (slight Loewe synergism) or a=5 (large Loewe synergism). Ideal continuous data were generated by inserting values of D1 and Dz for a specific design into Equation 1 and calculating E for each data point with a bisection root finder (15). The errors that were added lo the ideal data lo yield error-containing data were of four types: small absolute (SA); small relative with a 1 % 54A  (19). Initial parameter estimates for Econ were generated by laking an average of effects al the (0,0) design point. Initial parameter estimates for IC50.1. m1. and for JC50.2. m2, were generated by fitting single agent data for drug 1 and drug 2, respectively, with a linearized transformation of Equation 2 via weighted linear regression. Initial parameter estimates for a were generated by rearranging Equation 1 lo isolate a on the left side of the equals s ign, then for each combination point, plugging in raw data and the other five initial parameter values lo solve for a, and finally, calculating an average a for lhe set of combination points. The      software written by U1is group in the C programming language. and run on MSDOS-compatible 80486-and 80386-based microcomputers. Rival experimental designs: The 10 experimental des igns that were compared are shown in Tables 1 lo 10. All rival designs contain 32 data points. Table l shows a factorial or checkerboard design. with logarithmically spaced concentrations. Table 2 shows a central composite design. Table 3 shows a single fixed ratio or single ray design , in which the single fixed ratio is 10: 1 for Di:JJ,;. (the ratio of their lC50 values). Table 4 shows a  Tables 5 to 10 show D-opUmal designs and are directly dependent upon tl1e predicted values of U1e six parameters of Equation 1 and on the four different error structures. For each D-op1.imal design, the first five optimal concentration pairs (D1 . Lh) (those not associated with the estimation of the interaction parameter , a) were generated with formulas and concepts published by Bezeau and Endrenyi (20) for the precise estimation of the three single agent parameters for Equation 2. A discussion of the generation of 0 -optimal designs for tl1e sixth optimal concentration pair (the one associated willi the estimation of a) for combined action studies with numerical function m in imization methods has been published by the present group (4) . However, the D-op1.imal designs in Tables 5 to 10 we re generated with faster algorithms, which relied h eavily on partial analytical solu Lions of the determinants of llie variancecovariance matrix of parameters. for which the d etails will be published elsewhere. The D-optimal formulas were coded in the Mathematica mathematical programming language (21) . All concentrations calculated to be greater than 1000 were capped at 1000 to keep the D-optimal designs realistic. Figures 1 to 4 show box p lots of the distributions of the sets of estimated a parameters. Each figure is for a different error type. Figure 1 is for SA, Figure 2 is for SR, Figure 3 is for LR and Figure 4 is for LC. These figures were made with the Sigma Plot graphics package (22) . Table 11 lists the variances for the set of 500 a es timates for each design for each type of error (SA, SR. LR or LC) with each true a (1 or 5). For each type of error and each value of the tru e a, the designs are listed in the order of increasing variance for the sets of 500 a estimates. Each set of 500 estimates was also divided into 20 consecutive subsets of25 values, and the variance for each subset was calculated. Then these sets of 20 variance values were used to make pairwise comparisons of variance among all rival models in a group using llie Wilcoxon rank sum test via the SAS statistical package (23) .

RESULTS
Since each group of designs involved 10 comparisons, the type I error rate of 0.01 was chosen for making decisions of statistical significance for pairwise design differences , for an overall Bonferroni conservative type I error rate of 0 . 1 per group. For each group, llie upper designs willi the smallest variance form a subgroup, for  Table 11. First, for every group, the D-optimal designs were always among the design subgroup with the smallest variance. The superiority of the D-optirnal designs over a ll of the other designs is clearly seen for the SR, LR and LC error types, but is Jess clear for SA, especially for the case in which the true a= 1. Second. the FR4X design (Table 4) was always in the subgroup with the largest variance. Third , there was a tendency for the FR2X design (Table 3) to be a relatively precise design and for 5X5 (Table 1) (Table 7). OPTD (Table 8). OPrE (Table 9)  and LC error types, respectively. these designs resulted in a significant bias.

DISCUSSION
Among the five rival experiment.al designs, the seemingly counterintuitive D-optimal designs appear to result in the smallest. variance for the interact.ion parameter. a. From this result., one can infer that the D-optimal designs may be superior for assessing quantitatively the nature and intensity of combined action of two agents. The convenient single fixed ratio in duplicate design (FR2X, Table 3) was also very good. and ranked second.
Formal statistical experimental design often includes an interesting paradox: to design an experime nt well. one has to know the final answer well. However. if the final answer is well known, ie, both the correct model and the true model parameters. then one would not. h ave to conduct the experimenl. This paradox is solved with sequential experimentation; each experiment in a sequence provides betler information for the