1. Introduction

Combining drugs in biological systems is a common practice, and the effect can be one of three types: additive, synergistic, or antagonistic. Mathematical modeling is often used to describe the relationship of two drugs. One such model, referred to as the most suitable, is the Loewe additivity model under the assumption of no interactions in the model as in[1]. This model refers to a combination of two concentrations of drug A and drug B being additive and equivalent to an interaction effect. Reference[2] has established an extension of a model free test for synergy in multiple drug combinations. Reference[3] used additive models generated from stochastic differential equations for combinations of two aesthetic agents. In this paper, we use the later approach for making the model for data involving two isomers applied in combination. We propose the necessary and sufficient conditions to validate the reference[2] test.

It is well-known that the flow of a chemical in the human body and its concentration at various times can be described by differential equations. After an initial dose of the chemical is injected into the system, some of it will escapefrom the system and some will remain to decay over time. Observation of this phenomenon suggests that the concentration of the chemical in the system at time t can be modeled by a linear differential equation of the form

(1)

Here, is a positive constant that denotes the relative rate of elimination of the chemical from the body and is a generally decreasing function. The model is based on one- compartment linear model with infusion. Kinetic models of this type have been extensively studied (see e.g. the book of[4], papers of[5],[6] and many other works).

As pointed out in[7], the function may be subject to random fluctuations from a variety of physical and physiological sources. This has led to the introduction of stochastic versions of the model (1), where the function is assumed to contain a white noise component. In this case, equation (1) is more properly interpreted as a stochastic differential equation

(2)

where is a Wiener process and and are deterministic functions. Models of this type have been studied in[7],[8],[9],[10],[11] and others. An important feature of these works is the calculation of the internal variability of the system as defined by the variance of the process .

The objective of this paper is to develop a stochastic model to study the nature of the interaction of two isomers, S and R, acting in combination on a single individual. We assume that the combined concentration is an a priori unknown linear combination of the two isomers. We determine this linear combination using a synergy test.

The underlying methodology is as follows. The amounts of isomer S and isomer R in the system at any given time are assumed to follow diffusion processes and of the form (2), driven by the same Wiener process . We further assume that the combined concentration of the chemicals is given by the linear combination of the two processes, for some values of and to be determined numerically. In Section 2, we give formulae for the mean and variance of the process and identify the distribution of as Gaussian. The dataset is described in section 3, followed by the description of the synergy test and the statistical methodologies in the subsequent sections. We present results in section 6 and a simulation study comparison in section 7. In Section 8, we provide the discussion of the findings.

2. Distribution Theorems

We specialize the model (2). Assume that the level of concentration of the isomer in the system follows a linear stochastic differential equation of the form

(3)

Here is the initial amount of isomer, is its rate of absorption, is a diffusion constant, and is a standard Wiener process. The rates , , and the diffusion coefficient are considered to be constants. We assume that the initial concentration is zero. The following results will be needed in Section 6.

Theorem 1.

For each , the random variable has a Gaussian distribution with mean

(4)

and variance

(5)

We now describe the distribution of a linear combination of two processes satisfying equations of the form (3). To this end, consider the linear stochastic differential equations

(6)

with for i =1, 2. Since the chemicals are acting simultaneously, we assume the equations for and are driven by the same noise process . Set

(7)

and

(8)

In the sequel, let denote the process

(9)

Theorem 2.

The process is Gaussian with a mean

(10)

and a variance

(11)

Here and are initial amounts, and are the rates of eliminations, and are the rates of absorptions for two chemicals. and represent the variability coefficients within each process, and and are two constants as used in equation (9). We refer the and the as the main parameters for the model, as the variance parameters, and as the synergy parameters.

3. Description of Data

This study focuses on the concentrations of two stereoisomers, referred to as S and R, in the human bloodstream. Isomers are essentially identical chemical substances that differ only in their stereochemistry. That is, the elemental make-up of the substances is exactly the same, but the 3-dimensional positioning is different. Specifically for S and R isomerism, certain components of the compound are arranged in a different circular pattern about the center. S isomers indicate “left-hand” arrangements while R isomers indicate “right-hand” arrangements. These differences are often compared to the direction of spokes on a wheel as in[12].

