RxODEによる常微分方程式モデルに基づくPK/PD解析結果をShinyのGUIでリアルタイムに可視化する方法

はじめに

PK-PD解析に特化したRのパッケージとして，RxODEと呼ばれる微分方程式モデルを解くライブラリが存在します。これ自体はNONMEMと同様の記法でコンパートメント解析を行ったり，微分方程式を自身で記述することでQSP解析なども実施することが可能です。
こちらのライブラリの長所としては，ライセンスフリーな点や一度Cへとコンパイルしてから計算を行うため，実行速度の面で優れている，ということが挙げられます。
加えて，Shinyと呼ばれるRのWebアプリケーションとの互換性があり，作成したモデルのパラメータ調整をGUIでリアルタイムに結果を確認しつつ実行可能であるというメリットがあります。
解析者のパラメータ探索を効率化するのみならず，実際に解析結果を説明する際にリアルタイムでdose regimen等を変更できることが大きなアドバンテージになります。
GUIの実装自体はRxODEで解析用のスクリプトを作成した後，その延長線としてRのオブジェクトをShinyに渡すことでGUIでの可視化が可能になります。
その際に必要となるファイルは以下の3つで，いずれもRxODEとShinyを連結するgenShinyTemplate()でテンプレートを作成することが可能です。以下の3ファイルはすべて同じディレクトリに格納する必要があり，runApp("ディレクトリ名")でGUIを起動可能です。

1. RxODEの解析オブジェクト(.rda) : 必要となるオブジェクトはmod(モデルのデータフレーム), params(パラメータ), inits(初期値)の3つです。
   これらはRの標準機能であるsave()ファンクションで.rdaファイルとして保存することで作成可能です。特に追加の操作は必要ありません。
2. ui.R : GUIを規定するスクリプト
   shinyUI()中に，SidebarPanel()とmainPanel()を定義し，それらに調整したいパラメータのスライダー等(sliderInput,checkboxInput)を追加していくことで作成します。genShinyTemplate()でテンプレートを自動で作成できるので，目的に応じてそれぞれのパラメータを追加していく，といったイメージです。
3. sever.R : severでの動作を規定するスクリプト
   shinySeverfunction()中に，解析モデルを規定するreactive(), GUI上に表示する結果を規定するoutput()を追加していくことで作成します。reactive()自体はほとんど解析時に作成したRxODEのスクリプトを貼り付けることで作成でき，output＄summaryにrenderPrint()で表示するパラメータを，output$CpPlotにrenderPlot()で描画するプロットを格納していく，といったイメージです。

このように，通常の解析に，Server.Rとui.Rを追加で作成することでGUIwebサーバーでの可視化が可能になり，サーバーのURLやポート，GUIの微調整といった煩わしい記述を一切必要としないことが特徴です。

また，RxODEには非線形混合効果モデルによるpopPK解析を行うnlmixr2と呼ばれる拡張パッケージが存在します。
こちらでも同様にShinyによる可視化と探索，サーバー上での解析の実行が可能ですのでご確認ください。
https://qiita.com/kentaro_BIC/items/0aa9c53ac359ebcfe037

1. RxODEでの解析(.rdaファイルの作成)

RxODEでの解析は基本的に微分方程式モデルの作成(ode)/モデルのコンパイルコンパイル(RxODE)，インプットデータの作成(eventTable)，解析の実行の3ステップに分かれます。

1. 1. モデルの作成/コンパイル

モデルの作成/コンパイルは更に以下の3ステップに分類されます。

Step1. 微分方程式(ode <- "微分方程式")の定義

微分方程式の定義はオブジェクトに微分方程式を格納することから始まります。NONMEMにおける$DESに該当し，2コンパートメントモデルの場合，下記のような記述になります。

# Define 2-compartment model
ode <- "
   C2 = centr/V2;	
   C3 = peri/V3;
   d/dt(depot) = -KA*depot;
   d/dt(centr) = KA*depot - CL*C2 - Q*C2 + Q*C3;
   d/dt(peri) = Q*C2 - Q*C3;
   d/dt(eff) = Kin*(1-C2/(EC50+C2)) - Kout*eff;
"

