matlab-simulate-simbiology-model

模拟 SimBiology 模型

运行 SimBiology 模型模拟：确定性 ODE、随机 SSA、场景探索和灵敏度分析。

适用场景

• “模拟”“运行”“预测”模型行为
• “如果……会怎样”／“假如”（表示需要进行模拟或场景分析）
• 获取模型的时间过程结果
• 剂量-响应研究、参数扫描、析因设计
• 随机、SSA、Gillespie、噪声、基因表达变异
• “哪些参数最重要”、灵敏度、Sobol、Morris
• 关键词：“模拟”“运行”“预测”“假设分析”“随机”“灵敏度”

不适用场景

• 构建模型或调整图表布局（使用 matlab-build-simbiology-model）
• 根据数据估计参数（使用 matlab-fit-simbiology-model）
• 根据数据进行 NCA／AUC／Cmax 分析（使用 matlab-fit-simbiology-model）

必须遵循的规则

0. 首先将辅助脚本添加到 MATLAB 路径

每次会话开始时运行：

addpath(fullfile('<WORKSPACE_ROOT>', '.claude', 'skills', 'matlab-simulate-simbiology-model', 'scripts'));

1. 可观测量中的逐元素运算符

在可观测量表达式中混合使用随时间变化的物种和常量参数时，应使用 ./ 和 .*（逐元素运算）。普通的 / 和 * 会在模拟时导致尺寸不匹配。

2. 常量参数的 StatesToLog

当可观测量引用常量参数（例如 Drug ./ Vd）时，需将这些参数显式添加到 StatesToLog：

cs.RuntimeOptions.StatesToLog = [m.Species; sbioselect(m,'Type','parameter','Name','Vd')];

StatesToLog = 'all' 不会记录常量室或参数。

3. SSA 要求所有反应均为 MassAction

随机求解器不支持自定义速率表达式。每个反应都必须使用 addkineticlaw(rx, 'MassAction')。

4. 不要将 Scenarios 与 + 组合

+ 运算符不支持用于 SimBiology.Scenarios 对象。 始终使用 add() 追加条目。

5. 使用后重置局部灵敏度选项

局部灵敏度设置会保留在 configset 上，并影响后续模拟。务必重置：

cs.SolverOptions.SensitivityAnalysis = false;
cs.SensitivityAnalysisOptions.Inputs = [];
cs.SensitivityAnalysisOptions.Outputs = [];

6. 设置 MaximumWallClock，避免模拟卡住

拟合或扫描时，不合适的参数值可能导致单次模拟极其缓慢。请采取措施防止这种情况：

cs.MaximumWallClock = 60;  % seconds; default is Inf

这是 configset 属性（不是求解器或优化器选项）。 如果单次模拟超过墙钟时间限制，该属性会将其停止。

7. 单位转换需要设置 TimeUnits

当 cs.CompileOptions.UnitConversion = true 时，还必须将 cs.TimeUnits 设置为与 StopTime 单位一致（例如，'hour'）。 否则 SimBiology 会默认使用秒，24 个单位的模拟时长就会是 24 秒，而不是 24 小时：

cs.CompileOptions.UnitConversion = true;
cs.TimeUnits = 'hour';
cs.StopTime = 24;  % now correctly 24 hours

8. 场景结果是交错排列，而非分块排列

析因场景的结果会按照第一个维度交错排列。始终使用 generate(sc) 将结果索引映射到条件——切勿假设某个因子的所有条目都连续排列。

决策表

场景	方法
一次性模拟	sbiosimulate
参数扫描／蒙特卡洛	createSimFunction
剂量／变体／参数假设分析	SimBiology.Scenarios + createSimFunction
分子数量少／噪声	SSA 求解器（cs.SolverType = 'ssa'）
哪些参数最重要？	sbiosobol（Sobol）或 sbioelementaryeffects（Morris）
快速灵敏度检查	通过 configset 进行局部灵敏度分析

基本模拟（sbiosimulate）

优先返回 SimData（单个输出）——它包含状态名称、单位和元数据，并且可直接与 sbioplot 和 selectbyname 配合使用：

m = getModelByUUID(modelId);
cs = getconfigset(m, 'active');
cs.StopTime = 24;
cs.SolverType = 'ode15s';
simData = sbiosimulate(m);

