R语言实现非房室模型算法

一粒沙

发布于 2019-08-21 18:05:18

3.6K00

代码可运行

文章被收录于专栏：R语言交流中心R语言交流中心

运行总次数：0

代码可运行

药代动力学分析过程中房室模型和非房室模型成为两大主要分支。房室模型分析法的基础是把机体以类群形式分为几个不同的隔室或房室，然后根据药物在各房室间的转运或消除速率常数建立能够反应药物在机体内的变化规律的数学模型。其参数的估测都是依据房室模型而进行的。非房室方法不需要对药物或代谢物设定专门的房室。事实上，只要药物符合线性药物动力学，那不管它属于什么样的隔室模型，都能采用此法。同时非房室方法是处理药物在体内分布和消除不规则的药物动力学分析的主要手段。尽管非房室模型可以覆盖所有的房室模型，同时在用于药物浓度非特异性测定方面优于房室模型，但是目前房室模型已成为药代动的金标准。总之，两者各有所长。今天我们主要给大家介绍在R语言中如何实现非房室模型分析。我们需要用到R包PKNCA。

首先包的安装需要:install.packages(“PKNCA”)。当然想更加完美的可视化计算结果，需要加载其他的包：

suppressPackageStartupMessages({
 library(PKNCA)
 library(dplyr)
 library(cowplot)
 library(knitr)
 library(ggplot2)
})

数据的载入及可视化：

my_conc <- data.frame(conc=c(0, 2.5, 3,2, 1.5, 1.2, 1.1, 0, 0),
                      time=c(0:5, 8, 12, 24),
                      subject=1)
my_conc$BLOQ <- my_conc$conc == 0
my_conc$measured <- TRUE

ggplot(my_conc,
      aes(x=time,
          y=conc,
          shape=BLOQ,
          group=subject)) +
 geom_line() +
 geom_point(size=4) +
 scale_x_continuous(breaks=my_conc$time) +
 theme(legend.position=c(0.8, 0.8))

接下来就是包中数据的标准化：

1. PKNCAconc 构建模型输入的数据。我们直接看下构建实例：

conc_obj <- PKNCAconc(my_conc,conc~time|subject)

2. PKNCAdose 数据构建。实例：

conc_obj <-PKNCAconc(as.data.frame(datasets::Theoph), conc~Time|Subject)
 
## Dosing data needs to only have one rowper dose, so subset for
## that first.
d_dose <-unique(datasets::Theoph[datasets::Theoph$Time == 0,
                                 c("Dose", "Time", "Subject")])
knitr::kable(d_dose,
            caption="Example dosing data extracted from theophylline dataset")

dose_obj <- PKNCAdose(d_dose,Dose~Time|Subject)

3. PKNCAdata 主要是整合PKNCAconc以及PKNCAdose数据。看下实例：

data_obj <-PKNCAdata(data.conc=conc_obj,
                     intervals=data.frame(start=0,
                                          end=24,
                                          aucall=TRUE,
                                          auclast=TRUE,
                                          aucinf.pred=TRUE,
                                           aucinf.obs=TRUE))

通过上面的参数intervals可以说设置我们要计算的NCA参数以及对应的时间间隔。如果提供了对应的剂量数据，就可以直接进行整合：

data_obj <- PKNCAdata(conc_obj, dose_obj)

以上数据的具体分布查看可以通过如下的实例：

knitr::kable(PKNCA.options("single.dose.aucs"))

查看所有的subject子项的参数列表：

knitr::kable(data_obj_automatic$intervals)

接下来看下包中主要算法：

1. pk.calc.tfirst和pk.calc.tlast最低定量浓度第一个时间点和最后一个点。实例：

tfirst<- pk.calc.tfirst (conc=my_conc$conc,time=my_conc$time)
 
