Model the relationship between prion inoculum titer and incubation time

inoculum-incubationtime.r

# Eric Minikel
# 2014-01-12
# explore relationship between inoculated dose and incubation time
# this is all the code from this post: 
# http://www.cureffi.org/2014/01/12/prion-kinetic-models-relationship-inoculum-titer-incubation-time/

logdose = runif(n=100,min=0,max=10) # simulate 100 different experiments with diff titers
days = 10^( (logdose-26.66)/(-12.99) ) + 40 + rnorm(n=100,m=0,s=2) # apply Prusiner's model
plot(logdose,days,xlab="log10(LD50 units)",ylab="days to illness",pch=19,
     main="incubation time interval bioassay\nPrusiner 1982b - Sc237 in hamsters",
     sub="simulated data")

onset_dpi_Sc237 = function(titer) {
    dpi_min = 40
    return ( dpi_min + 10**((1/12.99)*(26.66-log10(titer))) )
}

# RML based on Bueler 1993 Table 1
source_titer = (1e7)*10 
onset_dpi = c(137,151,157,160,183,201) # Table 1 "inoculum" column, "disease" row
dilutions = 10^-(1:6)
titers = source_titer * dilutions
plot(titers,onset_dpi,log='x',pch=19,cex.main=.8,
     main='Inoculated titer vs. incubation time in RML-infected wt mice\nBueler 1993 Table 1')

logtiter = log10(titers)
m = lm(logtiter ~ onset_dpi)
summary(m)
#Coefficients:
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept) 17.34296    1.73838   9.977 0.000567 ***
#onset_dpi   -0.07791    0.01046  -7.449 0.001735 ** 
#Multiple R-squared:  0.9328,    Adjusted R-squared:  0.9159 

m = lm(logtiter ~ log10(onset_dpi))
summary(m)

m = lm(logtiter ~ log10(onset_dpi - 40))
summary(m) # adj R^2 = .9404


# so we can indeed do better with a log-log model
# can we do even better with some constant

find_best_adjustment = function(onset_dpi, titers) {
    # try every possible "adjustment" up to the lowest incubation time observed
    days_to_test = 1:(min(onset_dpi)-1)
    rsq_vals = numeric(length(days_to_test))
    for (dpi_adjustment in days_to_test) { 
        log_adjusted_dpi = log10(onset_dpi - dpi_adjustment)
        m = lm(logtiter ~ log_adjusted_dpi)
        rsq_vals[dpi_adjustment] = summary(m)$adj.r.squared
    }
    # figure out which model had the highest R^2
    return ( which(rsq_vals == max(rsq_vals)) )
}

find_best_adjustment(onset_dpi, titers)
# 100

dpi_adjustment = 100
log_adjusted_dpi = log10(onset_dpi - dpi_adjustment)
m = lm(logtiter ~ log_adjusted_dpi)
summary(m)
#Coefficients:
#                 Estimate Std. Error t value Pr(>|t|)    
#(Intercept)        25.670      2.185   11.75 0.000300 ***
#log_adjusted_dpi  -11.834      1.218   -9.72 0.000627 ***
#Multiple R-squared:  0.9594,    Adjusted R-squared:  0.9492 

# model: estimated_titer = 10^(25.670 - 11.834*log10(onset_dpi - 100))
# rearranging, predicted_onset = 10^((log10(titer) - 25.670) / (-11.834)) + 100

predicted_onset = function(titer) {
    return ( 10^((log10(titer) - 25.670) / (-11.834)) + 100 )
}
plot(titers,onset_dpi,log='x',pch=19,cex.main=.8,
     main='Incubation time bioassay model for RML in mice\nBased on Bueler 1993 Table 1 data',
     sub='model: predicted_onset = 10^((log10(titer) - 25.670) / (-11.834)) + 100')
points(titers, predicted_onset(titers), type='l', col='violet')

# the log-log relationship is irreconcilable with the log-linear relationship predicted by
# the exponential model


tga20_onset = c(69,69,87,111) # Fischer 1996 Table 2
source_titer = (10^7)*10
dilutions = 10^c(-1,-2,-4,-6)
titers = dilutions*source_titer
logtiter = log10(titers)
m = lm(logtiter ~ tga20_onset)
summary(m) # .9245

find_best_adjustment(tga20_onset, titers)
# 46

m = lm(logtiter ~ log10(tga20_onset - 46))
summary(m) # .9492