Case study: how long must PSA be followed to estimate regrowth?

library(pkident)

A treatment-sensitive PSA component declines while a resistant component grows: the ODE form of the biexponential decline-and-regrowth model used for PSA and tumour size (for example Stein et al., Clin Cancer Res 2008).

m <- pkpd_model(
  "d/dt(S) = -d*S
   d/dt(R) = g*R
   psa = S + R",
  outputs = "psa", initial = list(S = "s0", R = "r0")
)

The parameter values below are illustrative, not estimates from data: PSA 20 ng/mL at the start of treatment, 2.5% from the resistant component, sensitive half-life about 3.5 weeks (d = 0.8 per month), resistant doubling time about 4.6 months (g = 0.15 per month). PSA reaches its nadir at about 5.6 months.

values <- c(s0 = 19.5, r0 = 0.5, d = 0.8, g = 0.15)
error <- residual_error(add = 0.05, prop = 0.15)

Structure

From PSA alone, all four parameters are locally identifiable:

structural_identifiability(m)
#> [intdy -- t = 4.92728e-310 illegal]: 1 warning(s) for subject(s): 1
#> [intdy -- t = 4.92728e-310 illegal]: 1 warning(s) for subject(s): 1
#> [intdy -- t = 6.9532e-310 illegal]: 3 warning(s) for subject(s): 1
#> <pkident structural identifiability>
#> Decision is LOCAL (Lie derivatives at the known initial state up to order 3, checked by output sensitivities; 5 random points).
#> Methods agree at 5 of 5 points.
#>   d            identifiable
#>   g            identifiable
#>   r0           identifiable
#>   s0           identifiable

Length of observation

One patient sampled monthly from the start of treatment, for 3 to 24 months:

lengths <- c(3, 6, 9, 12, 18, 24)
candidates <- setNames(lapply(lengths, function(L) experiment(seq(0, L, 1))),
                       paste0(lengths, " months"))
cmp <- compare_designs(m, candidates, values = values, error = error, structural = FALSE)
cmp
#> <pkident design comparison>
#> Expected RSE (%) from the Fisher information; structural status for the observed outputs (LOCAL).
#>   parameter      3 months   6 months   9 months  12 months  18 months  24 months 
#>   d                 118.0       20.4       12.5       10.3        8.9        8.4 
#>   g                5298.3      200.3       47.1       20.6        7.9        4.4 
#>   r0               3177.7      182.9       55.5       29.5       15.2       10.7 
#>   s0                 78.1       12.7       12.6       12.5       12.3       12.3 
#>   identifiable        0/4        2/4        2/4        4/4        4/4        4/4 
#> Practically identifiable: RSE <= 30%.

p <- as_table(cmp, "parameters")
g <- p[p$parameter == "g", ]
d <- p[p$parameter == "d", ]
plot(lengths, g$rse_percent, log = "y", type = "b", pch = 19, ylim = c(3, 6000),
     xlab = "Months of follow-up", ylab = "Expected RSE (%)")
lines(lengths, d$rse_percent, type = "b", pch = 1, lty = 2)
abline(h = 30, col = "grey50")
abline(v = 5.6, col = "grey50", lty = 3)
legend("topright", c("regrowth rate g", "decline rate d", "RSE 30%", "nadir"),
       lty = c(1, 2, 1, 3), pch = c(19, 1, NA, NA), col = c(1, 1, "grey50", "grey50"), bty = "n")

The decline rate is determined by 6 months. The regrowth rate is not determined until PSA has been followed well past the nadir: its expected RSE falls below 30% only at about 12 months, twice the time to nadir.

Do the expected RSEs hold?

The expected RSEs come from a local, linear approximation. For 12 and 24 months, 1,000 patients were simulated and each was fitted by maximum likelihood (case-studies/psa-regrowth.R in the source repository):

chk <- read.csv(system.file("extdata", "psa-simulation-check.csv", package = "pkident"))
knitr::kable(chk[, c("months", "parameter", "expected_rse", "empirical_rse", "ratio")], digits = 2)

months	parameter	expected_rse	empirical_rse	ratio
12	d	10.29	10.83	1.05
12	g	20.62	21.81	1.06
12	r0	29.48	30.97	1.05
12	s0	12.50	12.98	1.04
24	d	8.40	8.89	1.06
24	g	4.36	4.42	1.01
24	r0	10.66	10.92	1.02
24	s0	12.27	12.73	1.04

The spread of the estimates is within 6% of the expected RSE.