18 葡萄糖模型

在這章節(jié)，我們實(shí)現(xiàn)了前一章描述的葡萄糖最小模型。我們將從run_simulation開始，用離散時(shí)間步長(zhǎng)解微分方程。這個(gè)方法可以很好地適用很多應(yīng)用，但是不夠準(zhǔn)確。在這章我們將探索一個(gè)更好的選擇—— ODE solver

18.1 實(shí)現(xiàn)

在開始之前，我們假設(shè)這個(gè)模型的參數(shù)都是已知的。我們將實(shí)現(xiàn)這個(gè)模型并使用它生成兩個(gè)時(shí)間序列G和X。然后我們將學(xué)到怎樣找到生成最符合數(shù)據(jù)的序列參數(shù)。

我們可以先用提前預(yù)估的值來(lái)假定這些參數(shù)：

params = Params(G0 = 290,
                k1 = 0.03,
                k2 = 0.02,
                k3 = 1e-05)

Params對(duì)象和System或者State對(duì)象類似，就是用來(lái)裝載變量的集合。

我們可以將params和data傳給make_system:

def make_system(params, data):
    G0, k1, k2, k3 = params
    
    Gb = data.glucose[0]
    Ib = data.insulin[0]
    I = interpolate(data.insulin)
    
    t_0 = get_first_label(data)
    t_end = get_last_label(data)

    init = State(G=G0, X=0)
    
    return System(params,
                  init=init, Gb=Gb, Ib=Ib, I=I,
                  t_0=t_0, t_end=t_end, dt=2)

make_system使用t=0時(shí)的狀態(tài)值為基準(zhǔn)，即Gb和Ib。t_0和t_end的值從數(shù)據(jù)中獲得。參數(shù)G0作為G的初始值。然后我們就可以把這些放進(jìn)對(duì)象System中。

下面是更新函數(shù)：

def update_func(state, t, system):
    G, X = state
    k1, k2, k3 = system.k1, system.k2, system.k3 
    I, Ib, Gb = system.I, system.Ib, system.Gb
    dt = system.dt
        
    dGdt = -k1 * (G - Gb) - X*G
    dXdt = k3 * (I(t) - Ib) - k2 * X
    
    G += dGdt * dt
    X += dXdt * dt

    return State(G=G, X=X)

通常來(lái)說(shuō)，更新函數(shù)會(huì)以State對(duì)象、時(shí)間和System對(duì)象作為參數(shù)。第一行使用了可變參數(shù)來(lái)獲取G和X的瞬時(shí)值。

之后的幾行把我們需要的參數(shù)從System對(duì)象中提取出來(lái)。

直接計(jì)算導(dǎo)數(shù)dGdt和dXdt；我們僅僅只是把數(shù)學(xué)方程式翻譯成了Python語(yǔ)句。

然后我們用每一個(gè)微分乘以時(shí)間間隔dt來(lái)運(yùn)行更新函數(shù)，dt在我們這個(gè)例子中是2分鐘。返回值是State對(duì)象和新的G和X的值。

在運(yùn)行模擬之前，先用初始條件運(yùn)行一下更新函數(shù)更好一些：

update_func(system.init, system.t_0, system)

如果運(yùn)行沒(méi)有報(bào)錯(cuò)而且結(jié)果也是正確的我們就可以開始運(yùn)行模擬了。我們將使用和之前極為相似的run_simulation的版本：

def run_simulation(system, update_func):
    init = system.init
    t_0, t_end, dt = system.t_0, system.t_end, system.dt
    
    frame = TimeFrame(columns=init.index)
    frame.row[t_0] = init
    ts = linrange(t_0, t_end, dt)
    
    for t in ts:
        frame.row[t+dt] = update_func(frame.row[t], t, system)
    
    return frame

我們可以這樣運(yùn)行：

results = run_simulation(system, update_func);

圖18.1展示了結(jié)果。上面的圖顯示的是該模型模擬的葡萄糖水平以及測(cè)量的數(shù)據(jù)。下面的圖顯示的是模擬的組織液中胰島素水平，但是單位不詳，不要和血液中測(cè)量的胰島素水平混淆。

18_1.png

圖18.1：葡萄糖最小模型運(yùn)行結(jié)果

我選擇的參數(shù)值可以很好的擬合真實(shí)數(shù)據(jù)，除了開始的極少數(shù)不需要特別精準(zhǔn)的數(shù)據(jù)點(diǎn)。

18_2.png

18.2 解微分方程

到目前為止我們已經(jīng)通過(guò)把它改寫成差分方程來(lái)解出了微分方程。在現(xiàn)在的例子中，微分方程是：

兩邊同乘以dt：

