Using R package drc for 4PL curve fitting

Using the 4 parameter logistic model (4PL) to study a dose-response curve and calculate the IC50/ED50/EC50 dose is often cumbersome in Excel or requires using commercial software.

Here shows how to use an R package drc to analyze the dose-response curve of the drug CB83 in a clonogenic assay.

The formula of the 4PL is given as \[y=\frac{a-d}{1+(x/c)^b}+d\] where y is the response and x is the concentration. The lower asymptote is a, the bottom of the curve or lower plateau (commonly referred to as the min) and the upper asymptote is d, the top of the curve or upper plateau (commonly referred to as the max). The steepness of the linear portion of the curve is described by the slope factor (also called the Hill slope), b . The parameter c is the concentration corresponding to the response midway between a and d. .

Let’s look at the raw data.

dose	surviving.fraction	cell.line
0.0	1.0000000	SCC-61
0.5	1.1126761	SCC-61
1.0	1.1971831	SCC-61
1.5	0.9859155	SCC-61
2.0	1.0704225	SCC-61
2.5	0.9436620	SCC-61
5.0	0.1760563	SCC-61
7.5	0.0028169	SCC-61
0.0	1.0000000	rSCC-61
0.5	0.3636364	rSCC-61
1.0	0.8961039	rSCC-61
1.5	0.4350649	rSCC-61
2.0	0.5930736	rSCC-61
2.5	0.1515152	rSCC-61
5.0	0.0025974	rSCC-61
7.5	0.0000000	rSCC-61

There are two options we can choose to fit the curve:

force the min to 0 and max to 1, in this case

#fit in a 4PL model with set lower and upper limits
model1 <- drm(surviving.fraction ~ dose,cell.line,
             fct = LL.4(names = c("Slope", "Lower", "Upper", "ED50")),
             data = data,
             lowerl = c(NA,0,NA,NA), #set the lower limit
             upperl = c(NA,NA,1,NA)) #set the upper limit

force the min and max to be the same between the two categorical variables (SCC-61 and rSCC-61), in this case

#fit in a 4PL model with common lower and upper limits
model2 <- drm(surviving.fraction ~ dose,cell.line,
             fct = LL.4(names = c("Slope", "Lower", "Upper", "ED50")),
             data = data,
             pmodels = data.frame(cell.line, 1, 1, cell.line)) #common lower and upper limits

But here we don’t want to choose neither of those options. So we run another model.

model3<- drm(surviving.fraction ~ dose,cell.line,
             fct = LL.4(names = c("Slope", "Lower", "Upper", "ED50")),
             data = data)

Quickly check the model summary and plot

## 
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
## 
## Parameter estimates:
## 
##                 Estimate Std. Error t-value   p-value    
## Slope:SCC-61   5.2843125  4.0852318  1.2935  0.231929    
## Slope:rSCC-61  5.1255797  7.2096275  0.7109  0.497325    
## Lower:SCC-61   0.0015613  0.2659521  0.0059  0.995460    
## Lower:rSCC-61 -0.0011548  0.1497806 -0.0077  0.994037    
## Upper:SCC-61   1.0829242  0.1053360 10.2807 6.901e-06 ***
## Upper:rSCC-61  0.7229636  0.1489842  4.8526  0.001268 ** 
## ED50:SCC-61    3.6345600  1.1125919  3.2668  0.011412 *  
## ED50:rSCC-61   2.1719768  0.4822210  4.5041  0.001991 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error:
## 
##  0.2060725 (8 degrees of freedom)

So the ED50 of CB83 is 3.63456 \(\mu\)M in SCC-61 and 2.1719768 \(\mu\)M in rSCC-61.

We can also obtain other effective doses using ED()

ED(model3, c(10, 90), interval = "delta")

## 
## Estimated effective doses
## 
##              Estimate Std. Error    Lower    Upper
## e:rSCC-61:10  1.41476    1.01821 -0.93324  3.76276
## e:rSCC-61:90  3.33448    1.84987 -0.93133  7.60028
## e:SCC-61:10   2.39812    0.92984  0.25392  4.54233
## e:SCC-61:90   5.50848    2.71951 -0.76272 11.77969

The following is to produce more sophisticated figures with ggplot2

# make new dose levels to calculate the predicted values based on the model
newdata <- data.frame(dose = rep(exp(seq(log(0.00001), log(10), length = 100)),2), cell.line = rep(c("SCC-61","rSCC-61"),each = 100))

# predictions and confidence intervals based on model3
pm <- predict(model3, newdata=newdata, interval = "confidence")

# append predictions to newdata
newdata$prediction <- pm[,1]
newdata$prediction.lower <- pm[,2]
newdata$prediction.upper <- pm[,3]

# plot with confidence intervals
ggplot(data, aes(x = dose, y = surviving.fraction,shape = cell.line)) +
  geom_point() +
  geom_ribbon(data = newdata, aes(x = dose, y = prediction,
                                  ymin = prediction.lower,
                                  ymax = prediction.upper), 
                                  alpha = 0.2) +
  geom_line(data = newdata, aes(x = dose, y = prediction)) +
  xlab(expression(paste("Dose (",mu,"M)"))) +
  ylab("Surviving fraction") +
  theme(legend.title=element_blank())

Normally we don’t plot that, so

p <- ggplot(data, aes(x = dose, y = surviving.fraction, shape = cell.line)) +
  geom_point(size = I(3)) +
  geom_line(data = newdata, aes(x = dose, y = prediction)) +
  xlab(expression(paste("Dose (",mu,"M)"))) +
  ylab("Surviving fraction")

q<- p + theme(legend.position = c(1,1), # place legend in the upper right corner
    legend.justification = c(1,1),
    legend.title = element_blank(),
    legend.background = element_rect(color = "black"), # add a box around the legend
    # remove panel background and grid
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(colour = "black")) # add back axis lines
q

If we want to use a logarithmic scale for the dose: (note the warning message about taking the log of 0)

q+scale_x_log10()

## Warning: Transformation introduced infinite values in continuous x-axis

Using R package drc for 4PL curve fitting

XC

9/14/2020