matlab-simulate-simbiology-model

SimBiology モデルをシミュレーションする

SimBiology モデルのシミュレーションを実行します。決定論的 ODE、確率論的 SSA、シナリオ探索、感度解析に対応しています。

使用する状況

• モデルの挙動を「シミュレーション」「実行」「予測」する場合
• 「もし〜なら」「〜したらどうなるか」（シミュレーションやシナリオを示します）
• モデルの時間経過に伴う結果を確認する場合
• 用量反応試験、パラメーター掃引、要因計画
• 確率論的、SSA、Gillespie、ノイズ、遺伝子発現のばらつき
• 「どのパラメーターが最も重要か」、感度、Sobol、Morris
• キーワード：「シミュレーション」「実行」「予測」「what-if」「確率論的」「感度」

使用しない状況

• モデルの構築やダイアグラムのレイアウト（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. + でシナリオを結合しない

+ 演算子は 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) を使用します。ある要因のすべてのエントリーが連続して並ぶとは想定しないでください。

判断表

シナリオ	方法
1 回限りのシミュレーション	sbiosimulate
パラメーター掃引／モンテカルロ	createSimFunction
用量／変異体／パラメーターの what-if	SimBiology.Scenarios + createSimFunction
分子数が少ない／ノイズがある	SSA ソルバー（cs.SolverType = 'ssa'）
どのパラメーターが重要かを調べる	sbiosobol（Sobol）または sbioelementaryeffects（Morris）
簡易的な感度確認	configset による局所感度

基本的なシミュレーション（sbiosimulate）

SimData を返す（出力を 1 つにする）方法を推奨します。状態名、単位、メタデータが含まれ、sbioplot および selectbyname で直接使用できます。

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

用量を指定する場合（用量を渡すときは、2 番目の引数として configset が必須です）。

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

結果のプロット

SimData を手早く可視化するには sbioplot を使用します。

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 を使わずにすべての状態へ数値として手早くアクセスするには、3 出力形式を使用します。

[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() の常に第 2 引数です。別の位置に置くとエラーになります。

用量スイープ

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

結果の順序（重要）

シナリオの結果は、ブロックごとではなく第 1 次元でインターリーブされます。 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 からアクセスします（構造体配列で、パラメーターごとに 1 つあります）

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 ではなく、モデルを直接渡します
• 境界行列：パラメーターごとに 1 行、列は [low high] です

参考資料

詳細なガイダンスは必要に応じて読み込んでください：

• references/stochastic-simulation-guidance.md — アンサンブルのプロット、分布の分析
• references/sensitivity-analysis-guidance.md — Sobol／Morris／Local の詳細なパターンと解釈

---

Copyright 2026 The MathWorks, Inc.

---