matlab-simulate-simbiology-model

name	matlab-simulate-simbiology-model
description	Simulate SimBiology models — ODE, stochastic (SSA), scenarios, and sensitivity analysis. Use when asked to run, simulate, predict, explore what-if, or identify influential parameters.
license	MathWorks BSD-3-Clause
metadata	{"author":"MathWorks","version":"1.0"}

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 matlab-build-simbiology-model)
• Parameter estimation from data (use matlab-fit-simbiology-model)
• NCA / AUC / Cmax from data (use matlab-fit-simbiology-model)

Must-Follow Rules

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:

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:

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:

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:

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

Decision Table

Scenario	Approach
One-off simulation	sbiosimulate
Parameter sweep / Monte Carlo	createSimFunction
Dose/variant/parameter what-if	SimBiology.Scenarios + createSimFunction
Low molecule count / noise	SSA solver (cs.SolverType = 'ssa')
Which parameters matter?	sbiosobol (Sobol) or sbioelementaryeffects (Morris)
Quick sensitivity check	Local sensitivity via configset

Basic Simulation (sbiosimulate)

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

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

With a dose:

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

Plotting Results

Use sbioplot for quick visualization of SimData:

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

For custom plots, extract numeric data first:

[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:

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:

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

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

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

Repeated Simulation (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();

• Compiles once, runs many — much faster than sbiosimulate in a loop
• 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):

% 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)

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

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)

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

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:

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:

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

Content Type	Example
Dose vector	SimBiology.Scenarios('DoseLevel', [d1, d2])
Variant vector	SimBiology.Scenarios('Pop', [v1, v2])
Parameter values	SimBiology.Scenarios('ke', [0.05 0.1 0.2])
Species values	SimBiology.Scenarios('Drug', [50 100 200])
Probability distribution	add(sc, 'elementwise', 'ke', makedist('Lognormal',...), 'Number', 50)

Virtual population via distribution sampling

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

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

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

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

Ensemble (multiple trajectories)

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

Gene expression template (all 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'};

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

Sensitivity Analysis

Sobol (global, quantitative)

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 (R2025b+)
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)

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

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

Normalization	Meaning
'None'	Raw dY/dp
'Half'	(p/y) dY/dp
'Full'	Dimensionless; both sides normalized

Conventions

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

  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
• 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.

---