當(dāng)dt非常小的時(shí)候或者無(wú)窮小的時(shí)候，方程就是精確的。但是在我們的模擬中dt是2分鐘，并不是很小。我們的模擬中假設(shè)導(dǎo)數(shù)dG/dt和dX/dt在每個(gè)2分鐘的時(shí)間間隔中都是常數(shù)。

這種在假設(shè)導(dǎo)數(shù)在離散的時(shí)間間隔中是常數(shù)的計(jì)算導(dǎo)數(shù)的方式，我們稱之為歐拉法（詳見http://modsimpy.com/euler）。

歐拉法可以很好地適用于一些簡(jiǎn)單的問(wèn)題，但是并不是很準(zhǔn)確。其他的一些方法更精確，但是大多數(shù)都復(fù)雜得多。

其中一個(gè)又簡(jiǎn)單又好的方法叫Ralston法。ModSim庫(kù)也提供了一個(gè)run_ode_solver函數(shù)來(lái)實(shí)現(xiàn)這種方法。

run_ode_solver中的”O(jiān)DE“代表”ordinary differential equation"即“一階微分方程”。我們解決的方程式是一階的因?yàn)樗械膶?dǎo)數(shù)都是關(guān)于同一變量；換句話說(shuō)就是沒(méi)有偏導(dǎo)數(shù)。

在運(yùn)行run_ode_solver之前我們需要構(gòu)建一個(gè)“slope function”，就像這樣：

def slope_func(state, t, system):
    G, X = state
    k1, k2, k3 = system.k1, system.k2, system.k3 
    I, Ib, Gb = system.I, system.Ib, system.Gb
    
    dGdt = -k1 * (G - Gb) - X*G
    dXdt = k3 * (I(t) - Ib) - k2 * X
    
    return dGdt, dXdt

slope_func和update_func類似；實(shí)際上，它以相同的順序包含相同的參數(shù)。但是slope_func更簡(jiǎn)單，因?yàn)槲覀儍H僅只需計(jì)算導(dǎo)數(shù)，也就是“斜率”。我們不需要構(gòu)建更新函數(shù)，因?yàn)?code>run_ode_solver中已經(jīng)包含了這些步驟。

我們可以這樣運(yùn)行run_ode_solver:

results2, details = run_ode_solver(system, slope_func)

圖18.2：歐拉法和拉斯頓法的結(jié)果

run_ode_solver和run_simulation類似：它把System對(duì)象和斜率函數(shù)作為參數(shù)，返回兩個(gè)值：一個(gè)包含結(jié)果的TimeFrame和一個(gè)包含附加信息的ModSimSeries。

ModSimSeries和System或者State類似；是一組變量及其值的集合。run_ode_solver中的ModSimSeries(我們將其分配給details)包含了有關(guān)求解器如何運(yùn)行的信息，包括成功代碼以及診斷結(jié)果。

我們分配給results的TimeFrame，行是代表時(shí)間間隔，列是代表每個(gè)狀態(tài)變量。在這個(gè)例子中，行是從0到182分鐘的時(shí)間間隔，列是兩個(gè)狀態(tài)變量，G和X。

圖18.2展示了run_ode_solver和run_simulation的結(jié)果，幾乎看不到兩者的差別。

我們可以計(jì)算一下兩者的百分差異值：

diff = results.G - results2.G
percent_diff = diff / results2.G * 100

最大的百分差異值小于2%，這已經(jīng)足夠小以至于我們可以忽略。在之后的案例中我們將看到拉斯頓法比歐拉法更準(zhǔn)確的例子。

在繼續(xù)之前，你可以閱讀一下本章chap18.ipynb的筆記本并完成練習(xí)。下載和運(yùn)行的代碼介紹見0.4節(jié)。

本書的中文翻譯由南開大學(xué)醫(yī)學(xué)院智能醫(yī)學(xué)工程專業(yè)2018級(jí)、2019級(jí)的師生完成，方便后續(xù)學(xué)生學(xué)習(xí)《Python仿真建模》課程。翻譯人員（排名不分前后）：薛淏源、金鈺、張?chǎng)?、張瑩睿、趙子雨、李翀、慕振墺、許靖云、李文碩、尹瀛寰、沈紀(jì)辰、迪力木拉、樊旭波、商嘉文、趙旭、連煦、楊永新、樊一諾、劉志鑫、彭子豪、馬碧婷、吳曉玲、常智星、陳俊帆、高勝寒、韓志恒、劉天翔、張藝瀟、劉暢。

整理校訂由劉暢完成，如果您發(fā)現(xiàn)有翻譯不當(dāng)或者錯(cuò)誤，請(qǐng)郵件聯(lián)系changliu@nankai.edu.cn。

本書中文版不用于商業(yè)用途，供大家自由使用。

未經(jīng)允許，請(qǐng)勿轉(zhuǎn)載。