使用剂量时（传入剂量时，configset 必须作为第 2 个参数）：

d = sbiodose('Bolus', 'schedule');
d.TargetName = 'Drug'; d.Amount = 100; d.Time = 0;
simData = sbiosimulate(m, cs, d);  % NOT sbiosimulate(m, d) — errors

绘制结果

使用 sbioplot 快速可视化 SimData：

simData = sbiosimulate(m, cs, d);
sbioplot(simData);

要自定义绘图，请先提取数值数据：

[t, x, names] = getdata(simData);
plot(t, x);
legend(names, 'Interpreter', 'none');
xlabel('Time'); ylabel('Amount');

从 SimData 中提取状态数据

使用 selectbyname 提取特定状态。它返回的是 SimData 对象，而非数值数组——进行数学运算前请先提取数值数据：

simData = sbiosimulate(m, cs, d);
result = selectbyname(simData, 'Central.Drug');  % returns SimData, NOT double
drugData = result.Data;   % numeric column vector
drugTime = result.Time;   % time column vector

或者使用 getdata() 获取数组：

[t, x, names] = getdata(selectbyname(simData, 'Central.Drug'));

如果要在不使用 SimData 的情况下快速访问所有状态的数值，请使用三输出形式：

[t, x, names] = sbiosimulate(m, cs, d);  % t, x are double arrays directly

加速（sbioaccelerate）

对于同一模型的重复 sbiosimulate 调用，请先加速一次：

sbioaccelerate(m);
simData = sbiosimulate(m);  % faster

规则：

• 仅在调用 sbiosimulate 前立即调用——不要在拟合或分析函数前调用
• 更改参数／物种值后仍然有效（例如，p.Value = 0.2）
• 结构发生变化（添加反应、物种、区室）后失效——必须重新加速
• 不要与 createSimFunction、Scenarios、敏感性分析或拟合搭配使用——这些功能通过 AutoAccelerate 在内部处理加速

重复仿真（createSimFunction）

% Signature: createSimFunction(model, params, observables, dosedSpecies)
simfun = createSimFunction(model, {'ke','ka'}, {'Drug'}, []);
r1 = simfun([0.1, 0.5], 24);           % single run
r2 = simfun([0.1, 0.5; 0.3, 1.0], 24); % multiple parameter sets (rows)
[t, x] = r1.getdata();

• 编译一次，运行多次——比在循环中调用 sbiosimulate 快得多
• **例外：**SSA（随机）需要在循环中调用 sbiosimulate，因为每次运行都需要全新的随机状态；createSimFunction 不支持随机求解器
• 兼容 parfor（Parallel Computing Toolbox）
• 返回 SimData 对象；使用 .getdata() 提取数组

带剂量的 SimFunction

createSimFunction 的第 4 个参数用于声明哪些物种接收剂量。执行时，请将剂量作为表传入（而不是剂量对象）：

% Create: specify dosed species names in 4th argument
simfun = createSimFunction(model, {'ke'}, {'Drug'}, {'Drug'});

% Execute: pass dose as a table with Time and Amount columns
doseTable = table(0, 100, 'VariableNames', {'Time', 'Amount'});
result = simfun(0.1, 24, doseTable);

% Multiple dose events
multiDose = table([0; 12], [100; 50], 'VariableNames', {'Time', 'Amount'});
result = simfun(0.1, 24, multiDose);

% Multiple dosed species: cell array of tables (one per species, same order)
simfun2 = createSimFunction(model, {'ke'}, {'Drug','Drug2'}, {'Drug','Drug2'});
doses = {doseTable1, doseTable2};
result = simfun2(0.1, 24, doses);

**常见错误：**将 sbiodose 对象传给 SimFunction——这会导致报错。务必将其转换为包含 Time 和 Amount 列的表。

场景（SimBiology.Scenarios）

系统地探索剂量、变体和参数的组合。

add() 签名（参数顺序至关重要）

add(sc, combination, name, values, ...)
%     ^^^^^^^^^^^^^
%     MUST be 2nd argument: 'cartesian' or 'elementwise'

组合类型（'cartesian' 或 'elementwise'）始终是 add() 的第二个参数。放在其他位置会导致报错。