tlast <-pk.calc.tlast(conc=my_conc$conc, time=my_conc$time)
ggplot(my_conc,
      aes(x=time,
          y=conc,
          shape=BLOQ,
           group=subject)) +
 geom_ribbon(data=my_conc[my_conc$include_auclast,],
              aes(ymin=0, ymax=conc),
              fill="lightblue") +
 geom_line() +
 geom_point(size=4) +
 scale_x_continuous(breaks=my_conc$time) +
 theme(legend.position=c(0.8, 0.8))

2. pk.nca 主要的非房室模型算法函数。实例：

results_obj <- pk.nca(data_obj)
knitr::kable(head(as.data.frame(results_obj)))

summary( results_obj)

如果利用添加intervals参数的数据进行建模。我们看下会有什么差异：

intervals_manual <- data.frame(start=0,
                               end=Inf,
                               cmax=TRUE,
                              tmax=TRUE,
                               aucinf.obs=TRUE,
                               auclast=TRUE)
data_obj_manual <- PKNCAdata(conc_obj,dose_obj,
                            intervals=intervals_manual)
knitr::kable(data_obj_manual$intervals)

 results_obj_manual <-pk.nca(data_obj_manual)
knitr::kable(head(as.data.frame(results_obj_manual)))

summary( results_obj_manual)

那么，从上面可以看出其实计算出来的结果跟是否选择的参数列表有关。但主要看数据是否满足。当然，他还有更高级的模式，那就是多组数据的引入。利用intervals可以设置多个时间间隔之间的数据计算。

intervals_manual <- data.frame(start=0,
                               end=c(24,Inf),
                               cmax=TRUE,
                              tmax=c(FALSE,TRUE),
                               aucinf.obs=TRUE,
                               auclast=TRUE)
data_obj_manual <- PKNCAdata(conc_obj,dose_obj,
                            intervals=intervals_manual)
knitr::kable(data_obj_manual$intervals)

最后，此包获取累积的药物浓度构建了对应的处理函数superposition。其中tau为时间间隔，n.tau剂量间隔的个数

但是此处需要注意的是药物浓度必须从零开始。也就是在时间点0时浓度为0。

## Correct nonzero concentrations at time 0to be BLQ.
theoph_corrected <- datasets::Theoph
theoph_corrected$conc[theoph_corrected$Time== 0] <- 0
conc_obj_corrected <-PKNCAconc(theoph_corrected, conc~Time|Subject)

1. 对从单剂量给药到稳态过程的累积浓度

## Calculate the new steady-stateconcentrations with 24 hour dosing
steady_state <-superposition(conc_obj_corrected, tau=24)
knitr::kable(head(steady_state, n=14),
            caption="Superposition at steady-state")

2. 对单剂量数据到指定剂量过程的累积浓度

## Calculate the unsteady-stateconcentrations with 24 hour dosing
unsteady_state <-superposition(conc_obj_corrected, tau=24, n.tau=2)
knitr::kable(head(unsteady_state, n=14),
            caption="Superposition before steady-state")

3. 自定义给药时间间隔的累积浓度计算

## Calculate the new steady-stateconcentrations with 24 hour dosing
complex_interval_steady_state <-superposition(conc_obj_corrected, tau=24, dose.times=c(0, 2, 4))
knitr::kable(head(complex_interval_steady_state,n=24),
            caption="Superposition at steady-state with complex dosing")

借助ggplotbao我们可以对我们的数据结果进行可视化展示：

ggplot(complex_interval_steady_state,
      aes(y=conc, x=time, colour=Subject)) +
 geom_line()

接下来为了确定达到稳定状态所需的浓度曲线，可以给出达到稳定状态所需的所有剂量时间。见实例：

up_to_steady_state <-superposition(conc_obj_corrected,
                                    tau=4*24,
                                    n.tau=1,
                                   dose.times=seq(0, 3*24, by=12))
ggplot(up_to_steady_state, aes(x=time,y=conc, colour=Subject)) +
 geom_line()