基于Python的稳态调控建模：多物种动态系统仿真优化求助

核心修改思路 #

针对你的需求，我们从多物种系统重构、昼夜节律项整合、参数批量测试三个维度扩展现有代码，以下是具体实现步骤：

1. 多物种系统搭建 #

将单物种的标量参数改为数组/矩阵，适配n个物种的相互作用：

• 参数定义：
  • q：n维数组，每个元素对应一个物种的生成率q_sigma_i
  • chi0：n维数组，每个元素对应一个物种的初始清除率chi_sigma_i_0
  • beta：n×n矩阵，beta[i][j]表示物种j对物种i清除率的阻塞系数beta_sigma_j
• 模型函数重构：输入的H变为n维数组，逐个计算每个物种的清除率和稳态压力变化率

2. 昼夜节律调控引入 #

根据论文设定，昼夜节律项C(t) = sin(ωt)（ω=2π/24，对应24小时周期）采用状态切换方式适配睡眠/清醒的参数差异：

• 当C(t) > 0时使用清醒状态参数（生成率高、清除率低）
• 当C(t) ≤ 0时使用睡眠状态参数（生成率低、清除率高）
• 若需直接将节律项融入清除率公式，可修改chi_i的计算逻辑为：chi_i(t) = chi0_i * (1 + amp*sin(2π*t/24)) / (1 + np.sum(beta[i] * H))（amp为节律幅度）

3. 参数测试与仿真实现 #

支持不同初始值、beta矩阵、q数组的组合测试，以下是完整代码示例：

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

# --------------------------
# 多物种与昼夜节律参数配置
# --------------------------
n_species = 3  # 物种数量
omega = 2 * np.pi / 24  # 昼夜节律角频率（24小时周期）

# 清醒状态（w）参数：生成率高，清除率低
q_w = np.array([7e-6, 65e-6, 140e-6])
chi0_w = np.array([0.04, 0.05, 0.06])
beta_w = np.array([
    [0, 0.1, 0.05],   # 物种1的清除率受物种2、3阻塞
    [0.05, 0, 0.2],   # 物种2的清除率受物种1、3阻塞
    [0.1, 0.05, 0]    # 物种3的清除率受物种1、2阻塞
])

# 睡眠状态（s）参数：生成率低，清除率高
q_s = q_w * 0.2
chi0_s = chi0_w * 1.5
beta_s = beta_w * 0.8

# --------------------------
# 多物种微分方程模型
# --------------------------
def multi_species_model(H, t):
    # 计算当前昼夜节律项
    C = np.sin(omega * t)
    
    # 根据节律切换睡眠/清醒参数
    q = q_w if C > 0 else q_s
    chi0 = chi0_w if C > 0 else chi0_s
    beta = beta_w if C > 0 else beta_s
    
    # 逐个计算每个物种的稳态压力变化率
    dH_dt = np.zeros(n_species)
    for i in range(n_species):
        chi_i = chi0[i] / (1 + np.sum(beta[i] * H))
        dH_dt[i] = q[i] - chi_i * H[i]
    
    return dH_dt

# --------------------------
# 仿真参数与测试组合
# --------------------------
total_time = 24 * 8  # 8天（1周）
time_points = np.linspace(0, total_time, 2000)  # 更密集的时间点

# 测试不同初始条件和beta组合
test_cases = [
    {"name": "初始H=0, 正常阻塞", "H0": np.zeros(n_species)},
    {"name": "初始H=10, 正常阻塞", "H0": np.full(n_species, 10)},
    {"name": "初始H=0, 无阻塞", "H0": np.zeros(n_species), "beta_scale": 0}
]

# 运行所有测试用例
results = []
for case in test_cases:
    if "beta_scale" in case:
        # 临时修改beta矩阵实现无阻塞测试
        beta_w_temp = beta_w * case["beta_scale"]
        beta_s_temp = beta_s * case["beta_scale"]
        
        def temp_model(H, t):
            C = np.sin(omega * t)
            q = q_w if C > 0 else q_s
            chi0 = chi0_w if C > 0 else chi0_s
            beta = beta_w_temp if C > 0 else beta_s_temp
            
            dH_dt = np.zeros(n_species)
            for i in range(n_species):
                chi_i = chi0[i] / (1 + np.sum(beta[i] * H))
                dH_dt[i] = q[i] - chi_i * H[i]
            return dH_dt
        
        res = odeint(temp_model, case["H0"], time_points)
    else:
        res = odeint(multi_species_model, case["H0"], time_points)
    results.append({"case": case["name"], "data": res})

# --------------------------
# 结果可视化
# --------------------------
fig, axes = plt.subplots(len(test_cases), 1, figsize=(12, 8), sharex=True)
species_labels = [f"物种{i+1}" for i in range(n_species)]

for idx, result in enumerate(results):
    ax = axes[idx]
    for i in range(n_species):
        ax.plot(time_points, result["data"][:, i], label=species_labels[i])
    ax.set_title(result["case"])
    ax.set_ylabel("稳态压力 H_i(t)")
    ax.legend()

axes[-1].set_xlabel("时间 (小时)")
plt.tight_layout()
plt.show()

关键说明 #

1. 多物种交互：通过n×n的beta矩阵实现物种间的阻塞效应，可根据论文需求调整矩阵元素（如对角阵表示物种自阻塞，全零矩阵表示无阻塞）
2. 昼夜节律适配：代码采用状态切换逻辑，贴合论文中睡眠时清除率升高的核心设定，可自由修改参数比例调整昼夜差异
3. 参数测试：通过test_cases列表批量测试不同初始条件和beta组合，可直接扩展添加更多q值、节律幅度的测试场景

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