## simulate the PK over the study

rm(list = ls()) # clear the environment

graphics.off()

options(stringsAsFactors=FALSE)

library(ggplot2)

library(dplyr)

library(mrgsolve)

library(magrittr)

select <- dplyr::select

plot.me <- TRUE

set.seed(-1)

auc.f = function(xIn, yIn)

## set up the 1-compartment model

code <- '$PARAM TVCL = 1.6, TVV=28, TVKA=0.6, DOSE=-1

$PKMODEL ncmt=1, depot=TRUE

$CMT GUT CENT

$SET req="GUT CENT"

$GLOBAL

$MAIN

double KA = TVKA;

double CL = TVCL*exp(IIVCL);

double V = TVV*exp(IIVV);

$OMEGA @annotated

IIVCL: 0.1 : CL

IIVV: 0.1 : V

$TABLE

double IPRED = CENT/V*1e3;//convert mg/L to ng/mL

double DV = IPRED;//DV is ng/mL

double DAY = floor(TIME/24)+1;

$CAPTURE IPRED DV DAY KA CL V DOSE'

## load the PK model

mod <- mcode("pk_1cmt", code)

## generate a PK data set

tbl.sim <- expand.ev(amt=c(0,10,30,100),id=1:10,addl=27,ii=24) %>% mutate(DOSE=amt)

## here is the time grid to 30 days

des <- tgrid(0,24*30-1e-10,0.1)

## simulate the PK model

out <- mod %>% mrgsim(data=tbl.sim,tgrid=des)

## calculate PK metrics by day

tbl.pk <- out %>%

as.data.frame %>%

group_by(ID,DAY,DOSE) %>%

summarise(

AUC = auc.f(xIn=time, yIn=IPRED),

CMAX = max(IPRED),

CMIN = min(IPRED)

) %>%

mutate(CAVG=AUC/24)