The RT-qPCR process leads to an exponential growth phase of the viral load, such that after some number of cycles \(n\), the viral load \(R_n\) will be a function of the efficiency of the reaction \(E\), the dilution factor and the starting viral load in the sample \(R_0\):
\[ R_n = \frac{R_0}{D}(1 + E)^n \]
Let the probability of a positive case being detected be 0.98 for an undiluted sample, and 0.8 for a 10x diluted sample. Allow up to \(C_q = 35\) cycles to reach some fluourescence detection threshold $R_{C_q}. Let \(R_0\) be distributed in the population such that
\[ R_0 \sim 10^{U(0, 3)} \]
and assume that efficiency is some distribution between 0.65 and 0.9 e.g.
\[ E ~ \sim Beta(\alpha, \beta) \times (0.9 - 0.65) + 0.65 \]
We need to find some distributions of \(R_0\) and \(E\) that will result in a sensitivity of 0.98 at \(D=1\) and 0.8 at \(D=0.8\)
Simulate some values of \(E\) and \(R_0\)
set.seed(12345)
library(scales)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
n <- 10000
R0 <- 10^runif(n, 0, 3)
E <- rbeta(n, 2, 5) %>% rescale(., to=c(0.65, 0.9))
hist(R0)
hist(E)
Calculate the \(R_{C_q}\) values for different dilutions
R35_1 <- R0 * (1 + E)^35
R35_10 <- R0/10 * (1 + E)^35
R35_30 <- R0/30 * (1 + E)^35
Are the sensitivities at the different dilutions concordant?
(quantile(log(R35_1), 0.02) - quantile(log(R35_10), 0.2))^2
## 2%
## 0.1844758
Quite similar. What sort of dilution curve does this give us?
D <- seq(1,30)
Rcq <- mean(c(
quantile(log(R35_1), 0.02),
quantile(log(R35_10), 0.2)
))
dat <- tibble(
D=D,
sensitivity=sapply(D, function(x)
{
sum(log(R0/x * (1 + E)^35) > Rcq)/n
})
)
plot(sensitivity ~ D, dat)
Use optimisation to find parameters that will give concordant sensitivities at D = 1 and 10
fn <- function(param, R0)
{
n <- length(R0)
E <- rbeta(n, param[1], param[2]) * (0.9-0.65) + 0.65
R35_1 <- R0 * (1 + E)^35
R35_10 <- R0/10 * (1 + E)^35
q1 <- quantile(log(R35_1), 0.02, na.rm=TRUE)
q10 <- quantile(log(R35_10), 0.2, na.rm=TRUE)
return((q1 - q10)^2)
}
fn(c(2,5), R0)
## 2%
## 0.3064208
n <- 100000
R0 <- 10^runif(n, 0, 3)
o <- optim(par = c(2, 5), fn=fn, R0=R0)
## Warning in rbeta(n, param[1], param[2]): NAs produced
o
## $par
## [1] 2.5880577 0.9795588
##
## $value
## [1] 4.684844e-08
##
## $counts
## function gradient
## 299 NA
##
## $convergence
## [1] 10
##
## $message
## NULL
Now that we have parameter values for the distribution of \(E\) we can generate the dilution curve
D <- seq(1,100)
E <- rbeta(n, o$par[1], o$par[2]) * (0.9-0.65) + 0.65
Rcq <- quantile(log(R0 * (1 + E)^35), 0.02)
dat <- tibble(
D=D,
sensitivity=sapply(D, function(x)
{
sum(log(R0/x * (1 + E)^35) > Rcq)/n
})
)
plot(sensitivity ~ D, dat)
The above leads to unrealistic Cq distributions. Change approach - generate viral load by simulating Cq values first. Need to find distributions of Cq and E that are reasonable, and give appropriate sensitivities at dilutions of 1x and 10x
See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7454307/ Figure 3 for distribution of Cq values
set.seed(12345)
n <- 500000
rcq <- 10
cqthresh <- 35
e_max <- 1.4
e_min <- 0.9
e_alpha <- 2
e_beta <- 4
sensitivity_1 <- 0.983
tpr1 <- 0.983
tpr10 <- 0.8
est_rcq <- function(r0, cqmax, E)
{
r0 * (1 + E)^cqmax
}
est_r0 <- function(cq, E, rcq)
{
rcq / (1+E)^cq
}
est_cq <- function(r0, E, rcq)
{
log(rcq/r0) / log(1+E)
}
get_cq <- function(n, uci, lci, cq_alpha=1.1, cq_beta=1.1)
{
rbeta(n, cq_alpha, cq_beta) * (uci-lci) + lci
}
get_E <- function(n, uci, lci, perc=0.99)
{
mu <- mean(c(uci, lci))
s <- (mu - lci) / qnorm(perc)
rnorm(100000, mu, s) %>% pmax(., 0) %>% pmin(., 1)
}
fn <- function(param, n, rcq, e_max, e_min, cqthresh, tpr1=tpr1, tpr10=tpr10)
{
cqd <- get_cq(n, param[1], param[2], param[3], param[4])
Ed <- rbeta(n, param[5], param[6]) * (e_max-e_min) + e_min
r0 <- est_r0(cqd, Ed, rcq)
r0_10 <- r0 / 10
cqd_10 <- est_cq(r0_10, Ed, rcq)
s1 <- sum(cqd < cqthresh)/n
s10 <- sum(cqd_10 < cqthresh)/n
return((s1 - tpr1)^2 + (s10 - tpr10)^2)
}
(efficiency_params <- optim(par = c(35, 15, 1.1, 1.1, 1, 1), fn=fn, n=n, rcq=rcq, e_max=e_max, e_min=e_min, cqthresh=cqthresh, tpr1=tpr1, tpr10=tpr10, control=list(maxit=1000)))
## Warning in rbeta(n, cq_alpha, cq_beta): NAs produced
## Warning in rbeta(n, cq_alpha, cq_beta): NAs produced
## Warning in rbeta(n, cq_alpha, cq_beta): NAs produced
## Warning in rbeta(n, param[5], param[6]): NAs produced
## Warning in rbeta(n, cq_alpha, cq_beta): NAs produced
## Warning in rbeta(n, param[5], param[6]): NAs produced
## Warning in rbeta(n, param[5], param[6]): NAs produced
## $par
## [1] 35.457963 15.822604 1.647962 1.234565 1.747571 1.976305
##
## $value
## [1] 1.0132e-08
##
## $counts
## function gradient
## 767 NA
##
## $convergence
## [1] 10
##
## $message
## NULL
cqd <- get_cq(n, efficiency_params$par[1], efficiency_params$par[2], efficiency_params$par[3], efficiency_params$par[4])
Ed <- rbeta(n, efficiency_params$par[5], efficiency_params$par[6]) * (e_max-e_min) + e_min
r0 <- est_r0(cqd, Ed, rcq)
r0_10 <- r0 / 10
cqd_10 <- est_cq(r0_10, Ed, rcq)
s1 <- sum(cqd < cqthresh)/n
s10 <- sum(cqd_10 < cqthresh)/n
sum(cqd < cqthresh) / n
## [1] 0.983202
sum(cqd_10 < cqthresh) / n
## [1] 0.801198
tc <- seq(0,21,by=0.001)
vl <- dgamma(tc, shape=2, scale=2)
par(mfrow=c(3,2))
plot(vl ~ tc, xlab="Time (days)", ylab="Viral load (arbitrary units)", yaxt='n')
hist(cqd, main="", xlab="PCR Cq,max values")
hist(Ed, main="", xlab="PCR Efficiency")
hist(cqd_10, main="", xlab="PCR Cq,max values (10x dilution)")
hist(r0, breaks=100, main="", xlab="Maximum viral load (arbitrary units)")
Change in sensitivity over dilutions
dat <- tibble(
D=seq(1,100),
sensitivity=sapply(D, function(x)
{
r0_x <- r0 / x
cqd_x <- est_cq(r0_x, Ed, rcq)
sum(cqd_x < cqthresh)/n
})
)
plot(sensitivity ~ D, dat)
Save the parameter values for the model
save(efficiency_params, Rcq, file="../data/efficiency_params.rdata")