16.5 미분 방정식 형태로 주어진 모형의 적합
미분 형태로 주어진 모형에서 수치적분(numerical integration)을 하면서 목적함수의 필수 요소인 \(\frac{\partial F}{\partial \eta_i}\) 를 구하는 것은 상당히 번거롭다. 하지만, Leibniz integral rule이라고도 불리는 differentiaion under the integral sign을 활용하면 시간이 상당히 단축된다.
즉, 적분 기호 밖의 미분 기호가 적분 기호 안으로 들어갈 때 편미분 기호로 바뀌어야 하는데, 의학에서의 문제들은 다행히 적분 구간(\(t_1\)에서 \(t_2\))이 적분 변수(time t)에 무관한 것이어서 상대적으로 단순한 형태이다.
좀 더 구체적으로 설명하면 다음과 같다.
관찰 구획(Observation compartment)에서 약물(drug, substance)의 양(amount)을 A라 하고 그 농도를 F라 하면(즉, \(A = F \cdot V\)) 다음과 같은 관계를 유도할 수 있다.
\[\frac{\partial F}{\partial \eta_i} = \left( \frac{dA}{d\eta_i} - F\frac{\partial V}{\partial \eta_i} \right) / V\]
여기에서 \(\frac{dA}{d\eta_i}\) 를 얻는 것이 일견 쉽지 않아 보이는데, 이것은 주어진 연립 미분 방정식(coupled differential equation)을 먼저 \(\eta\)들로 편미분(partial differentiation)한 뒤에, 이를 다른 \(\frac{dA}{dt}\) 들과 함께 시간에 대해 수치적분(numerical integration)을 하면 된다.
저장 구획(Depot compartment)의 양(amount)을 \(A_1\)이라 하고, 중심 구획(central compartment)의 양(amount)을 \(A_2\)라 하면 다음과 같은 미분 방정식을 세울 수 있다.
\[\begin{align*} \frac{d A_1}{dt} &= -k_a A_1 \\ \frac{d A_2}{dt} &= k_a A_1 - k A_2 \\ \end{align*}\]
\(A_2\)부분에 대한 것만 나타내면 다음과 같다.
\[\begin{align*} \frac{d A_2}{d \eta_i} &= \frac{d}{d \eta_i} \int \frac{d A_2}{dt} dt = \int \frac{d}{dt} \left( \frac{\partial A_2}{\partial \eta_i} \right) dt \\ \frac{d}{dt} \left( \frac{\partial A_2}{\partial \eta_i} \right) &= \frac{\partial k_a}{\partial \eta_i} A_1 + k_a \frac{\partial A_1}{\partial \eta_i} - \frac{\partial k}{\partial \eta_i} A_2 - k \frac{\partial A_2}{\partial \eta_i} \\ \end{align*}\]
\(\frac{\partial A_1}{\partial \eta_i}\)를 \(B_{1i}\)라 하고, \(\frac{\partial A_2}{\partial \eta_i}\)를 \(B_{2i}\)라 하면 다음과 같이 식을 단순화 할 수 있다.
\[\frac{dB_{2i}}{dt} = \frac{\partial k_a}{\partial \eta_i} A_1 + k_a B_{1i} - \frac{\partial k}{\partial \eta_i} A_2 - k B_{2i}\]
위에서 \(\frac{\partial k_a}{\partial \eta_i}\), \(\frac{\partial k}{\partial \eta_i}\) 등은 다음의 식을 편미분(partial differentiation)해서 구할 수 있다.
\[\begin{align*} k_a &= \theta_1 e^{\eta_1} \\ V &= \theta_2 e^{\eta_2} \\ k &= \theta_3 e^{\eta_3} \\ \end{align*}\]
만약 구획(compartment)이 2개이고, \(\eta\)가 3개라면 위와 같은 식이 6개가 나온다. 이것을 원래 주어진 식 2개와 함께, 즉 8개를 한꺼번에 수치적분하면 원하는 것(특정 시간에서 각 구획의 약물 양과 그것의 \(\eta\)에 대한 derivative)을 빠르게 구할 수 있다.
아래는 일구획 모헝(one-compartment model)이 미분 방정식으로 주어졌을 때, 시간(time) t에서의 양(amount) (\(A_i\)) 뿐 아니라, \(\frac{dA}{d\eta_i}\)도 함께 구하기 위한 수치적분(numerical integration) 예이다. 아래에서 앞의 두 개가 주어진 미분 방정식(differential equation)이며, 뒤의 6개는 앞의 2개식을 3개의 \(\eta\)로 편미분(partial differentiation)해서 나온 방정식이다.
PKo1c = function(t, y, p)
{
dy1dt = -p["ka"] * y[1]
dy2dt = p["ka"] * y[1] - p["k"] * y[2]
dy3dt = -p["ka"] * y[1] - p["ka"] * y[3]
dy4dt = -p["ka"] * y[4]
dy5dt = -p["ka"] * y[5]
dy6dt = p["ka"] * y[1] + p["ka"] * y[3] - p["k"] * y[6]
dy7dt = -p["k"] * y[7]
dy8dt = -p["k"] * y[2] - p["k"] * y[8]
return(list(c(dy1dt, dy2dt, dy3dt, dy4dt, dy5dt, dy6dt, dy7dt, dy8dt)))
}
V = EstRes[[4]][2]
Para = c(ka = EstRes[[4]][1], k = EstRes[[4]][3])
Times = DATA[DATA$ID == 1, "TIME"]
InitVal = c(320, 0, 0, 0, 0, 0, 0, 0)
require(deSolve) #require(odesolve)
PRED = lsoda(y=InitVal, times=Times, func=PKo1c, parms=Para)
PRED = cbind(PRED, F=PRED[,"2"]/V) # F = A2/V
PRED = cbind(PRED,
G1 = (PRED[,"6"] - PRED[,"F"]*0)/V, # round F/ round eta1
G2 = (PRED[,"7"] - PRED[,"F"]*V)/V, # round F/ round eta2
G3 = (PRED[,"8"] - PRED[,"F"]*0)/V, # round F/ round eta3
H1 = PRED[,"F"], # round Y/ round epsilon1
H2 = 1) # round Y/ round epsilon2
PRED[, c("F", "G1", "G2", "G3", "H1", "H2")] # compare with TabStep() result F G1 G2 G3 H1 H2
[1,] 0.000 0.0000 0.000 0.0000 0.000 1
[2,] 4.511 2.9495 -4.511 -0.0667 4.511 1
[3,] 6.729 2.3362 -6.729 -0.2572 6.729 1
[4,] 7.444 0.6273 -7.444 -0.6497 7.444 1
[5,] 6.984 -0.1475 -6.984 -1.2453 6.984 1
[6,] 5.793 -0.1979 -5.793 -2.1254 5.793 1
[7,] 5.064 -0.1735 -5.064 -2.5389 5.064 1
[8,] 4.135 -0.1417 -4.135 -2.9112 4.135 1
[9,] 3.345 -0.1146 -3.345 -3.0643 3.345 1
[10,] 2.423 -0.0830 -2.423 -3.0010 2.423 1
[11,] 0.669 -0.0229 -0.669 -1.6901 0.669 1sdtab[sdtab$ID == 1, c("PRED", "G1", "G2", "G3", "H1", "H2")] PRED G1 G2 G3 H1 H2
1 0.000 0.0000 0.000 0.0000 0.000 1
2 4.511 2.9495 -4.511 -0.0667 4.511 1
3 6.729 2.3362 -6.729 -0.2572 6.729 1
4 7.444 0.6273 -7.444 -0.6497 7.444 1
5 6.984 -0.1475 -6.984 -1.2453 6.984 1
6 5.793 -0.1979 -5.793 -2.1254 5.793 1
7 5.064 -0.1735 -5.064 -2.5389 5.064 1
8 4.135 -0.1417 -4.135 -2.9112 4.135 1
9 3.345 -0.1146 -3.345 -3.0643 3.345 1
10 2.423 -0.0830 -2.423 -3.0010 2.423 1
11 0.669 -0.0229 -0.669 -1.6901 0.669 1위에서 \(\frac{\partial V}{\partial \eta_1} = 0\), \(\frac{\partial V}{\partial \eta_2} = V\), \(\frac{\partial V}{\partial \eta_3} = 0\)을 대입하였다.
수치 결과의 자리수를 늘리면 끝자리가 약간씩 다른데, 이는 분석적 해(analytic solution)와 수치적 해(numerical solution)의 계산 과정이 다르기 때문이다.