GC1118 sensitivity analysis of exponent

library(RxODE)
library(tidyverse)
library(knitr)

Model definition

# Define model ##########
ode <- "
TBP = VP*RB/CLB;
CONC = A1/VC*0.16;
CPLXP =  ((KD*VP+A2+TBP)-((KD*VP+A2+TBP)^2-4*A2*TBP)^0.5)/2;
FAP = A2 - CPLXP;
RO = 100*CPLXP / (VP*RB/CLB);
d/dt(A1) = - A1*CLA/VC - Q*A1/VC + Q*FAP/VP;
d/dt(A2) =  -CPLXP*CLB/VP + Q*A1/VC -  Q*FAP/VP;
"
# Compile model ##########
mod1 <- RxODE(model = ode, modName = "mod1")

Sensitivity parameter setting

time.interval <- 1
test.dose <- c(50,100,250,400,500)
test.exp <- seq(0.75,0.85, 0.01)

PK parameters & simulation

data <- data.frame(time = seq(0, 21*24, time.interval),
                     CONC = 0)
for (i in test.dose) {
  for (j in test.exp) {
    DOSEmm3 = i   # Dose 
    All_ex = j    # Exponent
  BW = 70                      # Average bodyweight of monkey in 12 mg/kg group
  CLA_est = 1.16                     # Separet variable for senstivity analysis
  V_est = 128                  # Separet variable for senstivity analysis
  VP_est = 126                   # Separet variable for senstivity analysis
    
  
  CLA = CLA_est*(BW/3)**All_ex    # CL_ex  (mL/h)
  VC  = V_est*(BW/3)                    # Central volume of distribution (mL)
  Q    = 8.75 * (BW/3)**All_ex    # Intercompartmental clearance (mL/h)
  VP  = VP_est*(BW/3)           
  RB = 335 * (BW/3)**All_ex        # RB_ex  ; Target synthesis rate (pmol/h)
  CLB = 2.39*(BW/3)**All_ex        # Target & Complex clearance (mL/h)
  KD  = 4.5                            # Equilibrium binding constant (nM)
  TBP = VP*RB/CLB                  # Total binding target in peripheral
  # Simulate a single dose ##########
  MWlig = 150000  # Molecular weight of mAb
  DOSE = DOSEmm3 * 1.73
  DOSEnM = DOSE*1e-3/MWlig*1e12  #Dose as nanomole
  duration = 2   # infusion duration
  infrate = DOSEnM / duration  # inusion rate"
  sampling.day = 21  # sampling time as day
  sampling.time = 24*sampling.day   # sampling time as hour
  
params <- c(CLA = CLA, VC = VC, Q = Q, VP = VP, CLB = CLB, RB = RB, KD = KD)
inits <- c(0, 0)   
# Initialize event table ##########
ev <- eventTable()
ev$add.dosing(dose = DOSEnM, nbr.doses = 1, dosing.to =1, rate = infrate)#, dosing.interval = 168)
ev$add.sampling(seq(0, sampling.time, time.interval))
#integration
data2 <- as.data.frame(mod1$run(params, ev, inits))
data2 <- data2[,c(1,5)]
colnames(data2) <- c("time", paste("dose",i,"exp",j, sep="-"))
data  <- merge(data, data2, by = "time")
  }}

Data management for NCA

data <- data[, -2]
data.long <-
  data %>%
  gather(2:ncol(data), key = "group", value = "CONC")
data3 <- separate(data.long,
                 group, into = c("a", "dose", "b","exp"), sep = "-")
data3 <- data3[,c(1,3,5,6)]
#colnames(test) <- c("time","dose", "exp","conc")
unique(data3[,c("dose", "exp")])

nrow(unique(data3[,c("dose", "exp")]))

[1] 55

data3$NID <- paste(data3$dose, data3$exp, sep = "-")
data3$ID <- factor(data3$NID, label=1:length(unique(data3$NID)))
# Simulated Conc-Time profile by dose
ggplot(data3, aes(time, CONC, group = ID, colour = factor(exp))) + 
  geom_line() +
  scale_y_continuous(trans = "log10") +
  facet_wrap(~dose, ncol = 2, scales="free_y")