Step2. 微分方程式のコンパイルRxODE(model=ode)

次に定義したモデルをRxODE()コマンドでコンパイルします。C言語にコンパイルすることで計算時間を向上させています。

mod1 <- RxODE(model=ode, modName="mod1")

Step3. システムパラメータの設定

次にPKパラメータと初期値を設定します。共変量や分散を定義したいときはこのステップでそれらも定義して下さい。下記の例ではMASSライブラリのmvrnorm()関数を使い，ランダム行列を作成しクリアランスとcentralコンパートメントの分布容積(V2)に指数誤差モデルを適用しています。

# Simulate 10o times
nsub <- 100

# Create a parameter matrix with correlated interindividual variability on CL and V2
library(MASS)
# initial estimates of the variance of residual errors
sigma= matrix(c(0.09,0.08,0.08,0.25),2,2)
mv = mvrnorm(n=nsub, rep(0,2), sigma)
CL = 7*exp(mv[,1])
V2 = 40*exp(mv[,2])

params.all <- cbind(KA=0.3, CL=CL, V2=V2, Q=10, V3=300, Kin=0.2, Kout=0.2, EC50=8)     
inits <- c(0,0,0,1)

1. 2. インプットテーブルの作成(eventTable())

RxODEではNONMEMのようにcsvのテーブルを読み込む方法に加え，eventTable()コマンドを用いて直接R上で投与テーブルを作成する方法があります。後者はeventTable()コマンドでオブジェクトを作成した後,それぞれadd.dosing()で投与量を，add.sampling()でサンプリング頻度を格納することが可能です。

# create an empty event table object
ev <- eventTable()

# specify dosing
ev$add.dosing(dose = 10000, nbr.doses=10, dosing.interval =12) # 5 days BID dosing

# specify sampling
ev$add.sampling(0:120) # 24*5 hours

1. 3. 解析の実行

解析の実行自体はモデル，パラメータ(params)，テーブル(ev)，初期値(inits)を持ちいて下記のように実行でき，結果はデータフレームとして出力されます。可視化したい列を選んでmatplotで特定のコンパートメントでの濃度経時推移や薬効変化などを確認できます。

# simulate
x <- mod1$run(params, ev, inits=inits)

# create the linear-log concentration-time plot in central conpartment
matplot(x[,"C2"],type="l",xlab="time(h)",ylab="Central Conc.")

シミュレーションを繰り返したいときは下記のようにforループやapply関数を用いて繰り返し上記の操作を実行してください。

# perform simulation 100 times
nsub <- 100

# create an empty matrix for storing results
res <- NULL

# Loop through each row of parameter values and simulate
for (i in 1:nsub){
    params <- params.all[i,] # get a given paramter
    x <- mod1$run(params, ev, inits=inits) #simulate
    res <- cbind(res,x[,"eff"]) # retrieve the efficacy result row
}

上記のシミュレーション結果はresに格納されているので，下記のように可視化が可能です。また，95%CIもquantile()関数を用いて計算可能です。

