PK_modelling.Rmd

---
title: "pharmacokinetics modelling"
author: "liuc"
date: "1/12/2022"
output: pdf_document
editor_options: 
  markdown: 
    wrap: 72
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## pharmacokinetics modelling by R

> <http://sia.webpopix.org/pharmacokinetics.html>
> <http://sia.webpopix.org/nlme.html>
>
> 熟悉nlmixr2包的文档
>
> <https://www.pkpd-info.com/PK_introduction.php>
> <http://opensource.nibr.com/xgx/index.html>
> https://cran.r-project.org/web/views/Pharmacokinetics.html 
> https://discuss.go-isop.org/
> https://cran.r-project.org/web/packages/ubiquity/vignettes/NCA.html



### 房室模型

Once a drug is administered, we usually describe subsequent processes
within the organism by the pharmacokinetics (PK) process known as ADME:
absorption, distribution, metabolism, excretion.
具体的关于PK的一些描述可以参见第四个链接。

*什么是PK的房室模型*

将人体简化为depot compartment。PK compartment models assume that drug
concentration is perfectly homogeneous in each compartment of the body
at all times.

PK models do not describe how a drug is administered. Administered doses
are source terms, i.e., inputs that dynamically modify the state of a
system. This is important to note because it means that the same PK
model can be used for different dose regimens.
典型的PK房室模型包括中心房室（central compartment）和外周房室（peripheral compartments），中心房室代表血浆或血液循环系统，而外周房室代表组织或器官。通过建立和求解一组微分方程，PK房室模型可以描述药物在这些房室之间的动态分布和消除过程。这些方程描述了药物在各个房室之间的转移速率、药物在每个房室中的浓度以及药物的消除速率。中心房室（central compartment）通常代表血浆或血液循环系统，它是药物在体内的主要分布区域。药物在中心房室中的浓度可以反映血液中的药物浓度。外周房室（peripheral compartments）代表体内的组织或器官。药物从中心房室进入外周房室，通过血液循环分布到组织或器官中。外周房室的数量可以根据需要进行扩展，以更准确地描述药物在不同组织或器官之间的分布和转移。



*ordinary differential equations (ODEs): 常微分方程*
ODEs在PK房室模型中用于建立药物在中心房室（通常为血浆）和外周房室（代表组织或器官）之间的药物转移过程。通过这些方程，我们可以模拟和预测药物在体内的吸收、分布和消除过程，以及药物在不同组织和器官中的浓度变化。

利用ODEs，我们可以计算和推断出药物在给定剂量和给药方案下的浓度-时间曲线，从而了解药物在体内的动态行为。此外，ODEs还可以用于优化给药方案、评估药物的疗效和毒性，以及预测个体化的药物剂量调整。


*常见的一些指标：*
Vd（分布容积）：表示药物在体内分布的广度。它是药物总体量与血浆中药物浓度之间的比率。Vd可以反映药物在体内的分布程度，通常以升或升/公斤为单位。
Cl（清除率）：表示单位时间内从体内清除药物的速率。清除率可以用来评估药物的消除速度，通常以升/小时或升/分钟为单位。
t1/2（半衰期）：表示药物浓度减少一半所需的时间。半衰期是描述药物消除速度的重要参数，可以用来评估药物在体内的停留时间和剂量频率。
AUC（曲线下面积）：表示药物浓度-时间曲线下的面积，反映了药物在体内的总曝露程度。AUC可用于评估药物的吸收、分布和消除过程，并与药物的疗效和毒性相关。
Cmax（峰浓度）：表示药物在给药后达到的最高浓度。Cmax可以用来评估药物的吸收速度和最大效应。