NCA analysis using NonCompart package and data management

nca_data <- NonCompart::tabNCA(data3, "ID", "time", "CONC")
nca_data <- as.data.frame(nca_data)
nca_AUC <- nca_data %>% select(ID, AUCIFO) %>% arrange(ID)
nca_AUC %>% knitr::kable()

ID	AUCIFO
1	2880.2896
2	2791.0786
3	2704.6087
4	2620.8024
5	2539.5781
6	2460.8623
7	2384.5782
8	2310.6513
9	2239.0111
10	2169.5871
11	2102.3122
12	13672.2490
13	13249.2327
14	12842.4774
15	12446.4734
16	12062.4019
17	11689.8248
18	11328.4733
19	10978.0702
20	10638.2945
21	10308.8649
22	9989.5197
23	28778.0455
24	27784.3324
25	26844.6364
26	25954.3595
27	25111.8013
28	24312.3893
29	23548.3455
30	22814.1826
31	22105.8150
32	21420.6236
33	20757.2407
34	1059.1385
35	1026.3169
36	994.5086
37	963.6836
38	933.8118
39	904.8641
40	876.8120
41	849.6285
42	823.2866
43	797.7608
44	773.0256
45	40699.8388
46	39226.1353
47	37816.8056
48	36471.3114
49	35193.2682
50	33973.7577
51	32819.6687
52	31728.0985
53	30694.9231
54	29714.2862
55	28779.7284

data_sub <- 
  data3 %>% 
  select(ID, dose, exp) %>% 
  arrange(ID)
data_sub <- data_sub[!duplicated(data_sub$ID),]
data_sub <- merge(data_sub, nca_AUC, by="ID")
#AUC obs
AUC_obs <- read.csv(file="Cetuxmab_AUC.csv", header = T, stringsAsFactors = F, nrow = 5)
colnames(AUC_obs) <- c("dose","AUCOBS")
## merge AUCobs and AUCsim
data_all <- merge(data_sub, AUC_obs, by="dose") %>%
  select(ID, dose, exp,  AUCSIM=AUCIFO, AUCOBS) %>%
  arrange(ID)
knitr::kable(data_all)

ID	dose	exp	AUCSIM	AUCOBS
1	100	0.75	2880.2896	2380
2	100	0.76	2791.0786	2380
3	100	0.77	2704.6087	2380
4	100	0.78	2620.8024	2380
5	100	0.79	2539.5781	2380
6	100	0.8	2460.8623	2380
7	100	0.81	2384.5782	2380
8	100	0.82	2310.6513	2380
9	100	0.83	2239.0111	2380
10	100	0.84	2169.5871	2380
11	100	0.85	2102.3122	2380
12	250	0.75	13672.2490	14600
13	250	0.76	13249.2327	14600
14	250	0.77	12842.4774	14600
15	250	0.78	12446.4734	14600
16	250	0.79	12062.4019	14600
17	250	0.8	11689.8248	14600
18	250	0.81	11328.4733	14600
19	250	0.82	10978.0702	14600
20	250	0.83	10638.2945	14600
21	250	0.84	10308.8649	14600
22	250	0.85	9989.5197	14600
23	400	0.75	28778.0455	19000
24	400	0.76	27784.3324	19000
25	400	0.77	26844.6364	19000
26	400	0.78	25954.3595	19000
27	400	0.79	25111.8013	19000
28	400	0.8	24312.3893	19000
29	400	0.81	23548.3455	19000
30	400	0.82	22814.1826	19000
31	400	0.83	22105.8150	19000
32	400	0.84	21420.6236	19000
33	400	0.85	20757.2407	19000
34	50	0.75	1059.1385	795
35	50	0.76	1026.3169	795
36	50	0.77	994.5086	795
37	50	0.78	963.6836	795
38	50	0.79	933.8118	795
39	50	0.8	904.8641	795
40	50	0.81	876.8120	795
41	50	0.82	849.6285	795
42	50	0.83	823.2866	795
43	50	0.84	797.7608	795
44	50	0.85	773.0256	795
45	500	0.75	40699.8388	30870
46	500	0.76	39226.1353	30870
47	500	0.77	37816.8056	30870
48	500	0.78	36471.3114	30870
49	500	0.79	35193.2682	30870
50	500	0.8	33973.7577	30870
51	500	0.81	32819.6687	30870
52	500	0.82	31728.0985	30870
53	500	0.83	30694.9231	30870
54	500	0.84	29714.2862	30870
55	500	0.85	28779.7284	30870

