问题描述
我正在尝试使用 deSolve 包中的函数 ode 来求解一阶微分方程.问题如下:药物在某些时间(输注时间)以恒定输注速率给药，并以一级速率消除.因此，该过程可以描述为:
if(t %in% Infusion_times){Infusion 其中 t 是时间，Infusion_times 是包含给药时间的向量，C 是药物的数量，Ke是它的消除常数，Infusion是一个变量，有输液时取输液速度的值，否则取0.因此，假设我们要进行 3 次给药，从 0、24 和 40 次开始，每次输注持续两个小时，并且我们希望 deSolve 每 0.02 个时间单位计算一次答案.例如，我们希望 deSolve 以 0.02 次单位的步长求解 0 到 48 次之间的微分方程.因此，我们指定 ode 函数的时间的向量将是:
times 输液时间由下式给出:
Infusion_times 起初，我担心问题可能出在浮点数的比较上.为了防止它，我将两个向量四舍五入到小数点后两位:
times 现在，希望所有 Infusion_times 都包含在 times 向量中:
(sum(Infusion_times %in% times)/length(Infusion_times))*100[1] 100如您所见，Infusion_times 中的所有值(100%)都包含在向量times 中，因此变量Infusion应该在指定的时间取值 Infusion_rate.然而，当我们解方程时，它不起作用.让我们证明一下，但首先，让我们指定 ode 函数所需的其他值:
参数 现在，让我们根据需要编写一个函数来说明变化率:
OneComp 对于那些不熟悉 deSolve 包的人，函数的参数 t 将说明方程应该被积分的次数，amounts 将指定 C 的初始值，parameters 将给出参数的值(在这种情况下，只是 Ke).现在，让我们解方程:
out <- ode(func = OneComp, y = amount, parms = parameters, times = times)让我们绘制结果:
库(ggplot2)ggplot(as.data.frame(out)) + geom_line(aes(x = times, y = C))如果 Infusion 总是等于 0，我们会得到完全相同的结果.但是，请注意，如果我们只模拟单一剂量，并且我们要尝试类似的方法，它会起作用:
OneComp 这里，我们让变量Infuse只在时间小于2小时时取Inf_rate的值，它起作用了！因此，我对这些行为感到非常困惑.这不是改变变量值的问题，也不是浮点数之间比较的问题......知道这些是什么吗?谢谢
我一直在为同样的问题苦苦挣扎一段时间.我试图使用 deSolve 包而不是原始模型中使用的 RxODE 包复制 IV 输注，然后进行 PO 给药.我的解决方案是列出输液时间，然后提取最大值:
tmp 此处，Num.Doses 设置为 5 次 IV 输注.Tau(给药间隔)为 24 小时，0 是输液开始时间，2 是输液结束时间，单位为小时.
接下来我构建了一个模型，其完整版本是来自 RxODE 的多室 PKPD 模型(
I am trying to solve a first-order differential equation using the function ode from the deSolve package. The problem is as follows: a drug is administered by a constant infusion rate at some times (infusion times) and eliminated in a first-order rate. Thus, the process can be described by:
if(t %in% Infusion_times){Infusion <- Infusion_rate} else{Infusion <- 0}
dC <- -Ke*C + Infusion
where t is the time, Infusion_times is a vector containing at what times the drug is administered, C is the quantity of the drug, Ke is its elimination constant and Infusion is a variable taking the value of the infusion rate when there is infusion, and the value 0 otherwise.So, let's say we want to administer 3 doses, starting at times 0, 24 and 40, with each infusion lasting for two hours, and let's say we want deSolve to compute the answer every 0.02 units of time.We want deSolve to solve the differential equation for the times between 0 and 48 with steps of 0.02 times unit, for instance. So our vector specifying the times to the ode function would be:
times <- seq(from = 0, to = 48, by = 0.02)
The infusion times are given by:
Infusion_times <- c(seq(from = 0, to = 2, by = 0.02), seq(from = 24, to = 26, by = 0.02),
seq(from = 40, to = 42, by = 0.02))
At first, I was concerned that the problem could be in the comparison of floating-points. To prevent it, I rounded both vectors to two decimal places:
times <- round(times, 2)
Infusion_times <- round(times, 2)
So now, hopefully, all Infusion_times are included in the times vector:
(sum(Infusion_times %in% times)/length(Infusion_times))*100
[1] 100
As you can see, all the values in Infusion_times (100%) are contained in the vector times, and thus the variable Infusion should take the value Infusion_rate at the specified times. However, when we solve the equation, it doesn't work.Let's prove it, but first, let's specify the other values needed by the ode function:
parameters <- c(Ke = 0.5)
amounts <- c(C = 0) #Initial value for drug is 0
Inf_rate <- 5
And now, let's write a function stating the rate of changes, just as required:
OneComp <- function(t, amounts, parameters){
with(as.list(c(amounts, parameters)),{
if(t %in% Infusion_times){Infuse =Inf_rate} else{Infuse = 0}
dC <- -Ke*C + Infuse
list(c(dC))})
}
For those who are not familiar with the deSolve package, the argument t of the function will state the times for which the equation should be integrated, amounts will specify the initial value of C and parameters will give the value of the parameters (in this case, just Ke).And now, let's solve the equation:
out <- ode(func = OneComp, y = amounts, parms = parameters, times = times)
Let's plot the results:
library(ggplot2)
ggplot(as.data.frame(out)) + geom_line(aes(x = times, y = C))
This is exactly the same result we would get if Infusion was always equal to 0.However, note that if we were to mimic only a single dose, and we were to try a similar approach, it would work:
OneComp <- function(t, amounts, parameters){
with(as.list(c(amounts, parameters)),{
if(t < 2){Infuse =Inf_rate} else{Infuse = 0}
dC <- -Ke*C + Infuse
list(c(dC))})
}
out <- ode(func = OneComp, y = amounts, parms = parameters, times = times)
ggplot(as.data.frame(out)) + geom_line(aes(x = times, y = C))
Here, we have made the variable Infuse take the value of the Inf_rate only when the time is less than 2 hours, and it works!Therefore, I am totally puzzled with these behaviour. It is not a problem of changing the value of a variable, it is not a problem of a comparison between floating-point numbers... Any idea of what could these be? Thanks
I have been struggling with the same problem for some time. I was trying to replicate IV infusion followed by PO dosing using the deSolve package rather than the RxODE package used in the original model. My solution was to make a list of infusion times and later extract a maximal value:
tmp <- c()
for (i in seq(from = 1, to = Num.Doses, by = 1)) {
tmp[i] <- (i * Tau) - Tau
tmp2 <- seq(from = 0, to = 2, by = 0.1)
Inf.Times <- round(unlist(lapply(tmp, function(x, Add) x + tmp2)), 3)}
Here, Num.Doses is set for 5 IV infusions. The Tau (dosing interval) is 24 hours, 0 is infusion start time, and 2 is infusion end time in hours.
Next I constructed a model, the full version of which is a multicompartment PKPD model from RxODE (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4728294/):
Model <- function(time, yini, parameters) {
with(as.list(c(yini, parameters)), {
Infusion <- ifelse(time %% Tau < Inf.Time & time <= max(Inf.Times), Dose / Inf.Time, 0)
C1 <- Central / V1
dDepot <- - (Ka * Depot)
dCentral <- (Ka * Depot) - (CL * C1) + Infusion
list(c(dDepot, dCentral))})}
I got the idea for the Infusion <- line from Andrew Booker as shown at https://github.com/andrewhooker/PopED/issues/11. I modified his suggestion from if else to ifelse and incorporated a hard stop after the end of the 5th infusion with the time <= max(Inf.Times) otherwise the model would continually re-infuse. That also allowed me to implement additional non-IV dosing using an events table when running the model using deSolve:
Dose.Events <- data.frame(var = "Depot", time = c(120, 144, 168, 192, 216), value = Dose2, method = "add")
Times <- seq(from = 0, to = Tot.Hours, by = 1)
out <- ode(y = Ini.Con, times = Times, func = Model2, parms = Parms, events = list(data = Dose.Events))
When run using the full model, the output is nearly the same as using the original code in RxODE, which is more straightforward and "clean". The differences, as judged by AUC, were minimal with sig figs the same out to 6 digits. I am able to replicate the IV infusions (first set of 5 peaks) and also recapitulate the PO dosing (second set of 5 peaks).
这篇关于R- ode 函数(deSolve 包):将参数值更改为时间的函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！