When the experiment begins, a blood sample is taken from each of seven volunteers (referred to as subjects). It is determined that the isomers are not present in any of the initial blood samples. That is, the concentration at time zero is 0 ng/mL as in[13]. A mixture containing equal parts of both isomers is applied via eye drops to the subjects at a rate of one drop per minute for three minutes. At the five minute time mark and eleven intervals thereafter, blood samples are taken from each subject and the concentration (in ng/mL) is recorded. We define the concentration of the single dose (either R or S) by the concentration amount of the isomer over eleven time points. We also define the concentration of the combined dose by using the sum of the S and R concentrations over the corresponding eleven time points.

4. A Study for Synergy

Reference[14] proposed a synergy test of two-drug combinations that does not require modeling a response surface. They have used a Model Free Test (MFT) to establish a sufficient condition for synergy at a combination dose. The paper discusses the study of two or more agents in combination of dose ratios and describes synergy as it applies to dose equivalents of the mathematical equation of the form below.

(12)

Here, is the dose (or concentration) of the ith drug given in a combination of n drugs and is the dose (or concentration) of the ith drug given individually that would produce e, the magnitude of the combination concentration. If the value of the left-hand side of the equation (12) is less than 1, the effect is synergistic; if the effect is greater than 1, the effect is antagonistic. The equation (12) defines the theoretical, or more precisely, the e-theoretical line of dose additivity, where lies on an e-isobole as in[14]. Here the authors describe the e-theoretical line as the line connecting the dose points of to. The e-isobole is the curve describing the different dose combinations between and . Depending on the position on the line, one can determine if the combination is synergistic or antagonistic. Their proposed method for establishing a sufficient condition for a synergy test uses an arbitrary line through given by where is a given dose point.

In the fixed ratio design, as in[15] of two drugs, individual compounds are combined together in amounts such that the proportion between them is constant. In other words for two levels of drugs we look at linear combinations of the form

or , which must either be greater than, less than, or equal to 1.

In the expression, is a function of two drugs, , and are the functions of one drug in the absence of the other and r is the correlation coefficient that depends on the ratio of the combination of the two. For the synergy, we use the following tests,

for isomer S, for a specified value of r. A necessary condition for this hypothesis test is that the power for the specified value of r is at least 0.5. Also, since correlation is time dependent and expected to change as concentration combination changes, we establish a bound on r for significant results of the synergy tests for all subjects.

5. Statistical Methodology

The NLIN procedure in SAS uses the mean and variance calculated from the stochastic model without any assumptions about the parametric form of the distribution. The procedure uses the least squares method to fit the curve to the observations and estimate the model parameters. The process requires that the first derivatives of the equation with respect to each model parameter to be estimated. Mean Squared Error (MSE) is calculated as a result of convergence of the NLIN procedure.

The NLMIXED procedure in SAS requires that the mean and the variance expressions for the equations be supplied as the initial input with the assumption that model has a normal distribution with the stated parameters. The procedure uses maximum likelihood estimates of the parameters while fitting the model to the observations. It includes information about the mean and variance for the combined data, and hence we estimate the combination coefficients and of the drug efficacies. We use individual data for later analyses where and are considered to be known constants.

Akaike’s information criterion (AIC) is a useful statistic for statistical model evaluation and has been widely accepted in some areas of statistics, eg. See[16]. It is calculated for each selected model as AIC = (n)ln(SSE/n) + 2k , where k is the number of parameters to be estimated and SSE stands for sum of squared errors. A low value for AIC indicates a better fit as described in[17]. The value of AIC is computed after the convergence of the NLMIXED procedure. The value of AIC is calculated for NLIN procedure from the respective MSE values.