Calculate predictive performance (PPE, RMSE, AFE, AAFE)

# Individual predictive performance
data_all <-
  data_all %>%
  mutate(PPEi = (AUCSIM-AUCOBS)/AUCOBS,
         RMSEi = (AUCSIM-AUCOBS)^2,
         AFEi = log10(AUCSIM/AUCOBS),
         AAFEi = abs(AFEi))

package 愼㸱愼㹥bindrcpp愼㸱愼㹦 was built under R version 3.3.3

knitr::kable(data_all)

ID	dose	exp	AUCSIM	AUCOBS	PPEi	RMSEi	AFEi	AAFEi
1	100	0.75	2880.2896	2380	0.2102057	2.502897e+05	0.0828592	0.0828592
2	100	0.76	2791.0786	2380	0.1727221	1.689856e+05	0.0691951	0.0691951
3	100	0.77	2704.6087	2380	0.1363902	1.053708e+05	0.0555275	0.0555275
4	100	0.78	2620.8024	2380	0.1011775	5.798580e+04	0.0418573	0.0418573
5	100	0.79	2539.5781	2380	0.0670496	2.546517e+04	0.0281846	0.0281846
6	100	0.8	2460.8623	2380	0.0339757	6.538704e+03	0.0145103	0.0145103
7	100	0.81	2384.5782	2380	0.0019236	2.095956e+01	0.0008346	0.0008346
8	100	0.82	2310.6513	2380	-0.0291381	4.809240e+03	-0.0128425	0.0128425
9	100	0.83	2239.0111	2380	-0.0592390	1.987788e+04	-0.0265207	0.0265207
10	100	0.84	2169.5871	2380	-0.0884088	4.427358e+04	-0.0401999	0.0401999
11	100	0.85	2102.3122	2380	-0.1166755	7.711049e+04	-0.0538797	0.0538797
12	250	0.75	13672.2490	14600	-0.0635446	8.607219e+05	-0.0285129	0.0285129
13	250	0.76	13249.2327	14600	-0.0925183	1.824572e+06	-0.0421621	0.0421621
14	250	0.77	12842.4774	14600	-0.1203783	3.088886e+06	-0.0557040	0.0557040
15	250	0.78	12446.4734	14600	-0.1475018	4.637677e+06	-0.0693065	0.0693065
16	250	0.79	12062.4019	14600	-0.1738081	6.439404e+06	-0.0829191	0.0829191
17	250	0.8	11689.8248	14600	-0.1993271	8.469120e+06	-0.0965449	0.0965449
18	250	0.81	11328.4733	14600	-0.2240772	1.070289e+07	-0.1101815	0.1101815
19	250	0.82	10978.0702	14600	-0.2480774	1.311838e+07	-0.1238269	0.1238269
20	250	0.83	10638.2945	14600	-0.2713497	1.569511e+07	-0.1374808	0.1374808
21	250	0.84	10308.8649	14600	-0.2939134	1.841384e+07	-0.1511420	0.1511420
22	250	0.85	9989.5197	14600	-0.3157863	2.125653e+07	-0.1648082	0.1648082
23	400	0.75	28778.0455	19000	0.5146340	9.561017e+07	0.1803077	0.1803077
24	400	0.76	27784.3324	19000	0.4623333	7.716450e+07	0.1650464	0.1650464
25	400	0.77	26844.6364	19000	0.4128756	6.153832e+07	0.1501039	0.1501039
26	400	0.78	25954.3595	19000	0.3660189	4.836312e+07	0.1354567	0.1354567
27	400	0.79	25111.8013	19000	0.3216738	3.735411e+07	0.1211243	0.1211243
28	400	0.8	24312.3893	19000	0.2795994	2.822148e+07	0.1070740	0.1070740
29	400	0.81	23548.3455	19000	0.2393866	2.068745e+07	0.0932068	0.0932068
30	400	0.82	22814.1826	19000	0.2007465	1.454799e+07	0.0794513	0.0794513
31	400	0.83	22105.8150	19000	0.1634639	9.646087e+06	0.0657529	0.0657529
32	400	0.84	21420.6236	19000	0.1274012	5.859419e+06	0.0520785	0.0520785
33	400	0.85	20757.2407	19000	0.0924864	3.087895e+06	0.0384160	0.0384160
34	50	0.75	1059.1385	795	0.3322497	6.976914e+04	0.1245856	0.1245856
35	50	0.76	1026.3169	795	0.2909646	5.350750e+04	0.1109143	0.1109143
36	50	0.77	994.5086	795	0.2509542	3.980367e+04	0.0972414	0.0972414
37	50	0.78	963.6836	795	0.2121806	2.845414e+04	0.0835673	0.0835673
38	50	0.79	933.8118	795	0.1746060	1.926872e+04	0.0698922	0.0698922
39	50	0.8	904.8641	795	0.1381939	1.207013e+04	0.0562162	0.0562162
40	50	0.81	876.8120	795	0.1029082	6.693200e+03	0.0425393	0.0425393
41	50	0.82	849.6285	795	0.0687150	2.984268e+03	0.0288619	0.0288619
42	50	0.83	823.2866	795	0.0355807	8.001340e+02	0.0151839	0.0151839
43	50	0.84	797.7608	795	0.0034726	7.621750e+00	0.0015055	0.0015055
44	50	0.85	773.0256	795	-0.0276407	4.828729e+02	-0.0121732	0.0121732
45	500	0.75	40699.8388	30870	0.3184269	9.662573e+07	0.1200561	0.1200561
46	500	0.76	39226.1353	30870	0.2706879	6.982500e+07	0.1040389	0.1040389
47	500	0.77	37816.8056	30870	0.2250342	4.825811e+07	0.0881482	0.0881482
48	500	0.78	36471.3114	30870	0.1814484	3.137469e+07	0.0724147	0.0724147
49	500	0.79	35193.2682	30870	0.1400476	1.869065e+07	0.0569230	0.0569230
50	500	0.8	33973.7577	30870	0.1005428	9.633312e+06	0.0416070	0.0416070
51	500	0.81	32819.6687	30870	0.0631574	3.801208e+06	0.0265976	0.0265976
52	500	0.82	31728.0985	30870	0.0277972	7.363331e+05	0.0119074	0.0119074
53	500	0.83	30694.9231	30870	-0.0056714	3.065192e+04	-0.0024701	0.0024701
54	500	0.84	29714.2862	30870	-0.0374381	1.335674e+06	-0.0165713	0.0165713
55	500	0.85	28779.7284	30870	-0.0677121	4.369235e+06	-0.0304499	0.0304499