k : rate constants (kinetic) -- movement from 'x' to 'y'. (消除速率常数）：k表示药物从体内清除的速率常数。它决定了药物浓度下降的速度。k的倒数（1/k）被称为药物的消除半衰期（t1/2）。消除速率常数k通常以小时的倒数（1/h）为单位。
Ka: absorption constant. 吸收速率常数）：ka表示药物从给药部位吸收到血液循环的速率常数。它用于描述药物的吸收过程。ka通常以小时的倒数（1/h）为单位。
Vc: distribution compartment (central). （中心房室容积）：Vc表示中心房室（通常为血浆）的容积。它是用来表示药物在中心房室的分布范围。Vc的单位通常是升或升/公斤。
Q: rate constant between Vc and Vp. （外周房室间转移速率常数）：Q表示不同外周房室之间药物转移的速率常数。它用于描述药物在不同组织或器官之间的分布和转移。Q通常以升/小时为单位。
Vp: distribution compartment (periph.) 
CL:elimination constant (clearance) Described by differential equations of change in compartment (A) over time (T): dA/dT
At:（时间到达最大浓度）：这个指标表示药物在体内吸收过程中所需的时间，通常以时间单位表示（例如小时）。At表示从给药开始到药物在体内达到最高浓度（Cmax）所经过的时间。

Ktp（血浆半衰期）：血浆半衰期是描述药物从体内消除一半所需的时间。它是指药物浓度减少到初始浓度的一半所需的时间。通常，血浆半衰期以时间单位（例如小时）表示。

Kpt（分布常数）：Kpt是PK房室模型中的一个参数，表示药物分布过程的速度常数。它描述了药物从血浆向组织的转移速度。Kpt的单位通常是1/时间（例如小时的倒数）。

tka（时间到达最低浓度）：tka表示药物在体内吸收过程中所需的时间，以达到最低浓度（Cmin）。类似于At（时间到达最大浓度），tka表示从给药开始到药物在体内达到最低浓度所经过的时间。

tcl（血浆消除半衰期）：tcl是指血浆中药物浓度减少到初始浓度的一半所需的时间。它表示药物在体内消除过程的速度。tcl通常以时间单位（例如小时）表示。

tv（中央分布容积）：tv是药物动力学模型中的一个参数，表示药物分布到中央组织（通常是血浆或中央血池）的容积。它代表药物在体内的分布范围。tv的单位通常是容积单位（例如升）。



*非线性混合效应模型*

对于非线性混合效应模型，The system will support but not be limited to the following estimation methods: FO, FOI, FOCE, FOCEI, Laplacian, Lindstrom and Bates, MCMC, MCPEM, SAEM, Gaussian quadrature, and nonparametric methods。



*One-compartment model（一室模型）：*

在一室模型中，药物被视为在体内分布均匀。这意味着药物在血液中的浓度与组织中的浓度是一致的。这个简化模型对于那些迅速在血浆和组织之间分布的药物较适用。但对于那些在体内分布不均匀的药物，一室模型可能不能准确描述其药代动力学特征。

*Two-compartment model（两室模型）：*

两室模型包括两个独立的"室"：中心室（central compartment）和外周室（peripheral compartment）。中心室通常包括血浆、心脏、肝脏和肾脏等高血流组织；外周室通常包括肌肉、皮肤和脂肪等较低血流组织。在这个模型中，药物首先在中心室分布，然后通过血流进入外周室。两室模型较好地描述了那些在体内分布不均匀的药物。

*Three-compartment model（三室模型）：*

三室模型进一步细分了药物在体内的分布过程，将外周室分为两个不同的室。这种模型适用于那些在体内存在多个不同分布相的药物。例如，药物可能首先在中心室分布，然后迅速进入一个外周室（比如较高血流的组织），再慢慢进入另一个外周室（比如较低血流的组织）。这种模型可以更精确地描述那些在体内具有复杂分布特征的药物。



```{r, include=FALSE}
library(tidyverse)
# library(saemix) # https://saemixr.github.io/
# library(nlmixr) # 主要的R包，不过和2的版本不兼容，二者先library一个

# nlmixr2 is an R package for fitting population pharmacokinetic (PK) and pharmacokinetic-pharmacodynamic (PKPD) models.
library(nlmixr2) # 目前看起来是一个整合了多种nlme的包
library(rxode2) # ode-based models
library(nlmixr2data)

library(ggPMX)
library(xpose)

# library('nonmem2R')
library(xgxr) # 诺华 提供的一个探索性R包

```


### 拟合PK数据的房室模型

```{r}
require(pmxTools)
require(nonmem2rx)
require(babelmixr2)
```


#### single-dose 单次给药

##### fitting a simple one-compartment PK model use `nlmixr2` package.

one-compartment指的是，In a one-compartment model, the body is considered as a single well-stirred compartment where the drug distributes instantaneously and uniformly. The model assumes that the drug is rapidly and completely absorbed into the systemic circulation and that the elimination occurs directly from this central compartment.

下面的模型中，和在使用`saemix`包时是一致的，包括模型和起始值。`nlmixr2`包的语法和ini的命名系统是模仿的NONMEM软件。在ini中设置值时一般的是一个值，或者三个值（c(lower, est, upper)）,`~`是方差或相关性，`<-`是参数赋值,`#`后面的内容是
label的内容会在Parameter展现，也可以用label函数。

其中，model部分的意思是：A mathematical model which best describes the data.

`theo_sd`是常见的theophylline数据集，是longitudinal repeated measures
data，每一个个体其随时间测得的血药浓度在PK 中都是非线性的。 The form of
the nonlinear model is determined by the pharmacokinetic theory, not
derived from the data.
但是依据PK理论所得到的每一个个体的PK参数是有变化的，不同的，所以还要考虑population，引入个体间的random
effect。
在`nlmixr2`中数据的列名是有要求的吗？毕竟文档的示例中直接就输入了数据名。在saemix中尚且需要指定数据。


```{r}
# https://nlmixr2.org/articles/running_nlmixr.html
# 数据集的介绍来自官网
theo_sd
```

> https://nlmixr2.github.io/rxode2/articles/rxode2-datasets.html

Subject Identification Columns
The subject identification column separates subjects for identification of random effects.

ID: A subject identifier that may be an integer, character, or factor.
Observation Columns
Observation columns are used to indicate the dependent variable and how to use or measure it.

DV: A numeric column with the measurement
CENS: A numeric column for indication of censoring, such as below the limit of quantification for an assay.
LIMIT: A numeric column for helping indicate the type of censoring, such as below the limit of quantification for an assay.
MDV: An indicator for missing DV values
CMT: The name or number of the compartment
DVID: The dependent variable identifier
EVID The event identifier
Dosing Columns
AMT: The amount of the dose
CMT: The name or number of the compartment
EVID: The event identifier
ADDL: The number of additional doses
RATE or DUR: The rate or duration of a dose



```{r}
# 罗列所有的Est方法
nlmixr2est::nlmixr2AllEst()
```


Once the initialization block has been defined, you can define a model in terms of the variables defined in the ini block. You can also mix rxode2() blocks into the model if needed.
```{r}
#  ordinary differential equations (ODEs) model
# 采用不同的算法需要不同的模型参数
# 同时注意文档中所介绍的命名格式的问题，比如rx_ or nlmixr_开头是不允许的，不过一般也都是模仿NONMEM的命名
one.compartment.saem <- function() {
    ini({
        tka <- .5   # Log Ka
        tcl <- -3.2 # Log Cl
        tv <- -1    # Log V
        eta.ka ~ 1 ## BSV Ka
        eta.cl ~ 2
        eta.v ~ 1
        add.err <- 0.1 # residual variability
    })
    model({
        ka <- exp(tka + eta.ka)
        cl <- exp(tcl + eta.cl)
        v <- exp(tv + eta.v)
        d/dt(depot) = -ka * depot
        d/dt(center) = ka * depot - cl / v * center
        cp = center / v 
        cp ~ add(add.err)
    })
}


# Checking model syntax
nlmixr2::nlmixr(one.compartment.saem)
```

fit the model:
fit模型有几种方法，此处的方法，在不提供data参数时，会打印one.compartment.saem的内容。
```{r}
# 数据集相关的变量啥也不用指定吗？CDISC?
fit <- nlmixr2(one.compartment.saem, 
               data = theo_sd, 
               est="saem" # nlmixr2est::nlmixr2AllEst()
               )
```



在终端打印更漂亮：
```{r}
print(fit)
```


base plot:

```{r}
## Standard nlmixr plots
plot(fit)
```
*interpret the result: *




除上述`plot`函数得到的诊断图之外，还提供了`xpose`包的使用，`xpose`包是aims
to reduce the post processing burden and improve diagnostics commonly
associated the development of non-linear mixed effect models. 检查模型的表现。

```{r}
library(xpose.nlmixr2)

xpdb <- xpose_data_nlmixr(fit) # 此object可以使用xpose包的函数

xpose::dv_vs_ipred(xpdb)
```


Another option is to use the ggPMX package. 
```{r}
library(ggPMX)
ctr = pmx_nlmixr(fit)
pmx_plot_dv_ipred(ctr)
```


##### 两房室模型

The “true” model is two compartmental, with first-order absorption and linear elimination.

```{r}
Bolus_2CPT %>% head(n = 100)
```


###### 利用`saemix`包

详情见`saemix`包的官方文档。

```{r}
saemix.data <- saemixData(
  name.data = warfa_data, header = TRUE, sep = " ",
  na = NA, name.group = c("id"), name.predictors = c("amount", "time"),
  name.response = c("y1"), name.X = "time"
)
```


```{r}
model1cpt <- function(psi, id, xidep) {
  dose <- xidep[, 1]
  tim <- xidep[, 2]
  ka <- psi[id, 1]
  V <- psi[id, 2]
  k <- psi[id, 3]
  CL <- k * V
  ypred <- dose * ka / (V * (ka - k)) * (exp(-k * tim) - exp(-ka * tim))
  return(ypred)
}

saemix.model <- saemixModel(
  model = model1cpt, description = "warfarin",
  type = "structural", psi0 = matrix(c(1, 7, 1, 0, 0, 0),
    ncol = 3, byrow = TRUE,
    dimnames = list(NULL, c("ka", "V", "k"))
  ), transform.par = c(1, 1, 1),
  omega.init = matrix(c(1, 0, 0, 0, 1, 0, 0, 0, 1), ncol = 3, byrow = TRUE),
  covariance.model = matrix(c(1, 0, 0, 0, 1, 0, 0, 0, 1),
    ncol = 3,
    byrow = TRUE
  )
)
```


###### 利用`nlmixr2`包进行

```{r}
# 两房室模型对应的起始值
# 模型的微分方程和参数

#定义两室模型公式

# Two-compartment model with additive residual error:
model.2cpt.ode <- function() {
  ini({
    tka <- log(1.05)
    tcl <- log(0.121)
    tv2 <- log(1.939)
    tv3 <- log(5.65)
    tq <- log(0.282)
    eta.ka ~ 0.1
    eta.cl ~ 0.1
    eta.v2 ~ 0.1
    eta.v3 ~ 0.1
    eta.q ~ 0.1
    add.err <- 75
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v2 <- exp(tv2 + eta.v2)
    v3 <- exp(tv3 + eta.v3)
    q <- exp(tq + eta.q)
    d / dt(depot) <- -ka * depot
    d / dt(center) <- ka * depot - cl / v2 * center + q / v3 * periph - q / v2 * center
    d / dt(periph) <- q / v2 * center - q / v3 * periph
    cp <- center / v2
    cp ~ add(add.err)
  })
}

model.2cpt.ode_2 <- function() {
  ini({
    tka <- log(1.05)           # 吸收速率常数的对数值
    tcl <- log(0.121)          # 清除速率常数的对数值
    tv2 <- log(1.939)          # 第二个组织的体积常数的对数值
    tv3 <- log(5.65)           # 第三个组织的体积常数的对数值
    tq <- log(0.282)           # 第一个组织到第二个组织的转移速率常数的对数值
    eta.ka ~ 0.1               # 吸收速率的个体差异（ETA）的标准差
    eta.cl ~ 0.1               # 清除速率的个体差异（ETA）的标准差
    eta.v2 ~ 0.1               # 第二个组织的体积的个体差异（ETA）的标准差
    eta.v3 ~ 0.1               # 第三个组织的体积的个体差异（ETA）的标准差
    eta.q ~ 0.1                # 第一个组织到第二个组织的转移速率的个体差异（ETA）的标准差
    add.err <- 75              # 额外的误差（残差误差）常数值
  })
  model({
    ka <- exp(tka + eta.ka)    # 计算吸收速率
    cl <- exp(tcl + eta.cl)    # 计算清除速率
    v2 <- exp(tv2 + eta.v2)    # 计算第二个组织的体积
    v3 <- exp(tv3 + eta.v3)    # 计算第三个组织的体积
    q <- exp(tq + eta.q)       # 计算第一个组织到第二个组织的转移速率
    d / dt(depot) <- -ka * depot                       # 第一个组织（depot）的变化率
    d / dt(center) <- ka * depot - cl / v2 * center + q / v3 * periph - q / v2 * center   # 第二个组织（center）的变化率
    d / dt(periph) <- q / v2 * center - q / v3 * periph  # 第三个组织（periph）的变化率
    cp <- center / v2          # 计算中央组织（center）中的药物浓度
    cp ~ add(add.err)          # 建立中央组织中药物浓度与残差误差之间的关系
  })
}


# Two-compartment model with proportional residual error:
model.2cptp.ode <- function() {
  ini({
    tka <- log(1.05)
    tcl <- log(0.121)
    tv2 <- log(1.939)
    tv3 <- log(5.65)
    tq <- log(0.282)
    eta.ka ~ 0.1
    eta.cl ~ 0.1
    eta.v2 ~ 0.1
    eta.v3 ~ 0.1
    eta.q ~ 0.1
    prop.err <- 0.075
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v2 <- exp(tv2 + eta.v2)
    v3 <- exp(tv3 + eta.v3)
    q <- exp(tq + eta.q)
    d / dt(depot) <- -ka * depot
    d / dt(center) <- ka * depot - cl / v2 * center + q / v3 * periph - q / v2 * center
    d / dt(periph) <- q / v2 * center - q / v3 * periph
    cp <- center / v2
    cp ~ prop(prop.err)
  })
}


# Two-compartment model with proportional residual error and weight on CL:
model.2cpt.ode.wtcl <- function() {
  ini({
    tka <- log(1.14)
    tcl <- log(0.0190)
    tv2 <- log(2.12)
    tv3 <- log(20.4)
    tq <- log(0.383)
    wteff <- 0.35
    eta.ka ~ 1
    eta.cl ~ 1
    eta.v2 ~ 1
    eta.v3 ~ 1
    eta.q ~ 1
    prop.err <- 0.075
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + wteff * lnWT + eta.cl)
    v2 <- exp(tv2 + eta.v2)
    v3 <- exp(tv3 + eta.v3)
    q <- exp(tq + eta.q)
    d / dt(depot) <- -ka * depot
    d / dt(center) <- ka * depot - cl / v2 * center + q / v3 * periph - q / v2 * center
    d / dt(periph) <- q / v2 * center - q / v3 * periph
    cp <- center / v2
    cp ~ prop(prop.err)
  })
}


# Two-compartment model with proportional residual error and sex on 𝑉𝑉2:
model.2cpt.ode.sexv2 <- function() {
  ini({
    tka <- log(1.14)
    tcl <- log(0.0190)
    tv2 <- log(2.12)
    tv3 <- log(20.4)
    tq <- log(0.383)
    sexeff <- -0.2
    eta.ka ~ 1
    eta.cl ~ 1
    eta.v2 ~ 1
    eta.v3 ~ 1
    eta.q ~ 1
    prop.err <- 0.075
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v2 <- exp(tv2 + sexeff * (SEX) + eta.v2)
    v3 <- exp(tv3 + eta.v3)
    q <- exp(tq + eta.q)
    d / dt(depot) <- -ka * depot
    d / dt(center) <- ka * depot - cl / v2 * center + q / v3 * periph - q / v2 * center
    d / dt(periph) <- q / v2 * center - q / v3 * periph
    cp <- center / v2
    cp ~ prop(prop.err)
  })
}
```


```{r}
# Checking model syntax
nlmixr2::nlmixr(model.2cpt.ode)
```

模型的拟合过程非常耗时，要注意数据的大小，不知能否多线程运行。
对于此数据集，在我的Mac上跑了得有一个多小时。
```{r}
# fit.2cpt.ode.saem <- nlmixr2(object = model.2cpt.ode,
#                              data = Bolus_2CPT,
#                              est = 'saem'
#                              )

# write_rds(fit.2cpt.ode.saem, './datasets/fit.2cpt.ode.saem.rds')

fit.2cpt.ode.saem <- readRDS('./datasets/fit.2cpt.ode.saem.rds')
```


```{r}
# 运行也很久
AIC(fit.2cpt.ode.saem)
```


利用`xpose`包对模型的拟合程度做检查：
```{r}
# require(xpose.nlmixr2)

xpdb <- xpose.nlmixr2::xpose_data_nlmixr2(fit.2cpt.ode.saem)

xpdb
```

```{r}
# 似乎有点小问题
summary(xpdb, problem = 1)
```

数据拟合的很差，对于这种拟合很差的情况该如何解决呢？
```{r}
dv_vs_ipred(xpdb)
```



##### three compartment model


```{r}
# Three-compartment model with proportional residual error:

model.3cpt.ode <- function() {
  ini({
    tka <- log(1.42)
    tcl <- log(0.044)
    tv2 <- log(2)
    tv3 <- log(10)
    tv4 <- log(10)
    tq2 <- log(0.5)
    tq3 <- log(0.5)
    eta.ka ~ 0.1
    eta.cl ~ 0.1
    eta.v2 ~ 0.1
    eta.v3 ~ 0.1
    eta.v4 ~ 0.1
    eta.q2 ~ 0.1
    eta.q3 ~ 0.1
    prop.err <- 0.075
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v2 <- exp(tv2 + eta.v2)
    v3 <- exp(tv3 + eta.v3)
    v4 <- exp(tv4 + eta.v4)
    q2 <- exp(tq2 + eta.q2)
    q3 <- exp(tq3 + eta.q3)
    d / dt(depot) <- -ka * depot
    d / dt(center) <- ka * depot - cl / v2 * center + q2 / v3 * periph1 - q2 / v2 * center + q3 / v4 *
      periph2 - q3 / v2 * center
    d / dt(periph1) <- q2 / v2 * center - q2 / v3 * periph1
    d / dt(periph2) <- q3 / v2 * center - q3 / v4 * periph2
    cp <- center / v2
    cp ~ prop(prop.err)
  })
}


# 定义模型
model.3cpt.ode <- function() {

  # 初始值
  ini({
    # 吸收速率常数的对数值
    tka <- log(1.42)  
    # 清除率常数的对数值 
    tcl <- log(0.044)   
    # 中央室体积的对数值
    tv2 <- log(2) 
    # 外周室1体积的对数值 
    tv3 <- log(10)  
    # 外周室2体积的对数值
    tv4 <- log(10) 
    # 从中央室到外周室1的间质清除速率常数的对数值
    tq2 <- log(0.5) 
    # 从中央室到外周室2的间质清除速率常数的对数值
    tq3 <- log(0.5)
    # 吸收速率的随机效应标准差
    eta.ka ~ 0.1 
    # 清除率的随机效应标准差 
    eta.cl ~ 0.1    
    # 中央室体积的随机效应标准差
    eta.v2 ~ 0.1
    # 外周室1体积的随机效应标准差
    eta.v3 ~ 0.1  
    # 外周室2体积的随机效应标准差
    eta.v4 ~ 0.1
    # 从中央室到外周室1的间质清除速率随机效应标准差
    eta.q2 ~ 0.1
    # 从中央室到外周室2的间质清除速率随机效应标准差
    eta.q3 ~ 0.1
    # 归一化残差的标准差
    prop.err <- 0.075
  })

  # 模型方程 
  model({
    # 吸收速率
    ka <- exp(tka + eta.ka)
    # 清除率
    cl <- exp(tcl + eta.cl) 
    # 中央室体积
    v2 <- exp(tv2 + eta.v2)
    # 外周室1体积
    v3 <- exp(tv3 + eta.v3)
    # 外周室2体积
    v4 <- exp(tv4 + eta.v4)
    # 从中央室到外周室1的间质清除速率 
    q2 <- exp(tq2 + eta.q2) 
    # 从中央室到外周室2的间质清除速率
    q3 <- exp(tq3 + eta.q3)

    # 释放部位量的微分方程
    d/dt(depot) = -ka * depot

    # 中央室量的微分方程
    d/dt(center) = ka * depot - cl / v2 * center +  
                   q2/v3 * periph1 - q2/v2 * center + 
                   q3/v4 * periph2 - q3/v2 * center

    # 外周室1量的微分方程
    d/dt(periph1) = q2/v2 * center - q2/v3 * periph1

    # 外周室2量的微分方程 
    d/dt(periph2) = q3/v2 * center - q3/v4 * periph2

    # 中央室浓度
    cp = center / v2

    # 浓度的误差模型
    cp ~ prop(prop.err)
  })

}

```


#### multi-dose 多次给药

```{r}
theo_md
```




### use PKNCA 包


*NCA的实验设计要点：*
此包主要是提供non-compartmental analysis (NCA) PK计算方法, NCA is used
as method of description of PK with minimal assumptions of the rates of
distribution of the drug through the body. NCA is typically used to
describe the PK of a drug in clinical studies with many samples per
subject on the same and sequential days.

complete sampling and sparse sampling designs具体指的是，

以我这两天查资料的情况来看，NCA和上文的modeling方法还是不一样的。顾名思义，一个把人体理解为compartment，一个没有。

PKNCA需要两个数据集，分别是浓度数据集和dose数据集，对于每个数据集中的分类变量，其的语法为：concentration, dose, and time must all be numeric.


*NCA的一些资料收集：*

只要药物符合线性药物动力学，不管它属于什么样的隔室模型，都能使用非房室方法。非房室模型是临床计算PK参数的主要方法。
WinNonlin中非房室模型分析，以概率论和数理统计学中的统计矩方法为理论基础，对数据进行解析，主要包括零阶矩和一阶矩。

其中不依赖末端消除速率λz的参数有药峰浓度（Cmax），达峰时间（Tmax），最后一个可定量点所对应的时间点（Tlast）和药物浓度（Clast），这些参数不依赖于计算，而是直接由观测值得出。还有从0时刻到最后一个可定量点所对应的曲线下面积（AUClast），一阶矩平均驻留时间（MRT），反应了药物在体内的平均停留时间。


*一些概念的理清：*
AUClast(Area Under the Curve)：AUClast是药物生物利用度描述指标，也就是逐点法计算得到的从给药开始到最后测量得到测定浓度（Clast）的药物暴露程度。即药物在体内的暴露程度，反映药物的吸收程度。AUClast是指药物在最后一个可测量时间点（通常为最后一个采样时间点）前的曲线下面积。它反映了药物在体内的总体曝露程度，即药物在体内的浓度随时间的积累情况。AUClast可以通过计算药物浓度与时间的曲线下的面积来获得。

Cmax(Maximum Concentration)：药物吸收后，体内浓度达到最大值即为Cmax。一般来说，Cmax越大，说明药物吸收越快越全。它是一个用来描绘给药后药物浓度-时间曲线（concentration-time curve）的重要参数。 Cmax通常用来评估药物的吸收速度和程度，对于药物的疗效和安全性具有重要意义。*Cmax似乎就是时间梯度中最大的那个值，且来自观测数据本身，不是模型拟合计算得来的。*

Clast：Clast是指在给药后的最后测定时刻所测得的药物浓度，也就是药物在体内吸收过程的最后测定点的血药浓度。clast通常用于评估药物在体内的消除速度和半衰期。*其也是直接提取的raw data中的数据*

Tmax(Time to Maximum Concentration)：这是药物达到最大浓度的时间。这通常是一个反映药物吸收速度的重要指标。*来自数据本身*

AUC0-∞：这是从给药开始到无穷远的预测药物全身暴露程度。平均来说，AUC0-∞相比AUClast能更好地反映药物的全部体内暴露程度。它反映了药物在体内的总体曝露程度，包括药物的吸收、分布、代谢和排泄过程。

T1/2(Half-life) ：这是药物半衰期，表示药物浓度减半所需的时间。这个数字反映了药物从体内消除的速度。

CL(Clearance)：清除率，衡量体内清除药物的能力。它可以反映药物的代谢和排泄情况。清除率越高，药物在体内的停留时间越短。
$CL=Dose/AUC$

Vd(volume of distribution)：分布容积，反映药物在体内分布的广度。

MRT(Mean residence time): $MRT=AUMC/AUC$, 表示药物在体内平均停留的时间

Vss(Steady state volume of distribution): $Vss = MRT × CL$ , 稳态分布容积）是一个用于描述药物在体内分布的参数，通常用符号Vss表示。它表示在稳态条件下，使血浆中的药物浓度达到平衡所需的药物总量与血浆中的药物浓度之间的关系。

"lambda" 通常用于表示药物的消除速率常数（elimination rate constant）。消除速率常数是衡量药物从体内清除的速度的指标，它表示在单位时间内药物浓度减少的比例。

"tlag" 在药物药代动力学中是指药物在给药后进入循环系统的延迟时间（lagtime）。它表示从给药到药物开始出现在血液中测量到的浓度的时间间隔。

AUMCinf（Area Under the Moment Curve at Infinity）是指药物的无限时刻矩时刻曲线下的面积。 它是PK（Pharmacokinetics，药物药代动力学）中的一个衡量参数，用于描述药物在体内的总体暴露程度和分布。AUMCinf是通过计算药物浓度-时间曲线的矩时刻曲线下的面积来获得的。矩时刻曲线是将时间的各个次幂乘以对应时间点上的药物浓度所得到的曲线。AUMCinf通常与AUCinf（Area Under the Concentration Curve at Infinity，无限时刻浓度曲线下的面积）一起使用，以提供更全面的药物暴露和药代动力学信息。

Span ratio（跨度比）是在非组分分析（NCA）中用于评估药物的体内分布特性的指标之一。它是通过计算药物在血浆和体液（如组织或尿液）中的浓度之间的比值来得到的。跨度比的计算通常采用AUC（Area Under the Curve，曲线下面积）作为血浆和体液浓度曲线的度量，并分别计算血浆和体液的AUC值。然后，通过计算体液AUC与血浆AUC的比值，得到跨度比。跨度比大于1表示药物在体液中的浓度高于血浆中的浓度，表明药物在体液中有较高的蓄积。跨度比小于1表示药物在体液中的浓度低于血浆中的浓度，表明药物在体液中的分布较少。

Interpolation & Extrapolation, 插入/推断。
> https://en.wikipedia.org/wiki/Curve_fitting


*Generally NCA will determine the following directly from the data:*
These properties are all based on observational data.

Cmax - Maximum observed concentration (units=concentration)
Tmax - The time where the maximum concentration was observed (units=time)
AUC - The area under the curve (units=time × concentration)
AUMC - The area under the first moment curve (units = time2 × concentration)

*即然是基于观察数据得到的值，其本身只是总体的一个近似，这里变凸显出来实验设计的重要性，合理的时间点的设置*



*NCA采用的非线性模型的具体方法：*
_以下内容来自chatGPT，不一定正确_
PK NCA（Noncompartmental Analysis）本身并不涉及非线性统计模型的拟合。它是一种基于数值积分和统计方法的非参数分析，用于计算和估计药物在体内的曝露和药代动力学参数，而无需假设具体的药物分布模型。

PK NCA主要基于药物的浓度-时间数据，通过计算面积下曲线（AUC）和其他参数来评估药物的曝露程度、消除速率、吸收速度等。这些参数的计算通常不需要采用非线性统计模型。

非线性统计模型通常在药物动力学建模和药物代谢动力学（PK/PD）建模中使用，用于描述药物在体内的动力学过程，考虑吸收、分布、代谢和排泄等因素，并拟合实际数据。这些模型可以是线性模型或非线性模型，如生理药动模型（Physiologically Based Pharmacokinetic Model，PBPK）、双指数模型、Michaelis-Menten模型等。这些模型需要更复杂的参数估计方法，通常使用非线性最小二乘拟合等技术进行参数拟合。

WinNonlin中的PK NCA计算基于数值积分和统计方法，以计算和估计药物在体内的曝露和药代动力学参数。曲线下面积（AUC）计算：WinNonlin采用梯形法则（Trapezoidal Rule）进行AUC的计算。它通过对浓度-时间数据进行数值积分，将浓度-时间曲线下的面积近似为一系列梯形的面积之和。


*三种数据类型：*
1. Sparse noncompartmental analysis (NCA)，理解起来似乎就是拼凑起来的数据，在这种情况下，由于缺乏每个个体的完整曲线数据， 传统的非组分分析方法可能无法直接应用。稀疏NCA方法旨在解决这个问题，以在有限的数据情况下对药物的药代动力学进行评估。
2. Single dose
3. Multiple dose



```{r}
 # 同样可以用来处理PK数据的R包,其所提供的NCA计算也是基于PKNCA的,而且其提供的一些文档还是不错的
# library(ubiquity)
# 其提供的sparse NCA分析脚本值得参考，待完成

library(PKNCA)


# 其他的一些有用的R包
library(PKData)

# 韩国一个机构出品的包
require(NonCompart)
require(ncar)
```


#### Single dose(单次给药)

single dose 可以理解为visit时期开始的时间，比如C1D1, C0D7等等。但是数据里面的dose还依旧是可以多个的。

当然按照PKNCA的文档，真正只有一个dose的时候，`PKNCA.options("single.dose.aucs")`

美博所提供的数据:

```{r}
my_pk <- readxl::read_excel('~/OneDrive/kintor/Daily_Work/meibo_PK_pooled/8411-751_Concentration Table_01Jul2019.xls',
                            sheet = 'Concentration_individual',
                            skip = 3
                            ) %>% 
  janitor::remove_empty() %>% janitor::clean_names()
```

```{r}
my_pk1 <- my_pk %>% select(1:4) %>% 
  slice(-1)

my_pk2 <- my_pk %>% select(nominal_time_h:x16) %>% 
  set_names(slice(., 1)) %>% 
  slice(-1)

```

BQL编码为0，missing value为NA，
```{r}
pk_data <- my_pk1 %>% bind_cols(my_pk2) %>% 
  fill(c(visit,analyte,dose_level_mg),.direction = 'down')

pk_GT0918 <- pk_data %>% filter(analyte == 'GT0918',
                                visit == 'C1D1',
                                dose_level_mg == '200'
                                ) %>% 
  pivot_longer(cols = -(visit:subject),
               names_to = 'Time',
               values_to = 'conc'
               ) %>% 
  mutate(Time = as.numeric(Time),
         BLQ = if_else(conc == 'BQL', TRUE, FALSE),
         conc = ifelse(conc == 'BQL', 0, conc),
         conc = as.numeric(conc)
         )
```



*12例患者，每一个患者都在同一个dose浓度下，测量11个时间点的血清浓度数据。*
```{r}
# 132个观测
head(datasets::Theoph)

Theoph <- Theoph %>% as.data.frame()
```

```{r}
# 好多好多可供选择的内容
PKNCA.options()

# 可以设置的intervals
get.interval.cols()

#
?pk.calc.auc
?pk.calc.aumc
?pk.calc.half.life
```


单位设置
```{r}
d_units_auto <- pknca_units_table(concu="ng/mL", doseu="mg", amountu="mg", timeu="hr")

# 其可以对单位进行自动的转换
# 参数传给PKNCAdata函数
d_units <-
  pknca_units_table(
    concu="mg/L", doseu="mg/kg", timeu="hr",
    # use molar units for concentrations and AUCs
    conversions=
      data.frame(
        PPORRESU=c("(mg/kg)/(hr*mg/L)", "(mg/kg)/(mg/L)", "mg/L", "hr*mg/L"),
        PPSTRESU=c("L/hr/kg", "L/kg", "mmol/L", "hr*mmol/L"),
        conversion_factor=c(NA, NA, 1/180.164, 1/180.164)
      )
  )
```



以下计算是按照ID group列分开计算的，Parameter calculation will automatically split the data by the grouping factor(s), subset by the interval, calculate all required parameters, record all options used for the calculations, and include data provenance to show that the calculation was performed as described. 
```{r}

# 不手动设置则会自动设置
# 自动设置比较稳健，在你还没有全部了解的情况下
# Intervals are subsets within a group by start and end time.
# interval可以理解为所有时间梯度中的一部分，当然也可以选择全部
intervals <-
  data.frame(start=0, end=c(24, Inf),
             cmax=c(FALSE, TRUE),
             tmax=c(FALSE, TRUE),
             auclast=TRUE,
             aucinf.obs=c(FALSE, TRUE),
             mrt.iv.obs = c(FALSE, TRUE),
             mrt.iv.last = c(TRUE, FALSE),
             aumclast = TRUE,
             half.life = TRUE
             )

# 不同采样时间点的药物浓度
# 没有考虑不同的dose
conc_obj <- PKNCA::PKNCAconc(pk_GT0918, conc~Time|subject)

#
# d_dose <- unique(datasets::Theoph[datasets::Theoph$Time == 0,
#                                   c("Dose", "Time", "Subject")])

# 本步的目的只是为了得到每个患者的dose信息而已
d_dose <- unique(pk_GT0918[pk_GT0918$Time == 0,
                                  c("dose_level_mg", "Time", "subject")])

# dose 和 时间点
# ?PKNCAdose
dose_obj <- PKNCAdose(d_dose, dose_level_mg~Time|subject,
                      route = 'extravascular',#Define the route of administration."extravascular" or "intravascular" 
                      )

# ?PKNCAdata
data_obj_automatic <- PKNCAdata(conc_obj, 
                                dose_obj,
                                intervals = intervals
                                )

intervals_manual <- data_obj_automatic$intervals
```


```{r}
results_obj_automatic <- pk.nca(data_obj_automatic)

# 输出的结果符合
# pharmacokinetic parameter (PP) TESTCD of CDISC SDTM
results_obj_automatic
# results_obj_automatic$result
```


```{r}
# Plot the concentration-time data and the interval
ggplot(pk_GT0918, aes(x=Time, y=conc)) +
  geom_point() + geom_line() +
  scale_y_continuous(limits=c(0, NA)) +
  labs(x="Time Since First Dose (hr)",
       y="Concentration\n(arbitrary units)")
```

```{r}
PKNCA.set.summary(
  name="half.life",
  description="arithmetic mean and standard deviation",
  point=business.mean,
  spread=business.sd,
  rounding=list(signif=3)
)

summary(results_obj_automatic) |> as.data.frame()

# summary提供的是对上表中多个样本的统计值，包括几何均值等。
# 比如以下只对auclast进行summary计算
results_obj_automatic$result %>% filter(PPTESTCD == 'auclast') %>% 
  summarise(across(PPORRES, c(PKNCA::geomean, PKNCA::geocv)))
```


*导出宽数据：*
```{r}
# as.data.frame(results_obj_automatic) == results_obj_automatic$result
results_obj_automatic$result %>% filter(end == 'Inf') %>% pivot_wider(id_cols = c(ru_zu_bian_hao, start, end),
                                                                      names_from = PPTESTCD,
                                                                      values_from = PPORRES
                                                                      )
```



*以下还可以单独计算：*

AUC to the Last Value Above the Limit of Quantification (AUC~last~),
AUClast是指，在tlast（最后一个非零值所对应的时间点）到时间点为0的曲线下面积。

```{r}
# 首先是找到tlast
pk_GT0918_1subject  <- pk_GT0918 %>% filter(subject == 1021)

tlast <- pk.calc.tlast(conc=pk_GT0918_1subject$conc,
                       time=pk_GT0918_1subject$Time)
tlast
```


AUC~all~
在tlast后的一个点进行linear interpolation to zero后，求得的从时间0开始的曲线下积分。

```{r}
first_after_tlast <- my_conc$time[my_conc$time > tlast][1]
first_after_tlast
```


AUC ~0-inf~

commonly used for single-dose data. It calculates the AUC0-last and then extrapolates to ∞ using the estimated half-life.
文档中的图很好的说明了这一AUC所对应的面积。
和AUC_all_相比（可以简单的理解为曲线下的面积），其更像是从half-life推测出的延长线而得到的曲线下面积，一般而言其比AUC~all~的值要大一些。

```{r}

```


Half-Life Calculation

> https://billdenney.github.io/pknca/articles/v06-half-life-calculation.html

半衰期的计算是对concentration取自然对数的log后和时间拟合曲线。The default calculation method is curve stripping。
似乎就是选择不同的点拟合不同的模型，然后选择一个best fit。

```{r}
PKNCA.options("min.hl.points")

PKNCA.options("allow.tmax.in.half.life")
```



*use another package:*
其的结果和WinNonlin 的结果保持一致，尤其是AUC的部分。
其对只有一个患者的数据，和多个患者的数据是用不同的函数分析的。
```{r}
NonCompart::tblNCA(Theoph, key = "Subject", colTime = "Time", colConc = "conc", 
                   dose = 320, # Dose对不同的患者不一致咋办呢？每个样本一个浓度是可行的
                   adm = "Extravascular", dur = 0, doseUnit = "mg", timeUnit = "h", concUnit = "mg/L", 
                   down = "Linear")
```


```{r}
# 注意目前需要data.frame格式，tibble的话需要改一下
# 而且如果数据质量不行的话，会计算出很多的NA值

Time = Theoph[Theoph$Subject == 8, "Time"]
Concentration = Theoph[Theoph$Subject == 8, "conc"]
Res = sNCA(Time, Concentration,dose = 320, concUnit = "mg/L")
IntAUC(Time, Concentration, t1 = 0.5, t2 = 11, Res)
```


```{r}
ncar::pdfNCA(fileName = "pdfNCA-Theoph.pdf", Theoph, key = "Subject", 
       colTime = "Time", colConc = "conc", dose = 320, doseUnit = "mg", 
       timeUnit = "h", concUnit = "mg/L", down = "Linear")

# ncar::pdfNCA(fileName = "~/OneDrive/kintor/Daily_Work/GT1708F_CN2001/pdfNCA-GT1708F.pdf", 
#              cal_pk_c0d1, key = "ru_zu_bian_hao", 
#              colTime = "ji_hua_shi_jian_dian", colConc = "pk_nong_du", 
#              dose = aa, 
#              doseUnit = "mg", timeUnit = "h", concUnit = "mg/L", 
#              down = "Linear"
#              )
# 
# ncar::rtfNCA(fileName = "~/OneDrive/kintor/Daily_Work/GT1708F_CN2001/pdfNCA-GT1708F.rtf", 
#              cal_pk_c0d1, key = "ru_zu_bian_hao", 
#              colTime = "ji_hua_shi_jian_dian", colConc = "pk_nong_du", 
#              dose = aa, 
#              doseUnit = "mg", timeUnit = "h", concUnit = "mg/L", 
#              down = "Linear"
#              )
```



#### Multi dose（多次给药）

对于Multi dose的理解就是，选择最后一期visit时期。

multiple-dose study个人目前倾向于理解为，一个患者采集了多个dose的数据。即是多次给药。

多次给药后的稳态作何理解呢？

```{r}
?superposition
```