661 微生物群落模拟建模¶

dXi/dt = Xi × (ri + Σj aij × Xj)

Xi = 物种i的丰度
ri = 物种i的内禀增长率
aij = 物种j对物种i的互作系数（正=促进，负=抑制）

# gLV模型模拟微生物群落动态
import numpy as np  # 数值计算
from scipy.integrate import odeint  # ODE求解器
import matplotlib.pyplot as plt  # 绑图

def glv_model(X, t, r, A):
    """
    广义Lotka-Volterra模型
    X: 物种丰度向量
    t: 时间
    r: 内禀增长率向量
    A: 互作矩阵（n×n）
    """
    dXdt = X * (r + A @ X)  # gLV方程
    dXdt[X < 1e-10] = 0  # 丰度接近0时停止衰减
    return dXdt

# 3个物种的示例
n_species = 3  # 物种数
r = np.array([0.5, 0.3, 0.4])  # 内禀增长率

# 互作矩阵：对角线为自抑制，非对角为种间互作
A = np.array([
    [-0.1, -0.05, 0.02],   # 物种1：自抑制-0.1，被物种2抑制，被物种3促进
    [0.03, -0.08, -0.04],  # 物种2：被物种1促进，自抑制，被物种3抑制
    [-0.02, 0.06, -0.12]   # 物种3：被物种1抑制，被物种2促进，自抑制
])

# 初始条件和时间
X0 = np.array([1.0, 0.5, 0.8])  # 初始丰度
t = np.linspace(0, 100, 1000)  # 时间范围

# 求解ODE
solution = odeint(glv_model, X0, t, args=(r, A))  # 数值求解

# 可视化
fig, ax = plt.subplots(figsize=(10, 5))
species_names = ["Bacteroides", "Faecalibacterium", "Bifidobacterium"]
for i, name in enumerate(species_names):
    ax.plot(t, solution[:, i], label=name, linewidth=2)  # 画每个物种的丰度曲线
ax.set_xlabel("时间")  # x轴标签
ax.set_ylabel("丰度")  # y轴标签
ax.legend()  # 图例
ax.set_title("gLV模型群落动态模拟")  # 标题
plt.tight_layout()
plt.savefig("glv_dynamics.png", dpi=150)  # 保存图片

# 用MIMIC包从时间序列数据推断gLV参数
# pip install mimic-microbiome（2025年新包）
from scipy.optimize import minimize  # 优化器

def infer_glv_params(time_series_data, time_points):
    """
    从纵向数据推断gLV参数
    time_series_data: [n_timepoints, n_species] 丰度矩阵
    time_points: 时间点数组
    """
    n_species = time_series_data.shape[1]  # 物种数
    n_params = n_species + n_species**2  # 参数数：r(n) + A(n×n)

    def objective(params):
        """目标函数：最小化预测与观测的差异"""
        r = params[:n_species]  # 提取增长率
        A = params[n_species:].reshape(n_species, n_species)  # 提取互作矩阵

        # 用gLV预测
        predicted = odeint(glv_model, time_series_data[0], time_points, args=(r, A))

        # 均方误差
        mse = np.mean((predicted - time_series_data)**2)
        # 加L1正则化促进稀疏互作
        reg = 0.01 * np.sum(np.abs(A))  # 稀疏性惩罚
        return mse + reg

    # 初始参数
    x0 = np.random.randn(n_params) * 0.01  # 随机初始化

    # 优化
    result = minimize(objective, x0, method='L-BFGS-B',
                      options={'maxiter': 1000})  # L-BFGS-B优化

    # 提取结果
    r_est = result.x[:n_species]  # 推断的增长率
    A_est = result.x[n_species:].reshape(n_species, n_species)  # 推断的互作矩阵

    return r_est, A_est

# COMETS群落空间模拟
import cometspy as c  # COMETS Python接口
import cobra  # COBRA工具

# 1. 加载代谢模型
ecoli = cobra.io.read_sbml_model("e_coli_core.xml")  # 大肠杆菌核心模型
ecoli_model = c.model(ecoli)  # 转为COMETS模型
ecoli_model.initial_pop = [4, 4, 1e-6]  # 初始位置(x,y)和生物量

