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)

Fitting realistic Cq values

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")