18_葡萄糖模型

18 葡萄糖模型

在这章节，我们实现了前一章描述的葡萄糖最小模型。我们将从run_simulation开始，用离散时间步长解微分方程。这个方法可以很好地适用很多应用，但是不够准确。在这章我们将探索一个更好的选择—— ODE solver

18.1 实现

在开始之前，我们假设这个模型的参数都是已知的。我们将实现这个模型并使用它生成两个时间序列G和X。然后我们将学到怎样找到生成最符合数据的序列参数。

我们可以先用提前预估的值来假定这些参数：

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

Params对象和System或者State对象类似，就是用来装载变量的集合。

我们可以将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时的状态值为基准，即Gb和Ib。t_0和t_end的值从数据中获得。参数G0作为G的初始值。然后我们就可以把这些放进对象System中。

下面是更新函数：

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)

通常来说，更新函数会以State对象、时间和System对象作为参数。第一行使用了可变参数来获取G和X的瞬时值。

之后的几行把我们需要的参数从System对象中提取出来。

直接计算导数dGdt和dXdt；我们仅仅只是把数学方程式翻译成了Python语句。

然后我们用每一个微分乘以时间间隔dt来运行更新函数，dt在我们这个例子中是2分钟。返回值是State对象和新的G和X的值。

在运行模拟之前，先用初始条件运行一下更新函数更好一些：

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

如果运行没有报错而且结果也是正确的我们就可以开始运行模拟了。我们将使用和之前极为相似的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

我们可以这样运行：

results = run_simulation(system, update_func);

图18.1展示了结果。上面的图显示的是该模型模拟的葡萄糖水平以及测量的数据。下面的图显示的是模拟的组织液中胰岛素水平，但是单位不详，不要和血液中测量的胰岛素水平混淆。

18_1.png

图18.1：葡萄糖最小模型运行结果

我选择的参数值可以很好的拟合真实数据，除了开始的极少数不需要特别精准的数据点。

18_2.png

18.2 解微分方程

到目前为止我们已经通过把它改写成差分方程来解出了微分方程。在现在的例子中，微分方程是：
$\frac{dG}{dt}=-k_1[G(t)-G_b]-X(t)G(t)$

$\frac{dX}{dt}=k_3[I(t)-I_b]-k_2X(t)$

两边同乘以dt：
$dG=[-k_1[G(t)-G_b]-X(t)G(t)]dt$

$dX=[k_3[I(t)-I_b]-k_2X(t)]dt$

当dt非常小的时候或者无穷小的时候，方程就是精确的。但是在我们的模拟中dt是2分钟，并不是很小。我们的模拟中假设导数dG/dt和dX/dt在每个2分钟的时间间隔中都是常数。

这种在假设导数在离散的时间间隔中是常数的计算导数的方式，我们称之为欧拉法（详见http://modsimpy.com/euler）。

欧拉法可以很好地适用于一些简单的问题，但是并不是很准确。其他的一些方法更精确，但是大多数都复杂得多。

其中一个又简单又好的方法叫Ralston法。ModSim库也提供了一个run_ode_solver函数来实现这种方法。

run_ode_solver中的”ODE“代表”ordinary differential equation"即“一阶微分方程”。我们解决的方程式是一阶的因为所有的导数都是关于同一变量；换句话说就是没有偏导数。

在运行run_ode_solver之前我们需要构建一个“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类似；实际上，它以相同的顺序包含相同的参数。但是slope_func更简单，因为我们仅仅只需计算导数，也就是“斜率”。我们不需要构建更新函数，因为run_ode_solver中已经包含了这些步骤。

我们可以这样运行run_ode_solver:

results2, details = run_ode_solver(system, slope_func)

图18.2：欧拉法和拉斯顿法的结果

run_ode_solver和run_simulation类似：它把System对象和斜率函数作为参数，返回两个值：一个包含结果的TimeFrame和一个包含附加信息的ModSimSeries。

ModSimSeries和System或者State类似；是一组变量及其值的集合。run_ode_solver中的ModSimSeries(我们将其分配给details)包含了有关求解器如何运行的信息，包括成功代码以及诊断结果。

我们分配给results的TimeFrame，行是代表时间间隔，列是代表每个状态变量。在这个例子中，行是从0到182分钟的时间间隔，列是两个状态变量，G和X。

图18.2展示了run_ode_solver和run_simulation的结果，几乎看不到两者的差别。

我们可以计算一下两者的百分差异值：

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

最大的百分差异值小于2%，这已经足够小以至于我们可以忽略。在之后的案例中我们将看到拉斯顿法比欧拉法更准确的例子。

在继续之前，你可以阅读一下本章chap18.ipynb的笔记本并完成练习。下载和运行的代码介绍见0.4节。

本书的中文翻译由南开大学医学院智能医学工程专业2018级、2019级的师生完成，方便后续学生学习《Python仿真建模》课程。翻译人员（排名不分前后）：薛淏源、金钰、张雯、张莹睿、赵子雨、李翀、慕振墺、许靖云、李文硕、尹瀛寰、沈纪辰、迪力木拉、樊旭波、商嘉文、赵旭、连煦、杨永新、樊一诺、刘志鑫、彭子豪、马碧婷、吴晓玲、常智星、陈俊帆、高胜寒、韩志恒、刘天翔、张艺潇、刘畅。

整理校订由刘畅完成，如果您发现有翻译不当或者错误，请邮件联系changliu@nankai.edu.cn。

本书中文版不用于商业用途，供大家自由使用。

未经允许，请勿转载。

18_葡萄糖模型