剂量扫描

d1 = sbiodose('Low','schedule'); d1.TargetName = 'Drug'; d1.Amount = 50; d1.Time = 0;
d2 = sbiodose('High','schedule'); d2.TargetName = 'Drug'; d2.Amount = 200; d2.Time = 0;
sc = SimBiology.Scenarios('DoseLevel', [d1, d2]);

全因子（剂量 x 参数）

sc = SimBiology.Scenarios('DoseLevel', [d1, d2]);
add(sc, 'cartesian', 'ke', [0.05 0.1 0.2]);  % 2 x 3 = 6 combinations
simfun = createSimFunction(model, sc, {'Drug'}, []);
results = simfun(sc, 24);

带剂量 SimFunction 的参数扫描

sc = SimBiology.Scenarios('ke', [0.05 0.1 0.2]);
simfun = createSimFunction(model, sc, {'Drug'}, {'Drug'});
doseTable = table(0, 100, 'VariableNames', {'Time', 'Amount'});
results = simfun(sc, 24, doseTable);  % dose table as 3rd argument

结果顺序（至关重要）

场景结果按第一个维度交错排列，而非分组排列。 对于 2 剂量 × 3 ke 的全因子，结果返回如下：

results(1): Dose1, ke1
results(2): Dose2, ke1
results(3): Dose1, ke2
results(4): Dose2, ke2
results(5): Dose1, ke3
results(6): Dose2, ke3

使用 generate(sc) 获取一张表，将每个结果索引映射到其条件：

genTable = generate(sc);  % table with one row per scenario
for i = 1:numel(results)
    [t, x] = results(i).getdata();
    fprintf('Dose=%s, ke=%.2f: Drug at t=end = %.2f\n', ...
        genTable.DoseLevel(i).Name, genTable.ke(i), x(end,1));
end

切勿假设结果按组排列（先是所有 Dose1，再是所有 Dose2）。 务必使用 generate(sc) 将结果映射到条件。

条目类型

内容类型	示例
剂量向量	SimBiology.Scenarios('DoseLevel', [d1, d2])
变体向量	SimBiology.Scenarios('Pop', [v1, v2])
参数值	SimBiology.Scenarios('ke', [0.05 0.1 0.2])
物种值	SimBiology.Scenarios('Drug', [50 100 200])
概率分布	add(sc, 'elementwise', 'ke', makedist('Lognormal',...), 'Number', 50)

通过分布采样生成虚拟种群

Scenarios 可以从概率分布中采样——可用此方法进行虚拟患者仿真，而不是手动生成参数矩阵：

pd = makedist('Lognormal', 'mu', log(0.1), 'sigma', 0.3);
sc = SimBiology.Scenarios;
add(sc, 'elementwise', 'ke', pd, 'Number', 50);
simfun = createSimFunction(model, sc, {'Drug'}, []);
results = simfun(sc, 24);

重复给药下的稳态

d = sbiodose('RepeatDose', 'repeat');
d.TargetName = 'Drug'; d.Amount = 100;
d.StartTime = 0; d.Interval = 12; d.RepeatCount = 50;
cs.StopTime = d.Interval * (d.RepeatCount + 1);
[t, x, names] = sbiosimulate(model, cs, d);

随机仿真（SSA）

适用于连续 ODE 不再适用的低分子数量系统。

单条轨迹

cs = getconfigset(model, 'active');
cs.SolverType = 'ssa';
cs.StopTime = 100;
simData = sbiosimulate(model);
[t, x, names] = getdata(simData);

集合（多条轨迹）

nRuns = 200;
allResults = cell(nRuns, 1);
for i = 1:nRuns
    allResults{i} = sbiosimulate(model);
end

基因表达模板（全部为 MassAction）

model = sbiomodel('GeneExpr');
comp = addcompartment(model, 'cell');
addspecies(comp, 'Gene', 1);
addspecies(comp, 'mRNA', 0);
addspecies(comp, 'Protein', 0);
addparameter(model, 'k_txn', 0.1);
addparameter(model, 'k_tln', 0.5);
addparameter(model, 'k_mdeg', 0.05);
addparameter(model, 'k_pdeg', 0.01);

