2016-09-06 47 views
-2

我試圖模擬R中的細胞攝取,並從伯克利麥當娜移植了一個模型。該模型由幾個常量和微分方程組成,用於計算量和濃度。代碼的一部分列出:如何基於多個組構建R中的ODE

library(deSolve) 

fb = 0.0510 
Km = 23.5 
Pdif = 0.429 
Vmax = 270 
Vol_cell = 9.33 
Vol_media = 150 
S = 10 #concentration of dosing media 

yini = c(Amt_media=(S*Vol_media)-(S*fb*Vol_cell), 
     Amt_cell=S*fb*Vol_cell, 
     Amt_total=S*Vol_media, 
     Con_media=S-(S*fb), 
     Con_cell=S*fb) 

Uptake = function(t, y, p){ 
    dy1 = (- (Pdif * y[1]) + (Pdif * y[2]) - ((Vmax * y[4])/(Km + y[4]))) 
    dy2 = (+ (Pdif * y[1]) - (Pdif * y[2]) + ((Vmax * y[4])/(Km + y[4]))) 
    dy3 = dy1 + dy2 
    dy4 = dy1/Vol_media 
    dy5 = dy2/Vol_cell 
    list(c(dy1, dy2, dy3, dy4, dy5))} 

times1 = seq(from=0, to=15, by=0.01) 
out1 = ode(y=yini, times=times1, func=Uptake, parms=NULL, method="rk4") 

其餘代碼用於輸出到數據框和繪圖。我的問題是如何讓代碼結構化爲使用「S」作爲幾個濃度的列表,以便每個濃度都可以應用到微分方程(基本上給我一個out1用於S1,out2用於S2等等,然後可以被傳遞到數據框)?在伯克利麥當娜,這是通過寫35多個微分方程實現的,儘管如果可能的話我想在R中使用簡化的方法。

+1

寫一個函數,'S'作爲參數,然後使用' lapply'來調用它以獲得許多不同的'S'值:myfun <-function(S){return S + 1};拉普利(1:10,myfun)'。這會給你一個對應於'S'的每個值的輸出列表。 – MrFlick

+0

謝謝MrFlick。我一直在翻閱SO和R文檔幾個星期,現在正在尋找示例代碼來整合到模型中。無可否認,我不太清楚如何實施您的解決方案。您的建議是否被納入「攝取」功能中作爲「with(parm ...)」,或作爲獨立功能存在,強制「攝取」進入循環以對所有「S」的值產生「out」? – Ron

+0

''uptake''似乎與'S'沒有任何關係,你只能在你的'yini'中使用它。你是否剛剛開始編程?是這個問題嗎?也許一個基本的R教程可能會幫助你我不確定我是否理解你陷入困境 – MrFlick

回答

0

使用S的唯一部分是yini值的初始化。基本上我們只需要將那部分和那些運行ode的部分移動到一個新函數中。然後,您可以調用該函數以獲得您想要的任何值。例如,

#set up 
library(deSolve) 

fb <- 0.0510 
Km <- 23.5 
Pdif <- 0.429 
Vmax <- 270 
Vol_cell <- 9.33 
Vol_media <- 150 

Uptake <- function(t, y, p){ 
    dy1 = (- (Pdif * y[1]) + (Pdif * y[2]) - ((Vmax * y[4])/(Km + y[4]))) 
    dy2 = (+ (Pdif * y[1]) - (Pdif * y[2]) + ((Vmax * y[4])/(Km + y[4]))) 
    dy3 = dy1 + dy2 
    dy4 = dy1/Vol_media 
    dy5 = dy2/Vol_cell 
    list(c(dy1, dy2, dy3, dy4, dy5))} 

times1 <- seq(from=0, to=15, by=0.01) 

# function with S as a parameter 
runConc <- function(S) { 
    yini <- c(Amt_media=(S*Vol_media)-(S*fb*Vol_cell), 
    Amt_cell=S*fb*Vol_cell, 
    Amt_total=S*Vol_media, 
    Con_media=S-(S*fb), 
    Con_cell=S*fb) 

    ode(y=yini, times=times1, func=Uptake, parms=NULL, method="rk4") 
} 

#run for concentrations 10,20,30 
out <- lapply(c(10,20,30), runConc) 

這將產生一個列表對象,其中包含每個濃度的結果。所以out[[1]]是S = 10,out[[2]]結果S=20等,我們可以看到每個結果的前幾行與

lapply(out, head, 3) 

# [[1]] 
#  time Amt_media Amt_cell Amt_total Con_media Con_cell 
# [1,] 0.00 1495.242 4.75830  1500 9.490000 0.510000 
# [2,] 0.01 1488.103 11.89710  1500 9.442408 1.275145 
# [3,] 0.02 1481.028 18.97216  1500 9.395241 2.033457 
# 
# [[2]] 
#  time Amt_media Amt_cell Amt_total Con_media Con_cell 
# [1,] 0.00 2990.483 9.51660  3000 18.98000 1.020000 
# [2,] 0.01 2976.550 23.44980  3000 18.88711 2.513377 
# [3,] 0.02 2962.739 37.26072  3000 18.79504 3.993646 
# 
# [[3]] 
#  time Amt_media Amt_cell Amt_total Con_media Con_cell 
# [1,] 0.00 4485.725 14.27490  4500 28.47000 1.53000 
# [2,] 0.01 4465.153 34.84653  4500 28.33286 3.73489 
# [3,] 0.02 4444.761 55.23920  4500 28.19690 5.92060