在 Python 中求解耦合常微分方程 - 振荡太多/太小

问题描述 投票：0回答：1

所以我一直在尝试对我在论文中找到的耦合一阶微分方程组进行建模。它们具有生物性质，因此负面解决方案没有任何意义。

因为我知道论文中的最终图形应该是什么样子（见下文），所以我知道我的解决方案是完全错误的。

纸上图：

特别是，到目前为止，我已经用我的代码成功创建了两个场景：

  1. 默认求解器：no 和 np 剧烈振荡，而 nc 和 ns 保持恒定 - 不应该是这种情况。
  2. BDF 求解器：0 秒后立即进入“稳态”（但不是系统的实际稳态）。基本没什么用。

这是我的代码（场景 1），我决定只使用四个 DE 中的三个来为 ntitin 应保持恒定的约束腾出空间。

import numpy as nump

# Constants
kp = 0.07   # M^-1 s^-1
kr = 6  # s^-1
ks = (10**(-8)) - (10**(-6))  # s^-1
kdns = 0.002  # s^-1
ATP = 7.2*10**(-3)  # M
w0 = 10**6  # s^-1
kb = 1.3806452*10**(-23)  # J/K
T = 310  # K
u0 = 1*10**(-9)  # m
umax = 28*10**(-9)  # m
delta_G = 35*kb*T
ntitin = 4.7 * 10**(-6)  # in molar (M)  

# Differential equations
def equations(t, y):
    nc, no, np = y
    
    # Calculate ns using the constraint
    ns = ntitin - (no + nc + np)
    
    # Calculate Eplus and Eminus for the current force value
    f = force_input(t)  # Call the previously defined force input function
    Eplus = delta_G - f * u0
    Eminus = f * (umax - u0)
    
    # Calculate rate constants for the current force value
    kplus = w0 * nump.exp((-Eplus) / (kb * T))
    kminus = w0 * nump.exp((-Eminus) / (kb * T))
    
    # System of differential equations (without dns_dt since ns given through constraint )
    dnc_dt = -kplus * nc + kminus * no
    dno_dt = kplus * nc - kminus * no - kp * no * ATP + kr * np + kdns * ns
    dnp_dt = kp * no * ATP - kr * np - ks * np
    #dns_dt = ks * np - kdns * ns

    return [dnc_dt, dno_dt, dnp_dt]

# Initial conditions
no0 = 0.0006*ntitin
np0 = 0.0025*ntitin
ns0 = 10e-12
nc0 = ntitin-no0-np0-ns0 #assuming ns starts at 10e-12 and ntitin needs to stay constant
y0 = [nc0, no0, np0]
# IF NS STARTS AT 10E-6, ALL OF THE SUDDEN no INCREASES DRAMATICALLY OVER THE TIME PERIOD
#, BUT NC, NP AND NS STAY STEADY

# Time span for the solver
t_span = [-1, 3 * 260]  # 3 full cycles of the input
t_values = nump.linspace(t_span[0], t_span[1], 1000)

# Solve the system
solution = solve_ivp(
    equations, t_span, y0, t_eval = t_values)

# Calculate ns for each time point using the constraint
ns_values = ntitin - (solution.y[0] + solution.y[1] + solution.y[2])
sum_values = (ns_values + solution.y[1] + solution.y[1])/ntitin

