Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Model the relationship between prion inoculum titer and incubation time
# 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.