% Transcription: Gene -> Gene + mRNA (Gene is catalyst)
rx1 = addreaction(model, 'Gene -> Gene + mRNA');
kl1 = addkineticlaw(rx1, 'MassAction'); kl1.ParameterVariableNames = {'k_txn'};
% Translation: mRNA -> mRNA + Protein
rx2 = addreaction(model, 'mRNA -> mRNA + Protein');
kl2 = addkineticlaw(rx2, 'MassAction'); kl2.ParameterVariableNames = {'k_tln'};
% Degradation
rx3 = addreaction(model, 'mRNA -> null');
kl3 = addkineticlaw(rx3, 'MassAction'); kl3.ParameterVariableNames = {'k_mdeg'};
rx4 = addreaction(model, 'Protein -> null');
kl4 = addkineticlaw(rx4, 'MassAction'); kl4.ParameterVariableNames = {'k_pdeg'};

SSA 之后，重置求解器：cs.SolverType = 'ode15s';

敏感性分析

Sobol（全局、定量）

bounds = [0.01 1; 0.1 5];  % [low high] per parameter
sobolResults = sbiosobol(m, {'ke','ka'}, {'Drug'}, ...
    'OutputTimes', 0:1:24, 'NumberSamples', 500, 'Bounds', bounds);
plot(sobolResults);

% Extract indices from struct array
for i = 1:numel(sobolResults.SobolIndices)
    Si  = mean(sobolResults.SobolIndices(i).FirstOrder, 'omitnan');
    STi = mean(sobolResults.SobolIndices(i).TotalOrder, 'omitnan');
    fprintf('%s: Si=%.3f, STi=%.3f\n', sobolResults.SobolIndices(i).Parameter, Si, STi);
end

• **一阶（Si）：**由参数单独引起的方差
• **总阶（STi）：**由参数及所有交互作用引起的方差
• STi - Si 差值较大 → 交互作用效应强
• 通过 sobolResults.SobolIndices(i).FirstOrder / .TotalOrder 访问（结构体数组，每个参数对应一个元素）

Morris 筛选（全局、排序）

bounds = [0.01 1; 0.1 5];
eeResults = sbioelementaryeffects(m, {'ke','ka'}, {'Drug'}, ...
    'OutputTimes', 0:1:24, 'NumberSamples', 50, 'Bounds', bounds);

• 平均效应高 → 参数影响较大
• 标准差高 → 存在非线性或交互作用效应

局部敏感性

cs.SolverOptions.SensitivityAnalysis = true;
cs.SensitivityAnalysisOptions.Normalization = 'Full';
cs.SensitivityAnalysisOptions.Inputs = sbioselect(m,'Type','parameter','Name',{'ke','ka'});
cs.SensitivityAnalysisOptions.Outputs = sbioselect(m,'Type','species','Name','Drug');
simData = sbiosimulate(m);
[t, R] = getsensmatrix(simData);
% IMPORTANT: Reset after use
cs.SolverOptions.SensitivityAnalysis = false;
cs.SensitivityAnalysisOptions.Inputs = [];
cs.SensitivityAnalysisOptions.Outputs = [];

归一化	含义
'None'	原始 dY/dp
'Half'	(p/y) dY/dp
'Full'	无量纲；两侧均已归一化

约定

• 加载 .sbproj 文件： sbioloadproject 返回一个结构体，其中模型名称是一个字段 — 动态提取：

  proj = sbioloadproject('file.sbproj');
  fn = fieldnames(proj);
  model = proj.(fn{1});
  

• 将模型对象（而非 UUID 字符串）传递给仿真函数
• 使用 getModelByUUID(uuid) 获取句柄（由此技能的 scripts/ 目录提供 — 在会话开始时将其添加到路径中）
• createSimFunction 返回 SimData；使用 .getdata() 提取数据
• 直接将模型传递给 sbiosobol/sbioelementaryeffects（而非 SimFunction）
• 边界矩阵：每个参数占一行，列为 [low high]

参考资料

按需加载，获取详细指导：

• references/stochastic-simulation-guidance.md — 集合绘图、分布分析
• references/sensitivity-analysis-guidance.md — 完整的 Sobol/Morris/Local 模式及解读

---

版权所有 2026 The MathWorks, Inc.

---