Python：使用odein实现阈值模型

Question

我正在处理的问题（显示的示例高度简化）似乎是一个常见的问题，但我还没有找到解决方案。我有三种不同的反应，v1，v2和v3，定义如下：

v1:R<；->；A+C；v1=k1*（R-A*C/5000）

v2:R<；->；B+C；v2=k2*（R-B*C/5000）

v3:A+B->；p；v3=k3*A*B

使用资源R，前两个反应分别生成a和C以及B和C，而第三个反应将a和B转化为产物p（k1、k2、k3是常数，这里它们被设置为1）。在

第三个反应只有在C超过一个称为Cthr的阈值（这里：Cthr=25）时才会发生，否则v3为0。所以我们的想法是C积累，一旦达到一定的浓度，就会产生产物P

我执行如下：

def thresholdmodel (yn,tvec,allpara,R):

    (A, B, C, P) = yn

    k1, k2, k3 = allpara['kv']    

    Cthr = allpara['Cthresh']

    if C <= Cthr:
        v3 = 0
    else:      
        v3 = k3*A*B
        C = 0 #does not(!) affect the ouput, why?

    v1 = k1*(R - A*C/5000.)
    v2 = k2*(R - B*C/5000.)

    dA = v1 - v3
    dB = v2 - v3
    dC = v1 + v2
    dP = v3

    return (dA, dB, dC, dP) 

模拟的输出如下：http://i50.tinypic.com/apdvkj.png

所以很明显，它一直工作到第一次达到阈值（v3为0，p不产生），但之后C没有设置为0，我不知道为什么。
我想得到的是：C产生到阈值，降到0，再次产生，降到0，以此类推，就像锯齿波。
P的时间进程应该看起来像楼梯（只有在超过Cthr时才会产生）。在

有人知道我要怎么做才能把C设置回0并得到预期的输出吗？谢谢！在

Answer 1

我不懂你的方程式，但我能告诉你为什么C没有降到0。在

因为C是thresholdmodel()中的局部变量，因此将C更改为任何值都不会更改odeint中的状态。odeint将始终集成由thresholdmodel()返回的dC。因此C将持续增加。在

编辑：

C将持续增加，因此每次C>；threshold都需要更改阈值，例如：

lastC = 0

def thresholdmodel (yn,tvec,allpara,R):
    global lastC
    (A, B, C, P) = yn

    k1, k2, k3 = allpara['kv']    

    Cthr = allpara['Cthresh']

    if C <= lastC + Cthr:
        v3 = 0
    else:      
        v3 = k3*A*B
        lastC = C

    v1 = k1*(R - A*C/5000.)
    v2 = k2*(R - B*C/5000.)

    dA = v1 - v3
    dB = v2 - v3
    dC = v1 + v2
    dP = v3

    return (dA, dB, dC, dP)