# Plotting the results
#plt.plot(solution.t, solution.y[0]/ntitin, label='nc')
plt.plot(solution.t, solution.y[1]/ntitin, label='no')
plt.plot(solution.t, solution.y[2]/ntitin, label='np')
plt.plot(solution.t, sum_values, label='ns+np+no', linestyle='--')  # Dashed line for ns
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('Fraction of open/phosphorylated Titin')
plt.title('System Dynamics Over Time with Constraint')
plt.show()`

结果是这样的：

My Figure

力输入是正确的 - 我已经用附加图验证了这一点。

编辑：这是力输入功能：

def force_input(t):
    """
    Generates a force input as a function of time with 10 on-off repetitions followed by a rest.
    
    :param t: Time in seconds
    :return: Force in N at the given time
    """
    high_force = 20 * 10**(-12)  # N for 'on' state
    low_force = 4.5 * 10**(-12) # N for 'off' state and rest
    
    # Period for one 'on-off' cycle (20 seconds), and full period for the entire set (260 seconds)
    on_off_period = 20
    full_set_period = 260
    
    # Determine the position in the current set
    position_in_set = t % full_set_period
    
    if t<0:
        return low_force
    # Check if we are in the rest period
    if position_in_set > 10 * on_off_period:
        return low_force
    else:
        # If not in rest, determine if we are in 'on' or 'off' within the on-off period
        position_in_on_off = position_in_set % on_off_period
        return high_force if position_in_on_off < 10 else low_force

此外，如果我绘制前 10 个时间点的浓度，结果如下：

nc:  [0.99689787 0.99689787 0.99689785 0.99689779 0.99689774 0.99689769
 0.99689764 0.99689758 0.99689753 0.99689748]
no:  [ 0.0006      0.19058381 -0.33508352 -0.02302532 -0.01392478 -0.14079521
  0.01207641 -0.07641563  0.0918557   0.04086486]
np:  [ 0.0025     -0.18748378  0.3381835   0.0261254   0.01702492  0.14389538
 -0.00897615  0.07951593 -0.08875532 -0.03776443]
ns:  [1.00000000e-11 9.83700948e-12 1.02291693e-11 9.97153525e-12
 9.94892783e-12 1.00318165e-11 9.89773534e-12 9.95090262e-12
 9.80492439e-12 9.82905055e-12]

看起来 no 和 np 只是通过振荡甚至进入负空间来相互补偿。

python matplotlib ode differential-equations
1个回答
0
投票

我不确定我是否可以将其称为答案 - 只是一项正在进行的工作。然而。

你的方程非常僵硬。 （它们包含类似指数衰减的项。）基本上，您需要一个可以处理刚性方程的求解器；我建议“Radau”或“BDF”。通常我会寻找 

k.dt << 1
，其中 
k
是最大的速率常数。要强制时间步 
dt
，请使用 
max_step=
选项。 （它需要相当低）。

您应该使用满足“稳态”解（时间导数为零且无强迫）的值来初始化 

nc
、
n0
、
np
和 
ns
。您当前的价值观没有。

我不确定你的速率常数是否正确。该论文（或其补充材料）表明 

kp.ATP
大致为 
5.kr
，所以我冒昧地强制执行这一点。

使用您使用的强制，速率常数 

kplus
和 
kminus
与其他速率常数相比太小，无法做很多事情（尽管我不完全确定 
kr
的给定值在以下方面是正确的）纸张；看起来太大了。）我现在随意增加了
w0
，但还是有点不满意。

我不认为对 

nc+np+no+ns
施加约束在这里特别有帮助。它导致方程之间出现不自然的不对称性，方程本身应该很好地满足约束（只需将它们相加即可）。

最终，我认为最好的选择是给撰写论文的剑桥学者发送电子邮件，并准确询问他们使用哪些参数来创建图形：它对使用的数字非常敏感。

import numpy as nump
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Constants
kp = 0.07           # M^-1 s^-1
kr = 6              # s^-1
ks = 1.0e-7         # s^-1
kdns = 0.002        # s^-1
ATP = 7.2e-3        # M
w0 = 1.0e6          # s^-1
kb = 1.3806452e-23  # J/K
T = 310             # K
u0 = 1.0e-9         # m
umax = 28.0e-9      # m
delta_G = 35 * kb * T
ntitin = 4.7e-6     # in molar (M)

# OVER-RIDDEN VALUES - please see paper and supplementary material
kp = 5 * kr / ATP             # Paper says that ATP.kp is about 5 kr
w0 = 100 * w0                 # Otherwise forcing is just too small


def force_input(t):
    """
    Generates a force input as a function of time with 10 on-off repetitions followed by a rest.
    :param t: Time in seconds
    :return: Force in N at the given time
    """
    high_force = 20e-12     # N for 'on' state
    low_force = 4.5e-12     # N for 'off' state and rest
    
    # Period for one 'on-off' cycle (20 seconds), and full period for the entire set (260 seconds)
    on_off_period = 20
    full_set_period = 260
    
    # Determine the position in the current set
    position_in_set = t % full_set_period
    
    if t<0:
        return low_force
    # Check if we are in the rest period
    if position_in_set > 10 * on_off_period:
        return low_force
    else:
        # If not in rest, determine if we are in 'on' or 'off' within the on-off period
        position_in_on_off = position_in_set % on_off_period
        return high_force if position_in_on_off < 10 else low_force

#   return high_force if int( t / on_off_period ) % 2 == 1 else low_force      # Simpler alternative



# Differential equations
def equations(t, y):
    nc, no, np, ns = y
    
    # Calculate Eplus and Eminus for the current force value
    f = force_input(t)  # Call the previously defined force input function
    Eplus = delta_G - f * u0
    Eminus = f * (umax - u0)
    
    # Calculate rate constants for the current force value
    kplus  = w0 * nump.exp( -Eplus  / (kb * T))
    kminus = w0 * nump.exp( -Eminus / (kb * T))
#   print( kminus, kplus, kp * ATP, kr, kdns )                      # Check the size of rate constants
    
    # System of differential equations (without dns_dt since ns given through constraint )
    dnc_dt = -kplus * nc + kminus * no
    dno_dt =  kplus * nc - kminus * no - kp * ATP * no + kr * np + kdns * ns
    dnp_dt =                             kp * ATP * no - kr * np              - ks * np
    dns_dt =                                                     - kdns * ns  + ks * np

    return [dnc_dt, dno_dt, dnp_dt, dns_dt]



# Initial conditions - SATISFY STEADY STATE WITH NO FORCING
no0 = 0.0006 * ntitin                      # Fix this value
np0 = no0 * kp * ATP / ( ks + kr )         # Then everything else follows for the steady-state solution
ns0 = np0 * ks / kdns
nc0 = ntitin - no0 - np0 - ns0             # To all intents and purposes nc0 = ntitin
y0 = [ nc0, no0, np0, ns0 ]

# Time span for the solver
t_span = [-50, 3 * 260 ]                   # steady, followed by 3 full cycles of the input
t_values = nump.linspace(t_span[0], t_span[1], 100000 )

# Solve the system
solution = solve_ivp( equations, t_span, y0, t_eval=t_values, method='Radau', max_step=0.02 )

sum_values = ( solution.y[1] + solution.y[2] + solution.y[3] ) / ntitin

# Plotting the results
plt.plot(solution.t, solution.y[1]/ntitin, label='no')
plt.plot(solution.t, solution.y[2]/ntitin, label='np')
plt.plot(solution.t, sum_values, label='ns+np+no', linestyle='--')  # Dashed line for ns
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('Fraction of open/phosphorylated Titin')
plt.title('System Dynamics Over Time')
plt.show()

输出（注意：相当慢）：

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.