# Mean by exponent
data_all_ex <-
  data_all %>%
  group_by(exp) %>%
  summarise(PPE=sum(PPEi)/5*100,
            RMSE = sqrt(sum(RMSEi)/5),
            AFE = 10^(mean(AFEi)),
            AAFE = 10^(mean(AAFEi)))
knitr::kable(data_all_ex)

exp	PPE	RMSE	AFE	AAFE
0.75	26.2394339	6219.593	1.2469790	1.280160
0.76	22.0837925	5459.607	1.2061644	1.253924
0.77	18.0975178	4754.587	1.1669799	1.228415
0.78	14.2664707	4110.035	1.1292703	1.203706
0.79	10.5913777	3536.351	1.0930525	1.179801
0.8	7.0596964	3044.422	1.0582117	1.156619
0.81	3.6659714	2653.234	1.0247062	1.134154
0.82	0.4008633	2383.715	0.9924537	1.125584
0.83	-2.7443107	2253.554	0.9613755	1.120680
0.84	-5.7777268	2265.092	0.9313956	1.127975
0.85	-8.7065656	2399.635	0.9024455	1.148009

Plot the result.

data_all_ex_long <-
  data_all_ex %>%
  gather(2:5, key = "evaluation", value = "performance")
ggplot(data_all_ex_long, aes(exp, performance)) + 
  geom_point() +
  facet_wrap(~evaluation, scales="free_y")