基于Scipy odeint实现ODE参数随状态变量阈值自动切换的需求

Hey there! Let's fix this parameter switching issue for your 4th-order ODE system. Instead of guessing when IIa will hit your threshold, we can dynamically check the IIa concentration at each integration step and flip the kIIaAP parameter automatically once the threshold is crossed. Here's how to adjust your code:

Key Changes Explained #

• We'll maintain a current_kIIaAP variable that starts at 0.0 (your initial requirement).
• After each integration step, we check if the current IIa value meets or exceeds IIath. If yes, we set current_kIIaAP to 0.5 and keep it there permanently.
• No more pre-defining the kIIaAP array—this makes the code adaptive to any initial conditions or parameter tweaks, no manual guesswork needed.

Modified Full Code #

from math import *
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

plt.close("all")

IIath = 4e-8  # Thrombin concentration threshold for plat. activation. (kg/kg)
PL = 0.0016   # Platelet concentration (kg/kg == 200e9pl/L)

def model(x, t, kIIaAP):
    IIa = x[0]
    II = x[1]
    AP = x[2]
    RP = x[3]
    
    kAPAP = 5.24e-2    # Platelet activation by activated platelets (s-1)
    kAPII = 0.73       # Thrombin generation by activated platelets (s-1)
    ksurf = 7.3e-6
    kin0 = 1.7e-2      # Inhibition of thrombin in interior cells
    kinadd = 3e-3      # Enhanced inhibition of thrombin due to ATH III in boundary cells

    dIIadt = -(kin0 + kinadd)*IIa + (ksurf + kAPII*AP)*II
    dIIdt = -(ksurf + kAPII*AP)*II
    dAPdt = kAPAP * AP * RP + kIIaAP*RP
    dRPdt = -kAPAP*AP*RP - kIIaAP*RP
    return [dIIadt, dIIdt, dAPdt, dRPdt]

# Initial conditions
y0 = [0, 9.51e-5, 0, 0.0016]
tend = 2000
nsteps = 4000
dt = tend / nsteps
t = np.linspace(0, tend, nsteps)

# Initialize arrays to store results
IIa = np.zeros(nsteps)
II = np.zeros(nsteps)
AP = np.zeros(nsteps)
RP = np.zeros(nsteps)
# Optional: Track kIIaAP over time to verify the switch
kIIaAP_history = np.zeros(nsteps)

# Start with kIIaAP = 0
current_kIIaAP = 0.0
kIIaAP_history[0] = current_kIIaAP
IIa[0] = y0[0]
II[0] = y0[1]
AP[0] = y0[2]
RP[0] = y0[3]

# Solve step-by-step with dynamic parameter switching
for i in range(1, nsteps):
    tspan = [t[i-1], t[i]]
    # Use the current kIIaAP value for this time step
    y = odeint(model, y0, tspan, args=(current_kIIaAP,))
    y0 = y[1]  # Update initial condition for next step
    
    # Store current state
    IIa[i] = y0[0]
    II[i] = y0[1]
    AP[i] = y0[2]
    RP[i] = y0[3]
    
    # Check if IIa has crossed the threshold: switch kIIaAP to 0.5 if yes
    if y0[0] >= IIath:
        current_kIIaAP = 0.5
    kIIaAP_history[i] = current_kIIaAP  # Log parameter changes

def kgkg2molL(value):
    molmassIIa = 70e3  # Molar mass IIa [Da]
    M = (value * 1000 / molmassIIa) / 1.06  # Blood density = 1060 kg/m³
    nM = M * 10**9
    return nM

# Convert IIa to nM for plotting
IIanM = [kgkg2molL(i) for i in IIa]

# Plot results (added kIIaAP history for verification)
fig, axs = plt.subplots(5, 1, figsize=(8, 12))

axs[0].plot(t, IIanM)
axs[0].axhline(y=kgkg2molL(IIath), color="r", linestyle="--", label=f"Threshold ({kgkg2molL(IIath):.2f} nM)")
axs[0].set_ylabel('[IIa] nM')
axs[0].legend()

axs[1].plot(t, II, 'b--')
axs[1].set_ylabel('[II] kg/kg')

axs[2].plot(t, AP)
axs[2].set_ylabel('[AP] kg/kg')

axs[3].plot(t, RP)
axs[3].set_ylabel('[RP] kg/kg')

axs[4].plot(t, kIIaAP_history, 'g-')
axs[4].set_ylabel('kIIaAP value')
axs[4].set_xlabel('time (s)')

plt.tight_layout()
plt.show()

What's Improved? #

1. Adaptive Parameter Switching: The current_kIIaAP variable updates automatically the moment IIa crosses IIath—no more manual guessing of when the switch should happen.
2. Verification Plot: I added a plot for kIIaAP_history so you can see exactly when the parameter flips, making it easy to confirm the logic works as expected.
3. Threshold Alignment: The red threshold line in the IIa plot is converted to nM to match the y-axis, so you can visually confirm the switch lines up with the threshold crossing.

This approach ensures your model responds in real-time to the state variable, making it robust to any changes in initial conditions or other parameters.

内容的提问来源于stack exchange，提问作者jacktat