Simbiology Simulate Model

Simulate SimBiology Models

Run simulations of SimBiology models: deterministic ODE, stochastic SSA, scenario exploration, and sensitivity analysis.

When to Use

• "simulate", "run", "predict" model behavior
• "what if" / "what happens if" (implies simulation or scenarios)
• Time-course results from a model
• Dose-response studies, parameter sweeps, factorial designs
• Stochastic, SSA, Gillespie, noise, gene expression variability
• "which parameters matter most", sensitivity, Sobol, Morris
• Keywords: "simulate", "run", "predict", "what-if", "stochastic", "sensitivity"

When NOT to Use

• Model construction or diagram layout (use simbiology-build-model)
• Parameter estimation from data (use simbiology-fit-model)
• NCA / AUC / Cmax from data (use simbiology-fit-model)

Must-Follow Rules

0. Add helper scripts to the MATLAB path first

Run at the start of every session:

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

1. Element-wise operators in observables

Use ./ and .* (element-wise) in observable expressions when mixing time-varying species with constant parameters. Plain / and * cause size mismatches at simulation time.

2. StatesToLog for constant parameters

When observables reference constant parameters (e.g., Drug ./ Vd), add those parameters explicitly to StatesToLog:

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

StatesToLog = 'all' does not log constant compartments or parameters.

3. All reactions must be MassAction for SSA

The stochastic solver does not support custom rate expressions. Every reaction must use addkineticlaw(rx, 'MassAction').

4. Do NOT combine Scenarios with +

The + operator is not supported on SimBiology.Scenarios objects. Always use add() to append entries.

5. Reset local sensitivity options after use

Local sensitivity settings persist on the configset and affect subsequent simulations. Always reset:

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

6. Set MaximumWallClock to prevent hung simulations

When fitting or scanning, bad parameter values can make individual simulations extremely slow. Protect against this:

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

This is a configset property (not a solver or optimizer option). It stops any single simulation that exceeds the wall clock limit.

7. Unit conversion requires TimeUnits

When cs.CompileOptions.UnitConversion = true, you MUST also set cs.TimeUnits to match your StopTime units (e.g., 'hour'). Otherwise SimBiology defaults to seconds and your 24-unit simulation covers 24 seconds, not 24 hours:

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

8. Scenario results are interleaved, not blocked

Factorial scenario results come back interleaved by the first dimension. Always use generate(sc) to map result indices to conditions — never assume all entries of one factor appear consecutively.

Decision Table

Basic Simulation (sbiosimulate)

Prefer returning SimData (single output) — it carries state names, units, and metadata, and works directly with sbioplot and selectbyname:

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

With a dose (configset is required as 2nd argument when passing doses):

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

Plotting Results

Use sbioplot for quick visualization of SimData:

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

For custom plots, extract numeric data first:

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

Extracting State Data from SimData

Use selectbyname to extract specific states. It returns a SimData object, not a numeric array — extract numeric data before doing math:

matlabCopy code
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

Or use getdata() to get arrays:

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

For quick numeric access to all states without SimData, use the three-output form:

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

Acceleration (sbioaccelerate)

For repeated sbiosimulate calls on the same model, accelerate once first:

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

Rules:

• Call only right before sbiosimulate — not before fitting or analysis functions
• Valid after changing parameter/species values (e.g., p.Value = 0.2)
• Invalidated by structural changes (adding reactions, species, compartments) — must re-accelerate
• Do NOT use with createSimFunction, Scenarios, sensitivity, or fitting — these handle acceleration internally via AutoAccelerate

Repeated Simulation (createSimFunction)

matlabCopy code
% 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();

• Compiles once, runs many — much faster than sbiosimulate in a loop
• Exception: SSA (stochastic) requires sbiosimulate in a loop because each run needs fresh random state; createSimFunction does not support stochastic solvers
• Compatible with parfor (Parallel Computing Toolbox)
• Returns SimData objects; use .getdata() to extract arrays

SimFunction with doses

The 4th argument to createSimFunction declares which species receive doses. When executing, pass doses as a table (NOT a dose object):

matlabCopy code
% 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);

Common mistake: passing a sbiodose object to a SimFunction — this errors. Always convert to a table with Time and Amount columns.

Scenarios (SimBiology.Scenarios)

Systematically explore combinations of doses, variants, and parameters.

add() signature (argument order is critical)

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

The combination type ('cartesian' or 'elementwise') is always the second argument to add(). Putting it elsewhere errors.

Dose sweep

matlabCopy code
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]);

Full factorial (dose x parameter)

matlabCopy code
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);

Parameter sweep with dosed SimFunction

matlabCopy code
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

Result ordering (critical)

Scenario results are interleaved by the first dimension, not blocked. For a 2-dose × 3-ke factorial, results come back as:

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

Use generate(sc) to get a table mapping each result index to its conditions:

matlabCopy code
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

Never assume blocked ordering (all of Dose1 first, then all of Dose2). Always use generate(sc) to map results to conditions.

Entry types

Virtual population via distribution sampling

Scenarios can sample from probability distributions — use this for virtual patient simulations instead of manually generating parameter matrices:

matlabCopy code
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);

Steady-state with repeat dosing

matlabCopy code
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);

Stochastic Simulation (SSA)

For low molecule count systems where continuous ODE breaks down.

Single trajectory

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

Ensemble (multiple trajectories)

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

Gene expression template (all MassAction)

matlabCopy code
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'};

After SSA, reset solver: cs.SolverType = 'ode15s';

Sensitivity Analysis

Sobol (global, quantitative)

matlabCopy code
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

• First-order (Si): variance due to parameter alone
• Total-order (STi): variance due to parameter + all interactions
• Large gap STi - Si → strong interaction effects
• Access via sobolResults.SobolIndices(i).FirstOrder / .TotalOrder (struct array, one per parameter)

Morris screening (global, ranking)

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

• High mean effect → influential parameter
• High standard deviation → nonlinear or interaction effects

Local sensitivity

matlabCopy code
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 = [];

Conventions

• Loading .sbproj files: sbioloadproject returns a struct with the model name as field — extract dynamically:

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

• Pass model objects (not UUID strings) to simulation functions
• Use getModelByUUID(uuid) to recover handles (provided by this skill's scripts/ directory — add to path at session start)
• createSimFunction returns SimData; extract with .getdata()
• Pass model directly to sbiosobol/sbioelementaryeffects (not a SimFunction)
• Bounds matrix: one row per parameter, columns [low high]

References

Load on demand for detailed guidance:

• references/stochastic-simulation-guidance.md — ensemble plotting, distribution analysis
• references/sensitivity-analysis-guidance.md — full Sobol/Morris/Local patterns and interpretation

Copyright 2026 The MathWorks, Inc.