As suggested in[18] we use the Wilcoxon Rank Test for the linear combination of hypotheses, that were described in the synergy test establishes a sufficient condition for rejecting the hypotheses at a .05 significance level. Since the distribution of the time data for each subject is unknown, the nonparametric min test is appropriate and widely used. We perform the power study using the WILCOX.TEST procedure in R for a simulation size of 5000 data sets. For the synergy hypotheses, we test seven pairs of hypotheses against one-sided alternatives, one pair per subject on each isomer type for combined and individual data.

For model (9), using R software we generate 20,000 data sets for each patient using the estimated parameters from both the NLIN and NLMIXED procedures. We supply the initial estimates of the parameters, and use the NLS function to check the convergence of the model parameters to their initial values.

6. Results

The initial dose amounts are considered known and used as constants throughout the analyses. We use the combined data for seven patients to estimate the parameters in equation (9). We use the equations (10) and (11) to input the mean and variance required in the procedure. NLMIXED was necessary to incorporate the subjects as blocking variables. This allowed us to analyze the data in its entirety and estimate the synergy parameters and as given in Table 1. This capability is not available from the NLIN procedure.

However, we use the individual data to estimate mean and variance parameters for seven subjects using both the NLIN and NLMIXED procedures. The results are presented in Tables 2 - 3, using the given coefficients of and from Table 1.

Table 1. Parameter estimates from PROC NLMIXED using the full model, using combined data (standard errors are in parenthesis)

Table 2. Parameter estimates for subjects 1-7 from PROC NLIN using individual data (standard errors are in parenthesis)

Table 3. Parameter estimates for subjects 1 - 7 from PROC NLMIXED using individual data (standard errors are in parenthesis)

Figure 1. Concentration versus time for individual subjects plotted with estimated curves by NLIN (dotted lines) and NLMIXED (solid lines) methods on the observed values

The above two procedures generated parameter estimates with reasonably low standard errors of estimation. An AIC comparison further shed insights on the nature of the analysis done. The AIC values in Table 4 show lower numbers for the NLIN procedure in every subject than the NLMIXED and points towards a better fit of the data by the former procedure.

Using the parameter estimates obtained from the data in previous tables, we fit the estimated equations to the individual data as shown in Figure 1. The fitted curves below indicate an extremely good fit of the model to the data by both the procedures. A careful look at the fitted curves by the NLIN procedure confirms that it follows the data for the individual subjects slightly better than the curves fitted by the NLMIXED procedure.

The results for power study of each are given in Table 5. The lower bound for r is the maximum (i.e. it is equal to 1.00). The upper bound for r is the estimate that generated a power, which we accept (i.e. it is above 0.5). Using these ranges, we perform the synergy tests on the observed data.

In table 5 the sixth subject shows the maximum strength of maintaining the power. The synergy tests for individuals based on the choice of correlation coefficient between .08 and .23 (over all data) gives p-values less than .05 for each subject when performed for both the isomers. The W-score for the Wilcoxon Rank-Sum statistics and the p-values are shown in Tables 6 and 7.

A description of the estimated values for the correlation coefficient is shown below in the isobologram in Figure 2. The theoretical isobole described as in[14] illustrates two hypothetical doses with a generic placement of the synergistic locations. The straight line represents the e-theoretical line of additivity connecting and. In Figure 2, we use the same line setup and place the correlation bounds to show the synergistic power and possible location for the dose combination of the data.

Table 4. AIC comparisons for the individual subjects using two SAS procedures SUBJECT NLMIXED NLIN 1 79.2 35.8 2 93.6 68.5 3 88.5 54.3 4 98.0 71.8 5 94.9 63.0 6 115.8 78.5 7 97.7 72.1

Table 5. Power study for the synergy tests using individual data Subject Result for r Lower Bound Upper Bound 1 0.08 0.23 2 0.17 0.26 3 0.19 0.27 4 0.20 0.29 5 0.12 0.23 6 0.34 0.71 7 0.21 0.28