# 2. 创建空间布局
layout = c.layout([ecoli_model])  # 创建含一个物种的布局
layout.grid = [10, 10]  # 10×10空间网格

# 3. 设置培养基
layout.set_specific_metabolite("glc__D_e", 0.001)  # 葡萄糖浓度
layout.set_specific_metabolite("o2_e", 1000)  # 氧气充足
layout.set_specific_metabolite("nh4_e", 1000)  # 氮源充足

# 4. 设置模拟参数
params = c.params()
params.set_param("maxCycles", 200)  # 最大模拟步数
params.set_param("timeStep", 0.1)  # 时间步长
params.set_param("spaceWidth", 0.02)  # 空间单元宽度（cm）
params.set_param("defaultDiffConst", 5e-6)  # 默认扩散系数

# 5. 运行模拟
sim = c.comets(layout, params)
sim.run()  # 运行

# 6. 分析结果
biomass = sim.total_biomass  # 总生物量随时间变化
media = sim.media  # 培养基成分随时间变化
print(f"最终总生物量: {biomass.iloc[-1]['e_coli_core']:.6f}")  # 最终生物量

# 7. 空间分布可视化
# COMETS输出的空间数据可以用matplotlib绘制热图

# BacArena个体水平代谢建模（R语言）
library(BacArena)  # 加载BacArena包
library(sybil)     # 加载代谢模型求解器

# 1. 加载代谢模型
mod1 <- readSBMLmod("e_coli.xml")    # 加载大肠杆菌模型
mod2 <- readSBMLmod("b_theta.xml")   # 加载拟杆菌模型

# 2. 创建细菌对象
bac1 <- Bac(mod1,
            deathrate = 0.01,         # 死亡率
            duprate = 0.5,            # 复制速率
            speed = 2)                # 移动速度
bac2 <- Bac(mod2,
            deathrate = 0.01,
            duprate = 0.3,
            speed = 1)

# 3. 创建模拟环境
arena <- Arena(n = 50, m = 50)  # 50×50网格
arena <- addOrg(arena, bac1, amount = 50)   # 添加50个大肠杆菌
arena <- addOrg(arena, bac2, amount = 50)   # 添加50个拟杆菌

# 4. 设置培养基
arena <- addSubs(arena, smax = 0.5,
                 mediac = "EX_glc__D_e",     # 添加葡萄糖
                 unit = "mM")
arena <- addSubs(arena, smax = 1.0,
                 mediac = "EX_o2_e",         # 添加氧气
                 unit = "mM")

# 5. 运行模拟
sim <- simEnv(arena, time = 20)  # 模拟20个时间步

# 6. 可视化结果
plotCurves(sim)    # 生长曲线
evalArena(sim, show_legend = TRUE)  # 空间分布

# MIMIC：2025年Bioinformatics发布的Python包
# pip install mimic-microbiome
# 集成gLV、高斯过程(GP)、向量自回归(VAR)等模型

# 概念使用示例
from mimic import gLV, GaussianProcess, VAR

# 选择模型
model = gLV()  # 使用gLV模型
model.fit(abundance_data, time_points)  # 拟合数据
predictions = model.predict(future_times)  # 预测未来

# 也可以用贝叶斯推断
model_bayes = gLV(inference="bayesian")
model_bayes.fit(abundance_data, time_points)

Virtual Colon特点：
- 基于BacArena框架扩展
- 模拟肠道空间结构（粘液层、肠壁）
- 空间依赖的扩散系数
- 宿主-微生物互作模拟
- 发表于mSystems，开源

# 模型选择决策
只关注种间互作强度 → gLV（最简单）
需要代谢机制解释 → COMETS（群体FBA+空间）
需要个体异质性 → BacArena（个体水平）
需要集成多种模型 → MIMIC (2025)
需要模拟肠道空间 → Virtual Colon (2025)
需要营养-微生物抑制 → BN-BacArena (2024)

# gLV参数含义
r > 0: 物种能自发增长
aij > 0: 物种j促进物种i（互利/偏利）
aij < 0: 物种j抑制物种i（竞争/捕食）
aii < 0: 自抑制（必须，保证稳定性）