A Nonparametric Method for Evaluating Results from Laboratory Antinociceptive Tests
JOHN D. TAULBEE AND GERALD3. KASTING
A method is presented in which Cox’s proportional hazards model, a survival analysis technique, is used to assess the results of hot-plate antinociceptive testing. The method appropriately handles censored data values and variable pretest latency times without making arbitrary assumptions about the distribution of the data. It may be used to characterize and compare dose-response curves or to examine the effect of agent or other treatment variables on the response. The technique is easily implemented using the SAS statistical software package. Due to the similar way in which data is obtained, we believe the method to be applicable to several other laboratory models of pain, including the tail-flick, tailimmersion, and paw-pressure (Randall-Selitto) assays.
Key Words:
Statistical methods; Antinociceptive
tests; Censored data
INTRODUCTION In the hot-plate antinociceptive assay, an animal, usually a rodent, is placed on a heated surface and observed until either a sign of discomfiture is exhibited or a predetermined cutoff time has elapsed. The antinociceptive effect of an agent or treatment is related to the time the animal remains on the pfate before exhibiting the endpoint. The data analysis for this type of experiment has been discussed recently by Kitchen and Crowder (1985). One major concern is the correct treatment of observations that are censored by the cutoff latency time. Other issues of interest are the form of the dose-response curve and the appropriate use of pretreatment latency times in the analysis. Kitchen and Crowder dealt with the first two of these issues by constructing a parametric model for the data based on the Weibull distribution and a linear log (response time) versus dose curve, The model could be solved iteratively to obtain maximum likelihood estimates of the Weibull parameters and, consequently, the dose required to produce a given mean latency time, such as the ADS0 value (the dose required to produce 50% of the maximum antinociceptive response). They were able to show that the model produced satisfactory fits to hot-plate data generated in mice with a variety of antinocic~ptive agents. We have also approached the problem from a survival analysis perspective, but have used a nonparametric regression model (Cox, 1972) rather than making ar-
FromThe Procter& GambleCompany,Miami Valley Laboratories, Cincinnati, Ohio. Address reprint requests to: Dr. John D. Taulbee, The Procter & Gamble Laboratory, Laboratories, P.O. Box 395707, Cincinnati, OH 45239-8707. Received December 9,1986; revised and accepted April 17,1988.
Miami Valley
197 Journal of Pharmacological 0 7988 Elsevier
Science
Methods Publishing
20,197-206 Co.,
Inc.,
0160-5402/1988BO3.50
KWJ) 655 Avenue
of the Americas,
New
York,
NY tMFl0
198
J. D. Taulbee and G. 8. Kasting
bitrary assumptions about the shape of either the survival curve or the dose-response curve. The method is compatible with Kitchen and Crowder’s analysis, and it has the further advantage of allowing pair-wise comparisons of treatment variables other than dose and of incorporating the pretreatment latency time as a covariate in the regression, if significant. The model and parameters required for these methods will be discussed first, with their interpretations described. Examples will then be given, followed by a discussion of results. Model
and Parameters
Our model is a nonparametric-time-until-event model. The rate of occurrence of the event of interest, e.g., hind-paw lick, is modeled as a function of treatment and predose latency time (the time until occurrence of the event prior to dosing). The rate of occurrence is just the inverse of average time until occurrence in the case of exponentially distributed data. For statistical reasons, it is more convenient to model the rate of occurrence than the time until occurrence. A higher rate of occurrence will be associated with shorter times until occurrence, while a lower rate of occurrence will be associated with longer times. We will denote the rate of occurrence of our event of interest as h(t), also known as the hazard function or hazard rate. The underlying models, derivations, and theory are found in the Appendix. In the hot-plate assay for dose-response, the dose administered, D, is an experimentally controlled variable, and the pretest latency time, TP, is an observed variable. Both D and TP may be important in predicting the actual postdose latency time, T. For our model, we would have that the rate of occurrence is h(t)
= h,(t)
exp[B,D
+ B,T,]
(1)
where h,(t) is an “underlying” rate that does not need to be estimated with the Cox (Cox, 1972) methodology. We would be interested in tests about B, and BZ. The interpretation for B, and BZ values is given below: B, < 0 B, = 0 B, > 0
Increasing dose results in lower hazard rate, i.e., longer time until response. Dose is not related to time until response, i.e., agent dosed does not demonstrate effectiveness. Increasing dose results in higher hazard rate, i.e., shorter time until response.
Analogous interpretation is true for &. So that we can interpret the magnitude of B,, we will take Dd, the “doubling dose,” as the dose necessary to reduce the rate of occurrence of the event of interest by a factor of 2 from that which occurs with zero dose. Holding predose latency
Evaluating Results of Antinociceptive
time
Tests
constant at Tp, and using Eq. (I), we find Dd =
-In
2
(2)
B
1
(See Appendix
for details.)
where 5, is estimated using Cox’s (Cox, 7972) method. The choice of reduction by a factor of 2 is not essential. We could pick any factor, and substitute that factor for 2 in Eq. (2). This method allows for easy comparison of treatments as well as dose-response. To compare treatments A and 5, we create a dummy variable, I, which indicates which treatment is being discussed. Let I = 0 for treatment A, I = 1 for treatment 5. Then for our model, the rate of occurrence is h(t) For interpretation, 5, < 0 5 = 0 5, > 0
= h,(t)
exp[&/
+ 5~7-J
(3)
we have
Treatment 5 has a lower hazard rate (longer time until response) than treatment A. Treatments A and 5 are not different in their effects on response time. Treatment B has a higher hazard rate (shorter time until response) than treatment A.
To summarize the relationship of A to 5, we define the effectiveness ratio of A to B, ERAs, as the ratio of the hazard rate of B to that of A. Using Eq. (3) at a fixed predose latency time, Tp, ERAS
(See Appendix
=
exp[&I
(4)
for details.)
If we use the interpretation of 51 above, for B, < 0, ERAs < 1, and A is less effective than B; for B1 = 0, ERAB = 1, and A and B are equally effective; and for B, > 0, ERAB> 1, and A is more effective than B. Note that ERAs = l/ERsA. Using Cox’s model, we can test whether predose latency time, Tp, is predictive for postdose latency time. If Tp is predictive, it should be included in the model for its reduction of error variability. If T, is not predictive, it can easily be dropped from the model. EXPERIMENTAL Male Sprague-Dawley rats, 100-120 g, are placed one at a time on a copper surface heated to 55°C. A Plexiglas cylinder, 20 cm in diameter by 30.5 cm high and open on both ends, is used to restrict the animals to a defined area of the surface. The animal is observed until either a jump or hind-paw lick endpoint is reached. A jumping response is reached when the animal either jumps to the top of the cylinder or
199
200
j. D. Taulbee and C. B. Kasting
when all four paws leave the surface in an attempt to jump. A hind-paw lick endpoint is reached when the animal’s mouth comes into contact with either hind paw through licking or biting. The latency time is recorded as the elapsed time from initial contact with the hot-plate surface to occurrence of the endpoint, or as 60 set if this behavior is not seen. The animal is removed from the plate immediately after exhibiting the endpoint or at 60 set, whichever occurs first. The antinociceptive agents discussed in the dose-response and treatment comparison examples were experimental agents, administered orally. Hot-plate latency times were recorded immediately prior to dosing (to obtain Tp) and at 1.5,3, 5, and 24 hr postdose. Dose-Response
Example
In this study, a compound was dosed at levels of 0, 75, 150, and 300 mg/kg with eight animals at each dose level. Predose latency times and 5-hr postdose latency times are available. The data are found in Table 1. Using Cox’s regression in PROC PHGLM from the SAS statistical package (Harrell, 1983), and the model of Eq. (I), we find these results:
STANDARD
REGRESSION COEFFICIENT
p-VALUEFORTESTINC WHETHER
ERROR
ESTIMATE
COEFFICIENT = 0
Bl
-0.0068
0.0027
0.013
&
-0.1861
0.1143
0.103
&
-0.0058
0.0026
0.022
(without f32)
We find 8, < 0, so increasing dose statistically significantly reduces the hazard rate, i.e., increases postdose latency time. Also, although it is not significant at the usual 0.05 level, BZ < 0, which suggests that longer predose latency times are associated with longer postdose latency times. The doubling dose estimated from B, in the model using predose latency times is, from Eq. (2), Dd = -In
2/( -0.0068
‘_’ 0.0027)
f67 = 102
mg/kg -29
where the errors reflect the standard error of B,. This doubling dose can be used to assess the relative potency of the compound used here, say compound A, with the Dd value for another compound, B, tested in this assay. Suppose the Dd value was 204 mg/kg for compound B. Then the relative potency of A to B would be 204/ 102 = 2. That is, it takes twice as much of compound B to produce the same reduction in hazard rates as produced by a given dose of compound A.
Evaluating Results of Antinociceptive
Tests
TABLE 1. Hot-plate latency Times Predose and 5 Hr Postdose for 4 Doses of an Experimental Antinociceptive Agent” DOSE LEVEL,mg/kg
LATENCY
-
TIME
300
150
75
0 PRE
POST
PRE
POST
PRE
POST
10.0
PRE
_~
POST 60.0*
60.0*
4.8
22.9
a.1
35.7
10.6
9.5
6.7
7.7
20.6
8.4
33.7
5.9
5.4
6.8
5.5
21.3
7.8
6.8
6.6
6.7
8.3
7.4
9.8
60.0'
6.6
37.6
6.9
60.0b
10.3
4.6
7.2
7.8
60.0b
6.8
60.0b
8.0
60.0b
6.7
9.9
9.0
19.5
5.6
60.0b
8.0
60.0'
9.6
5.6
11.5
60.0'
11.0
44.3
6.0
21.7
12.0
10.9
5.6
22.6
7.3
36.5
6.0
60.0b
a Eight animals per group.
’ Censored observation (animal had not yet responded,
was removed from plate).
To test whether or not the doubling doses for two compounds are different, let BA be the value of & for compound A and SEA its standard error, likewise BB and SEB for compound B. Then
= = v’(SE,)*
4 (SE,)’
has a standard normal distribution under the null hypothesis of no difference between compounds, and the usual test of treatment differences for normally distributed variables with known variances applies. As noted earlier, the doubling dose is a convenient way of expressing the magnitude of the parameter, B,, which is the fundamental quantity estimated by Cox’s model. One is in no way constrained to calculating a doubling dose; one can just as easily calculate a tripling dose or an n-tupling dose simply by replacing In 2 in Eq. (2) with In 3 or In n. Comparative
of Treatments
Example
In this study, two treatments, A and B, were given with eight animals per treatment. Latency times were measured predose and 3 hr postdose, The results are in Table 2. Using the model of Eq. (31, we find these results:
REGRESSION COEFFICIENT 51 52 51 (without 5,)
ESTIMATE -1.44 -0.63 -0.78
STANDARD ERROR 0.70 0.34 0.55
~-VALUEFORTESTING WHETHERCOEFFICIENT= 0 0.038 0.061 0.154
201
202
J. D. Taulbee and G. B. Kasting TABLE 2 Hot-plate Latency Times Predase Antinociceptive Agents A and B
and 3 Hr Postdose for Experimental
TREATMENT
_______^“. A __~I”Latency time
PREDOSE
._.
4.2 6.4 6.0 7.3 5.4 8.1 7.8 7.2
POWJOSE
-.... PREDOSE
4.2 6.6 5.8 60.0” 5.2 14.3 16.4 71.9
5.0 6.9 7.2 5.7 5.4 7.2 7.0 7.3
5 POS~DOSE 14.9 30.0 24.0 13.9 60.0a 13.4 22.5 21.8
a Censored observation (animal had not yet responded, was removed from plate).
We see that treatment 5 is associated with significantly longer postdose times on the plate than is treatment A, after adjusting for predose latency. The effectiveness ratio of A to B is, from Eq. (4) and 5, with 52 in the model, ER.+M= expt-I.441
= 0.24
The hazard rate for treatment 8 is about one fourth that for treatment A. Thus, B is a significantly more effective treatment than A (p = 0.038). If we use average postdose latency times as our endpoint, the result for treatment A is 15.5 -+ 18.5 set and for B is 25.1 + 15.2 set, which are not significantly different by the t-test (p = 0.28). If the predose values are subtracted prior to averaging the postdose values (to obtain change in latency values), one obtains 9.0 i: 18.0 set for A and 18.6 + 15.5 set for 5, which again are not significantly different by t-test (p = 0.27). Normaii~ing the latency values by the range of the test according to Percent Antinociception
= 100 x (Post -
Pre)/(GO -
Pre)
(5)
makes little difference: The percent antinociception values are 17.1 -r- 34.1 for A and 34.6 t 28.3 for B, not significantly different (t = 1.12, p = 0.28). Comparison
of Opioid
with Other Methods Analgesics
for Relative Potency Estimation
The EDSo is a commonly used statistic for summarizing antinociceptive effectiveness. In one setting, it is estimated to be the dose that would produce a significant antinociceptive effect in half of the animals tested at that dose. Significant antinociceptive effect is often taken to be a response time at least three standard deviations above the pretest response time. In another setting, the EDS0 is that dose estimated to produce 50% antinociception where the percent antinociception for each animal is calculated as in Eq. (5). In both cases, the Litchfield-Wilcoxon method is used to estimate EDso.
Evaluating Results of Antinociceptive
Tests
We have data from experinients on three opioid analgesics-codeine phosphate, oxycodone, and morphine sulfate-from several groups of animals. The /Id and both Ei& values are given in Table 3. Note that Dd cannot be directly compared with ED5O, nor is it technically correct to directly compare E5Ovalues from the two methods. However, relative potencies from the Dd and EDSOvalues can be compared for consistency. Table 4 contains relative potencies of morphine to codeine and oxycodone based on ratios of Dd values or EDso values. It can be seen that the relative potencies from the Dd values are in excellent agreement with those from the ED5Ovalues. This demonstrates that the methods presented here produce results consistent with those of methods currently used. DISCUSSION
The method discussed was developed during the course of a development program with a new class of experimental antinociceptive agents known as vanilloids (Berman et al., 1986). The objective of the program and the nature of the compounds being tested led to the use of the rat hot-plate assay as the primary antinociceptive screen. Potency estimation using a quantile grading system followed by the probit analysis (method of Litchfield and Wilcoxon [1949]) gave satisfactory results for morTABLE 3 EDso Estimated by Pretest +3SD Method and by % Antinociception the Dd Estimates for Three Opioid Analgesics”
COMPOUND Codeine phosphate Codeine phosphate Meanb Oxycodone Morphine sulfate Morphine sulfate Meanb
ROUTEOF ADMINISTRATION
No. OF RATS
p.0. p.0.
40 29
p.0. S.C. S.C.
32 39 36
&i
PRETEST + 3SD METHOD
112 130 120 7 1.7 2.1 1.9
216 329 267 20 4.0 5.1 4.5
Method and
PERCENT ANTINOCICEPTION METHOD 644 472 551 30 6.4 6.2 6.3
a All doses are expressed as m&kg. The doses of codeine and morphine refer to the salt. The morphine responses are 30 min postdose, and the codeine and oxycodone responses are 60 min postdose. b Geometric mean.
TABLE 4
Relative Potencies of Opioid Analgesics by Three Methods RELATIVE POTENCYOF MORPHINESULFATE TO COMPOUND USING:
ED50 COMPOUND Codeine phosphate Oxycodone
ROUTE p.0. p.0.
Dd
(PRETEST + 3SD)
65 3.9
59 4.4
EDso (% ANTINOCICEPTION) 88 4.7
203
204
J. D. Taulbee and G. B. Kasting
phine, but it failed to give reproducible EDso values for the experimental agents. Furthermore, pair-wise comparisons using the t-test were found to be unreliable due to censoring and to a very non-Gaussian distribution of response time. While these previously used statistical methods do enable valid distinction between antinociceptive agents with larger differences in potency, they lack statistical power for determining smaller, but possibly important differences. Thus, these statistical methods presented a problem in selecting compounds for further testing. The survival analysis method presented here provides a marked improvement in our ability to discriminate between antinociceptive agents that are closer to each other in potency. In addition, the proposed statistical method allows for detection of antinociceptive activity in compounds that might have gone undetected by the previously used statistical methods. Thus, while quantile grading or other standard methods of biological data analysis may work better for other classes of analgestics than they did in our case, the method presented here has the potential for increasing the sensitivity of the analysis for these agents as well, allowing high-quality decisions to be made on the basis of fewer animals tested. It appears to be general enough to apply to other assays in which a stimulus is applied up to a predetermined threshold or limit, such as the tail-flick and paw-pressure tests. In addition, the method presented here produces results similar to other methods as evidenced by the comparability of relative potencies from the various methods. CONCLUSIONS The method presented here has several advantages: like that of Kitchen and Crowder, it is a survival analysis method and so takes into account the different nature of censored observations, which is not done with simple averaging or other methods. In addition, it allows for the inclusion of predose latency time, which is, in some instances, an important source of variability that should be taken into account. It can be implemented using the SAS statistical package. Finally, it is a nonparametric approach that requires no assumptions about the distributional shape of response times. REFERENCES Kitchen I, Crowder M (1985) Assessment of the hotplate antinociceptive test in mice: a new method for the statistical treatment of graded data. / fharmacol Methods, 13:1-7. Cox DR (1972) Regression models and life tables (with discussion). / RoyalStat Sot 8,34:X37-220. Harrell
F (1983) PROC PHGLM.
In .X/C/
Supple-
mental Library User’s Guide, 1983 Edition. Cary, N.C.: SAS Institute, Inc. Berman EF, Schwen RJ, LaHann TR, Gardner JH, Janusz JM, Buckwalter BL, Brand LM (1986) NE19550, a novel non-narcotic agent. Fed Proc, 45:665. Litchfield JT and Wilcoxon F (1949) 1 Pharmacol Exp Ther, 96:99.
APPENDIX Theory and Underlying
Models
Survival analysis methods allow for the cutoff time to be acknowledged as a censored, rather than actually observed, measurement. The survival probability at a
Evaluating Resufts of Antinocicepti~e
Tests
given time t, expressed as S(t), is just the proportion of subjects that have yet to respond at time t. Typically, the functional form of S(t) is such that the distribution is characterized by h(t)
=
-d/n[S(t)] &
the so-called hazard function or hazard rate. The instantaneous probability of response in a small time interval [t, t + At1 is equal to h(t)At. For the Weibull distribution discussed by Kitchen and Crowder (I), h(t) A special case of this distribution
= hyPY-’
is where y = 1, so that
h(t)
= A, a constant.
This is the exponential distribution, and here the average time until response is just X-l. In general, the average time until response increases as the hazard rate decreases. Cox (2) introduced a nonparametric method for assessing the effects of experimentally controlled or observed variabies on time until event. For his model h(t)
= h,(t)
exp[&Z,
+ L&Z2 + +..I
where h,(t) is some underlying hazard rate that does not need to be estimated (hence the method is nonparametric) and where Zt , Z2, . . . are variables of interest and B,, B2, . . . are, respectively, the regression coefficients for the variables of interest. In the hot-plate assay for dose-response as described in the text, we would have, using Cox’s model, h(t)
= h,(t)
exp[B,D
+ B2TPl
(1)
where we would be interested in tests about 3, and &. So that we can interpret the magnitude of B,, we will take &, the “doubling dose,” as the dose necessary to reduce the hazard rate by a factor of 2 from that which occurs with zero dose. t-iolding predose latency time constant at TP, and using Eq. (I), I/2 =
MU
expl&Dd
h,(t) exp&NB
+ &T,l + B2TPl
exi,[&Ddl Solving for D4,
where 5, is estimated using Cox’s method. The choice of reduction by a factor of 2 is not essential. We could pick any factor, and substitute that factor for 2 in Eq. (2). The ratio of the Dd for agent A to Dd for agent B provides an estimate of the
205
206
J. D. Taulbee and C. B. Kasting
relative potency of B to A. This method allows for easy comparison of treatments as well as dose-response. To compare treatments A and B, we create a dummy variable, I, which indicates which treatment is being discussed. Let I = 0 for treatment A, / = 1 for treatment B. Then for Cox’s model, h(t) For interpretation,
= h,(t)
we have the effectiveness ER/,s = ERAB
+ B,T,l
exp[B,/
(3)
ratio of A to B,
h,(t)
exp[&U)
+ B2T,l
h,(t)
exp[&(O)
+ B2T,J
(4)
exp[&l
=
From Eq. (I), the ratio of the hazard rate at dose D to that at dose zero at a given predose latency, TP, is
For the exponential distributions,
h(Dose
= D)
h(Dose
= 0)
= exp[B.,Dl
(which is a member of the Weibull
TO -= TD
h(Dose
= D)
h(Dose
= 0)
(5) family) and other similar
(6)
where To is the average latency time for dose 0, and To is the average latency time for dose D. So, fromEqs. (5) and (6), To = To exp[-BID1 or, taking logarithms, In To = In To Eq* (7) is equivalent
B,D
to Eq. (2)(i) of Kitchen and Crowder
(7)
(1):
E(ln t ) d) = a + pd where d is dose. (Note that B, will be negative for an effective compound, so that Eq. (7) implies a positive correlation between latency time and dose.) Thus, in the special case of the exponential distribution, the present model gives the form of the dose-response curve assumed by Kitchen and Crowder.