Table 6. Synergy tests for hypothesis      for subjects 1 – 7 using the isomer R Subject Wilcoxon Rank-Sum Results Lower Bound Upper Bound 1 W-score 121 97 p-value <0.0001 0.009025 2 W-score 121 92 p-value <0.0001 0.01999 3 W-score 118 101 p-value <0.0001 0.003317 4 W-score 121 97 p-value <0.0001 0.009025 5 W-score 119 89 p-value <0.0001 0.03291 6 W-score 121 90 p-value <0.0001 0.02837 7 W-score 121 97 p-value <0.0001 0.009025

Figure 2. A hypothetical isobologram with the estimated bounds for the correlation coefficient and the corresponding power is shown here. Note that the isobole curve represents only a sampling of the possible dose combinations

Table 7. Synergy test for hypothesis      for subjects 1 – 7 using the isomer S Subject Wilcoxon Rank-Sum Results Lower Bound Upper Bound 1 W-score 121 110 p-value <0.0001 0.0006443 2 W-score 119 104 p-value <0.0001 0.002369 3 W-score 121 104 p-value <0.0001 0.002369 4 W-score 113 96 p-value <0.0001 0.009615 5 W-score 121 110 p-value <0.0001 0.0006443 6 W-score 121 114 p-value <0.0001 <0.0001 7 W-score 121 97 p-value <0.0001 0.009025

7. Simulations and Efficiency Comparison for Parameter Estimations

A Bootstrapping Monte Carlo simulation study was conducted for individual subjects to compare the inferential performance of the NLIN and NLMIXED procedures described in section 5. We considered the parameter estimates from the data as the initial starting points. After 20,000 iterations, the estimates converged to assigned criteria and produced the average of each parameter value with the corresponding standard deviation. We used the known estimates of and in the simulation to be consistent with the synergy test.

In most cases in Tables 8 – 9 the estimates are very close to the true parameters in Tables 2 - 3. Differences shown in Tables 10-11 indicate that some of the variance parameters were incorrectly estimated by the data, though the main parameter estimates were quite close. There was only one indication of a large difference that was detected by the simulation. Table 11, the NLMIXED difference table shows a large difference for two main parameter estimates for Subject 1.

Table 8. Simulation results for subjects 1 – 7 using parameters of PROC NLIN procedures (standard errors are in parenthesis)

Table 9. Simulation results for subjects 1 – 7 using parameters of PROC NLMIXED procedures (standard errors are in parenthesis)

Table 10. Difference of true parameter values from NLIN procedure from the estimated values from the simulation results for subjects 1 – 7

Table 11. Difference of true parameter values from NLMIXED procedure from the estimated values from the simulation results for subjects 1 – 7

The NLIN procedure seems to converge faster than the NLMIXED procedure and the parameter estimates are also closer to the actual values obtained from the experimental data. We encountered problems of occurring singular values using the NLMIXED method. The variance estimates seemed not to converge well to the true values for the either procedure. However the bulk of the convergence was smooth, giving precise estimates with very small standard errors.

8. Conclusions

In this paper we have addressed one very important question of drug combination. Pharmacology studies deal with combined data that are interactive in nature. In mathematical terms an interaction is usually represented by multiplication of terms. It is not quite clear for drug interactions whether the combined concentrations are additive, multiplicative, or inhibitive in nature. Combination of isomers in eye doses gives a reason to test for synergy to assert the researchers’ effort to show an additive effect.

The proposed model uses very simple stochastic differential equation techniques to solve for the model mean and variance in closed form. The two different statistical procedures are used here to confirm that our model has the flexibility for use by practitioners.

The proposed model is simple in nature and uses a limited number of dose combinations. We were faced with the challenge of using the existing data. But the data collection can be done based on the recommendations in the paper[14].

While many other methods of drug assessment and model building are available, our model tests for a mathematical synergy instead of a biological synergy. A biological synergy may or may not happen in the experiment; but, if data are collected based on the combined dose levels as suggested by other researchers, a mathematical synergy can be tested as shown in this paper.