# Plot results of Efficacy(supression rates of a given receptor)
par(mfrow = c(1,2), mar = c(4,4,1,1))
matplot(res, type="l", ylab = "Effect", xlab="Time(hrs)"

# Calculate and plot quantiles
res.q <- apply(res,1,quantile,prob=c(0.05,0.5,0.95))# calculate 5%, mean(50%), 95%
res.q.t <- t(res.q) # transpose the result
matplot(res.q.t, type="l", lty=c(2,1,2), col=c(2,1,2), ylab="Effect", xlab=Time(h)") # plot the result
legend("topright",c("median,"5 and 95%),lty=c(1,2),col=c("black","red"),cex=0.8) # add the legend to the figure

2. Shinyテンプレートの作成

RxODEの解析はgenShinyApp()関数を用いてRshinyと連結させることが可能です。下記のコマンドで現在のディレクトリ以下にShinyで可視化するために必要なファイルとServer.R, ui.Rのテンプレートが作成されます。

library(RxODE)
library(shiny)

# generate the directory for shiny template and needed files
setwd("path/to/directory")
genShinyApp.template(appDir="Shiny",verbose=TRUE)

この操作でカレントディレクトリの中に指定した名称でフォルダが作成されます。フォルダの中身としては以下のとおりです。

1. sever.R : サーバーでの動作を規定するスクリプトテンプレート
2. ui.R : uiを規定するスクリプトテンプレート
3. rx_shiny_data.rda : RxODEのオブジェクトを格納したファイル
   こちらが存在しなかった場合，下記の手順で作成して下さい
4. mod1.d(フォルダ) : その他，Shinyの実行に必要なファイル群
   この後はそれぞれsever.R, ui.Rを書き換えて行くことでGUIの実装を行っていきます。
   また，rx_shiny_data.rdaファイルが作成されていない場合，save()関数を用いて，rdaオブジェクトを作成して下さい。この時，Rのワークスペースに重複や無関係のデータが存在するとうまく動作しない場合がありますので，その場合は改めRstudioを再起動し，ワークスペースをリセットした上でこちらの操作を行って下さい。

setwd("Shiny")
save(mod1, params, inits, file="rx_shiny_data.rda")

3. ui.Rファイルの記述

ui.Rファイルでインターフェイスの外観，調整するパラメータ，表示する結果を規程します。genShinyApp.template()ですと下記のようなテンプレートが作成されますので，必要に応じて書き換えて行って下さい。

#
# Dosing regimen template automatically generated by
# RxODE::genShinyApp.template()
#

library(RxODE)
library(shiny)

# Define UI for dataset viewer application
shinyUI(pageWithSidebar(

  # Application title.
  headerPanel("Regimen simulator"),

  # Sidebar with controls to select a dataset and specify the number
  # of observations to view. The helpText function is also used to
  # include clarifying text. Most notably, the inclusion of a
  # submitButton defers the rendering of output until the user
  # explicitly clicks the button (rather than doing it immediately
  # when inputs change). This is useful if the computations required
  # to render output are inordinately time-consuming.
  sidebarPanel(
    sliderInput("Dose", "Dose:",
                  min=100, max=1000, value=400, step=50),

    selectInput("regimen", "Regimen:",
                choices = c("QD", "BID"),
                selected="QD"),

    selectInput("cycle", "Cycle:",
                choices = c(
                        "continous",
                        "1wkon 1wkoff",
                        "2wkon 1wkoff",
                        "3wkon 1wkoff"),
                selected="2wkon 1wkoff"),

    sliderInput("ncyc", "# cycles:",
                  min=1, max=5, value=3, step=1),

    radioButtons("compartment", "Select PK compartment for plotting:",
                  choices=c("depot", "centr", "peri", "eff"),
                  selected = "depot")
  ),

  # Show a summary of the dataset and an HTML table with the requested
  # number of observations. Note the use of the h4 function to provide
  # an additional header above each output section.
  mainPanel(
    h4("State variables"),
    verbatimTextOutput("summary"),
    h4("Concentration time course"),
    plotOutput("CpPlot")
  )
))

このスクリプトの本体はshinyUI()中にあり，こちらの中身を編集することでGUIを設定していきます。今回は下図と同様，サイドバーパネルを持ったページshinyUI(pageWithSidebar())の関して説明を行っていきます。
大まかに分けてShinyUI(pageWithdiSidebar())はサイドバーの表示を設定するheaderPanel() (下図タイトル), sidebarPanel() (下図左グレー箇所参照)とmainPanel() (下図右参照)の2要素で構成されており，それぞれsidebarPanel中に調整したいパラメータを，mainPanel上に表示したいアウトプットを格納します。

概略としては以下のとおりです。

shinyUI(pageWithSidebar(
    headerPanel("タイトル"),
    sidebarPanel(<sidebarのオブジェクト(sliderInput,CheckboxInputなど)>),
    mainPanel(<メインパネルの出力(数値テキスト，図)>)
))

3.1. sidebarPanel()の設定

基本的にSidebarpanelは以下の要素で成り立っています。これらをカッコの中に","で並列して記述することでそれぞれ要素を付け足すことが可能です。

• sliderInput(inputId(入力オブジェクト),label(表記名),min,max,value(平均),step(刻み幅))
  sliderInputで左上図のようにスライダーを用いて調整するパラメータとその範囲，および刻み幅を設定することが可能です。以下の例ではそれぞれDose,日数とサンプリング頻度，改行とサブタイトルを挟み，Effect(受容体占有率)の有効域幅をパラメータとして設定しています。

shinyUI(pageWithDisebar(
    sidebarPanel(
        sliderInput("Dose", "Dose:",min=0, max=20, value=10, step=.5),
        sliderInput("nDays", "# of Days:",min=5, max=50, value=20, step=1),
        sliderInput("sampFreq","Sampling Frequency (days):", min=1,max=10,value=1,step=1),
        br(), #空白行
        h4("Define limits for target range:"),
        sliderInput("lowerEffectLimit","Lower Effect Limit:",min = 0, max=0.6,value=0.4,step=0.05),
        sliderInput("upperEffectLimit","Upper Effect Limit:",min = 0.4, max=1,value=0.6,step=0.05),
    	br(),    

• br() : 空白行, h4("文字列") : 水平文字列
  オブジェクトとオブジェクトの間に空白行や文字をはさみたいときには上記のようにそれぞれbr(),h4("文字列")などと追加することができます
• checkboxInput(inputId, label, value=TRUE/FALSE(初期設定))
  また，スライダーだけでなく，チェックボックスでの条件設定も同様に可能です。刻み幅でなく，value=TRUE/FALSEを設定することで，それぞれの条件を変更できます。下記の例では有効域(40%<Eff<60%)中に収まるように投与量を調節する際にトラフで調節するかピークで調節するか(初期値=トラフ)を設定しています。

        checkboxinput("troughcheck", label="Control trough effect (24hrs) within target range", value=TRUE),
        checkboxinput("peakcheck", label="Control peak effect (tmax~12hrs) within target range")
        submitButton("Update View")

• submitButton("Update View"), その他
  こちらのボタンは上記のcheckboxinputの変更に伴い，図表をアップデートするボタンを追加する関数になります。他にも下記のような関数がありますので，詳細はこちらのShinyのwebページ(https://shiny.posit.co/r/reference/shiny/1.3.1/submitbutton) からご確認下さい。
  actionButton, checkboxGroupInput, dateInput, dateRangeInput, fileInput, numericInput, passwordInput, radioButtons, selectInput, textAreaInput, textInput, varSelectInput

3.2. mainPanel()の設定

mainPanelは基本的にsummaryパラメータを出力するverbatimTextOutput()と，図を出力するplotOutput()の２要素からなっています。出力要素(outputId)としてserver.ui側で規程したオブジェクトを出力することが可能です。なお，br()やh4()といった要素はsidebarPanel()と共通です。

• verbatimTextOutput(outputId)
  こちらは実際の解析で計算されたサマリーパラメータを文字列として表示することができます。下記の例では有効域(40%<Eff<60%)から外れている時間の割合を格納した"summary"オブジェクトを表示しています。
• plotOutput(outputId, width, height, etc..)
  こちらの関数ではserver.uiで描画したplotオブジェクト(R base graphics)を表示することができます。ggplot2やlatticeなどのgrid-based graphicsと呼ばれるライブラリは使用できないので，server.Rでの設定時に注意が必要です。
  下記の例では，縦軸に有効性,横軸に時間を取ったプロット("CpPlot")を出力しています。

    mainPanel(
        h4("Summary of therapeutic monitoring:")
        verbatimTextOutput("summary"),
        br(),
        h4("Concentration time course"),
        plotOutput("CpPlot")

4. server.Rファイルの設定

ui.Rファイルでインターフェイスの外観，調整するパラメータ，表示する結果を規程します。genShinyApp.template()ですと下記のようなテンプレートが作成されますので，必要に応じて書き換えて行って下さい。

#
# Dosing regimen template generated by RxODE::genShinyApp.template()
#

debug = TRUE
#wd = sprintf("%s/../", getwd())
#setwd(wd)

# Server inputs: Dose, dosing regimen, dosing frequency,
# dosing cycle definition, number of dosing cycles

library(shiny)
library(RxODE)

# read objects from "rx_shiny_data.rda" in the  AppDir folder,
# objects include, mod1, params, inits, method, atol, rtol.]

load("./rx_shiny_data.rda")
if (!rxDynLoad(mod1)) mod1 <- RxODE(mod1, modName="mod1")
# Define server logic
shinyServer(function(input, output) {

  get.cp <- reactive({
    ds <- input$Dose
    reg <- switch(input$regimen, "QD"=1, "BID"=2)
    cyc <- switch(input$cycle,
        "continous"=c(7,0),
        "1wkon 1wkoff"=c(7,7),
        "2wkon 1wkoff"=c(14,7),
        "3wkon 1wkoff"=c(21,7)
    )
    cyc <- rep(1:0, cyc)
    ncyc <- input$ncyc
    lcyc <- length(cyc)

    ev <- eventTable()
    for (i in 1:ncyc) ev$add.dosing(
        dose=ds,
        nbr.doses=sum(cyc)*reg,
        dosing.interval=24/reg,
        start.time=(i-1)*lcyc*24
    )
    ev$get.EventTable()
    ev$add.sampling(0:(ncyc*lcyc*24))

    mod1$solve(params, ev, inits, method=method, atol=atol, rtol=rtol)
  })


  output$summary <- renderPrint({
    x <- get.cp()
    print(x[1:4,])
    if (debug) print(getwd())
  })

  output$CpPlot <- renderPlot({
    x <- get.cp()
    cmp <- input$compartment
    plot(x[,c("time", cmp)], xlab = "Time", ylab = "Drug amount",
          main = cmp, type = "l")
  })
})

こちらのスクリプトはload("./rx_shiny_data.rda")の保存したRxODEオブジェクトを読み込む部分とサーバー側の動作を規定するshinyServer(function(input,ouput){})から構成されています。上部のコメントアウト箇所はディレクトリを移動するコマンドなので，runApp()コマンドを正しいディレクトリで実施する場合には特に気にしなくて大丈夫です。
なお，ShinyServer()の概略としては以下の通りです。それぞれ，入力データの処理を行うreactive()箇所，summary paramterの出力を行うrenderPrint(),オブジェクトの描画と出力を行うrenderPlot()の3要素に分けられています。output$~以下がui.Rでの出力時に指定するoutputIdになります。

shinyServer(function(input, output) {
    get.cp <- reactive({}) #save results as get.cp
    output$summary <- renderPrint({}) #store summary parameters in output$summary
    output$CpPlot <- renderPlot({})

4.1. サーバー側の処理:reactive({})の設定

こちらの{}の中に実際にServer上で行う解析内容を記載します。基本的にはRxODEの実行時のスクリプトをそのまま記載できるため，"rx_shiny_data.rda"に保存したObjectから必要なデータをinput$データ名，として読み込み，同様に解析を行って下さい。解析結果を格納したオブジェクトはreturn(オブジェクト名)で出力することができ，下記の場合get.cpオブジェクトに格納され，その後のrenderPrint()やrenderPlot()で処理することが可能です。
以下が#1で行った解析での例になります。

shinyServer(function(input, output) {
  get.cp <- reactive({

  #define decision Rule, using inputs from UI
	effect.limits <- c(0, input$lowerEffectLimit, input$upperEffectLimit, 9)  
	dose.multipliers <- c(.5, 1, 2)   #decision rule effects  

	#Define simulation parameters, using inputs from UI
	nDays <- input$nDays				   #number of days to simulate
	unit.dose <- input$Dose*1000			   #starting dose	   
	start.dose <- 1
	sampling.frequency <- input$sampFreq   	   #sample every day

	#Initialize other variables
	vars <- c("depot", "centr", "peri", "eff")
	res <- NULL 		#results vector
	timevec <- NULL	#time vector
	doses <- NULL	#doses vector

	#Simulate each day. At the end of each day, test the effect  level, and adjust the dose level according to the decision rule
	for (i in seq(1,nDays,by=sampling.frequency)) {
		if (i==1)   #initialize on first day
		{
			inits <- c(0, 0, 0, 1)      
			last.multiplier <- start.dose
			this.multiplier <- 1
		} else	#use end of previous day as initial conditions for next day
		{
			inits <- x[dim(x)[1], vars]    							#use last value of state variables as new initial conditions
			if (input$troughcheck==TRUE) {
				wh <- cut(inits["eff"], effect.limits)					#compare trough effect with decision rule limits
				this.multiplier <- dose.multipliers[wh]					#determine dose multiplier accordingly
			}
			if (x[13,"eff"]<effect.limits[2] & input$peakcheck==TRUE) {
				this.multiplier <- dose.multipliers[1]
			}
		}
	#adjust dose and store new dose		
	this.multiplier <- this.multiplier*last.multiplier	
	last.multiplier <- this.multiplier			

	#specify dosing and sampling	
	ev <- eventTable()						
	ev$add.dosing(dose=this.multiplier*unit.dose, dosing.interval=24, nbr.doses=sampling.frequency)	
	ev$add.sampling(1:(24*sampling.frequency))	

	#simulate			
	x <- mod1$solve(params, ev, inits)			
	
	#compile output matrix
	timeTotal <- ev$get.EventTable()[,"time"]+(i-1)*(24)
	doses <- rep(last.multiplier,length(timeTotal))
	x <- cbind(x,timeTotal,doses)
	res <- rbind(res,x)
}
	return(res)

上記の例ではパラメータの設定とeventTable()の作成，有効域に基づいた投与量の調節(Eff > 60% : 用量半減， Eff < 40% : 用量倍増), シミュレーションの実施(mod1$solve), アウトプット行列のコンパイル(res <- rbind(res,x)), 出力(return(res))を実施し，時間，投与量，有効性からなる出力結果(res)をget.cpオブジェクトに格納しています。

4. 2. ui.Rに返すoutputの設定

上記のreaction({})で出力したget.cp()をもとに，それぞれ出力するsummaryパラメータとプロットオブジェクトを作成していきます。
outputオブジェクトが初期値として与えられており，これに対してデータフレームのようにoutput$と格納していくことでui.RからGUIに出力することが可能です。
以下が概要になります。

    # generate summary of output to be displayed in GUI window
    output$summary <- renderPrint({})
    
    # generate summary of figure py plot{} function
    output$CpPlot <- renderPlot({})

上記の解析例ですとそれぞれの記載は以下の通りになります。

    # generate summary of output as string
    output$summary <- renderPrint({
        res = get.cp() #store the returned object from reactive({})
        lo = mean(res[,"eff"]<input$lowerEffectLimit)*100
        hi = mean(res[,"eff"]>input$upperEffectLimit)*100
        cat(sprintf("%% time below lower limit: %5.1f%%\n", lo))
        cat(sprintf("%% time above upper limit: %5.1f%%\n", hi)) # output the .str result with cat() function
  })
    
    # generate summary figure as plot() object
    output$CpPlot <- renderPlot({
        res <- get.cp()
        plot(res[,"timeTotal"]/24, res[,"eff"],type="l",ylim=c(0,1),xlim=c(0,input$nDays),xlab="Time (days)", ylab="Effect / Dose Multiplier")
        rect(-200, input$lowerEffectLimit, 1500, input$upperEffectLimit, col=rgb(0.5,0.5,0.5,1/4)) # add the effective area as a gray box
        points(res[,"timeTotal"]/24, res[,"doses"], type="s",xlab="Time (days)", ylab="Dose",col="red",lty=2) # add the points
        legend("topright",c("Effefct","Dose Multiplier"),lty=c(1,2), col=c("black","red")
    })
})

5. runApp()によるShinyの実行

以上の手順でそれぞれrx_shiny_data.rda, server.R, ui.R，その他必要ファイルを作成した後は当該のディレクトリでのRコンソール上でrunApp()コマンドを実施することでGUIが開かれ，ui.Rで規程したパラメータをスライダ-やコマンドボックスで調整しリアルタイムで描画することが可能です。

setwd("path/to/dir")
runApp("Shiny")

参考文献

1. A Tutorial on RxODE: Simulating Differential Equation Pharmacometric Models in R https://ascpt.onlinelibrary.wiley.com/doi/full/10.1002/psp4.12404
2. Shiny Basics https://shiny.posit.co/r/getstarted/shiny-basics/lesson1/index.html
3. nlmixr2 https://qiita.com/kentaro_BIC/items/0aa9c53ac359ebcfe037