Matlab Fit Simbiology Model

Fit SimBiology Models

Estimate parameters from data using fitproblem, fit population models

with NLME, generate virtual patients, and compute NCA metrics.

When to Use

• "fit", "estimate", "calibrate" model parameters
• Parameter estimation from experimental/observed data
• Population PK/PD, NLME, mixed effects, inter-individual variability
• Virtual patients, virtual cohorts
• NCA, AUC, Cmax, Tmax, half-life, clearance
• Keywords: "fit", "estimate", "calibrate", "population", "NCA", "AUC"

When NOT to Use

• Model construction or diagram (use matlab-build-simbiology-model)
• Simulation without fitting (use matlab-simulate-simbiology-model)
• Sensitivity analysis (use matlab-simulate-simbiology-model)

Must-Follow Rules

1. Use fitproblem for parameter estimation

Always use fitproblem instead of calling sbiofit or sbiofitmixed

directly. fitproblem provides a unified, declarative interface:

prob = fitproblem;
prob.Model = model;
prob.Data = data;
prob.ResponseMap = "Species = DataColumn";
prob.Estimated = estimatedInfo({'param'}, 'Bounds', [lo hi]);
results = fit(prob);

Do NOT call sbiofit(model, data, ...) or sbiofitmixed(model, data, ...)

directly — their positional argument signatures are error-prone.

2. Fitting requires groupedData, NOT a plain table

Always wrap data:

data = groupedData(table(...));
data.Properties.IndependentVariableName = 'Time';

3. ResponseMap maps model outputs to data columns

Format is always "ModelOutput = DataColumnName":

% Single compartment — use species name on the left
prob.ResponseMap = "Drug = DrugConc";

% Multi-compartment — use qualified name to disambiguate
prob.ResponseMap = "Central.Drug = DrugConc";

% When species name matches data column name, still use the = format
prob.ResponseMap = "Drug = Drug";

Use the unqualified species name unless the same species name exists

in multiple compartments (then qualify with Compartment.Species).

4. Always set bounds

Prevent non-physical values (negative rates, etc.):

estimParams = estimatedInfo({'ke','ka'}, ...
    'InitialValue', [0.2, 1.0], ...
    'Bounds', [0.01 1; 0.1 5]);

5. Use log transform for rate constants

Parameters spanning orders of magnitude (clearances, rate constants)

benefit from log-transform estimation. Set .Transform after creation:

ei = estimatedInfo({'ke','ka'}, 'InitialValue', [0.1, 0.5], 'Bounds', [0.01 1; 0.1 5]);
ei(1).Transform = 'log';
ei(2).Transform = 'log';

Alternative: use 'log(param)' name syntax (equivalent result):

ei = estimatedInfo({'log(ke)','log(ka)'}, 'InitialValue', [0.1, 0.5], 'Bounds', [0.01 1; 0.1 5]);

Important: InitialValue and Bounds are always in the

untransformed (natural) domain. Do NOT pass log(value).

Available transforms: 'log', 'logit', 'probit'

Do NOT pass 'Transform' as a name-value pair to the estimatedInfo

constructor — it errors. Always set the .Transform property after.

6. Error models for population fitting

Choose the error model that matches the noise structure:

• 'constant' — absolute noise uniform
• 'proportional' — noise scales with magnitude (most PK data)
• 'combined' — both constant and proportional
• 'exponential' — log-normal residual

7. NCA requires sbioncaoptions object

Do not use name-value pairs. Column names are camelCase.

EVDose column uses NaN for non-dose rows.

Decision Table

| Scenario | Approach |

|----------|----------|

| Single subject or pooled fit | fitproblem with FitFunction="sbiofit" |

| Individual fits per subject | fitproblem with Pooled=false |

| Population NLME (IIV, random effects) | fitproblem with FitFunction="sbiofitmixed" |

| Model-independent PK metrics | sbionca |

fitproblem Workflow (Preferred)

Use fitproblem for all parameter estimation. It provides a unified,

declarative interface that replaces direct calls to sbiofit/sbiofitmixed:

% 1. Prepare data
data = groupedData(table(tSample, yData, 'VariableNames', {'Time','Drug'}));
data.Properties.IndependentVariableName = 'Time';

% 2. Define parameters with bounds
estimParams = estimatedInfo({'ke','ka'}, ...
    'InitialValue', [0.2, 1.0], ...
    'Bounds', [0.01 1; 0.1 5]);

% 3. Build the fit problem
prob = fitproblem;
prob.Model = model;
prob.Data = data;
prob.ResponseMap = "Drug = Drug";
prob.Estimated = estimParams;
prob.Doses = dose;                  % optional
prob.FunctionName = 'scattersearch';
prob.ProgressPlot = true;           % show live progress

% 4. Fit
results = fit(prob);

% 5. Inspect
disp(results.ParameterEstimates);
plot(results);

Key fitproblem properties

| Property | Purpose |

|----------|---------|

| Model | The SimBiology model object |

| Data | groupedData table |

| Estimated | estimatedInfo object (not EstimatedParameters) |

| ResponseMap | Maps model species to data columns |

| Doses | Dose object(s) (not Dose) |

| FitFunction | "sbiofit" (default) or "sbiofitmixed" |

| FunctionName | Algorithm: 'scattersearch', 'nlinfit', 'fminsearch', 'lsqnonlin', 'particleswarm' |

| ProgressPlot | true to show live fitting progress |

| UseParallel | true for parallel evaluation |

| Pooled | true/false/"auto" (sbiofit only) |

| ErrorModel | "constant", "proportional", "combined", "exponential" |

| Variants | Variants to apply during fitting |

Common property name mistakes: prob.Estimated (not EstimatedParameters),

prob.Doses (not Dose), prob.FunctionName (not Algorithm or Method).

Estimation algorithms

| Method | Use Case |

|--------|----------|

| 'scattersearch' | Built-in global search, no extra toolbox — start here |

| 'nlinfit' | Default local; smooth problems |

| 'lsqnonlin' | Bounded least squares (Optimization Toolbox) |

| 'fminsearch' | Derivative-free, simple problems |

| 'particleswarm' | Global search (Global Optimization Toolbox) |

Dosing from multi-subject data

When subjects receive different doses, use createDoses to extract

per-subject dose objects from the data. The dose column must have NaN

on non-dosing rows:

% Data format: dose amount only at administration time, NaN elsewhere
%   ID  Time  Dose  DrugConc  Group
%   1   0     50    0         LowDose
%   1   1     NaN   2.05      LowDose
%   ...
%   3   0     200   0         HighDose

% Create template dose targeting the depot species
tempDose = sbiodose('StudyDose');
tempDose.TargetName = 'Depot.Drug';   % match your model's dose target

% Extract per-subject doses from groupedData
doseArray = createDoses(gData, 'Dose', '', tempDose);

% Pass to fitproblem
prob.Doses = doseArray;

Critical: If all rows have the dose value (not just dosing times),

createDoses will treat every row as a dose event. Use NaN on

non-dosing rows.

Population fitting (pooled vs individual)

data.Properties.GroupVariableName = 'SubjectID';

% Pooled — one parameter set for all
prob.Pooled = true;

% Individual — separate per subject
prob.Pooled = false;

Category-based pooling (per-group estimates)

To estimate parameters separately per category (e.g., dose group), use

CategoryVariableName on the estimatedInfo object — not on

fitproblem or sbiofit:

estimParams = estimatedInfo({'ke'}, 'InitialValue', 0.1, 'Bounds', [0.01 1]);
estimParams.CategoryVariableName = 'DoseGroup';  % column in data table
% Do NOT set prob.Pooled — leave it at the default

Warning: Do NOT set prob.Pooled when using CategoryVariableName.

Setting Pooled=false triggers per-subject individual fitting that

ignores CategoryVariableName (MATLAB issues a warning). Leave

Pooled unset to let the category-based pooling work correctly.

NLME Population Fitting

For inter-individual variability and random effects estimation,

set FitFunction to "sbiofitmixed":

% 1. Load & tag grouped data
data = groupedData(readtable('pop_pk_data.csv'));
data.Properties.IndependentVariableName = 'Time';
data.Properties.GroupVariableName = 'SubjectID';

% 2. Define parameters (Bounds ignored by sbiofitmixed — use InitialValue only)
estimParams = estimatedInfo({'CL','Vd','ka'}, ...
    'InitialValue', [5, 50, 1.2]);

% 3. Build the fit problem
prob = fitproblem;
prob.Model = model;
prob.Data = data;
prob.ResponseMap = "DrugConc = Concentration";
prob.Estimated = estimParams;
prob.FitFunction = "sbiofitmixed";
prob.ErrorModel = "proportional";
prob.ProgressPlot = true;

% 4. Fit
results = fit(prob);

% 5. Inspect
results.FixedEffects
results.RandomEffectCovarianceMatrix
results.IndividualParameterEstimates

When to use NLME vs sbiofit

| Criterion | FitFunction="sbiofit" | FitFunction="sbiofitmixed" |

|-----------|-----------|----------------|

| Single subject | Yes | |

| Multiple subjects, no IIV | Yes (pooled) | |

| Inter-individual variability | | Yes |

| Random effects estimation | | Yes |

| Covariate modeling | | Yes |

| Small datasets (< 5 subjects) | Yes | May not converge |

| Bounds on parameters | Yes (enforced) | Ignored — use good InitialValue instead |

NLME with covariates (CovariateModel)

When covariates (e.g., weight, age) influence parameters, use a

CovariateModel instead of estimatedInfo:

covModel = CovariateModel;
covModel.Expression = {
    'CL = theta1 + theta2*WT + eta1'
    'Vd = theta3 + theta4*WT + eta2'
    'ka = theta5 + eta3'
};
initVals = covModel.constructDefaultFixedEffectValues;
initVals.theta1 = 5; initVals.theta2 = 0.1;
initVals.theta3 = 50; initVals.theta4 = 0.5;
initVals.theta5 = 1.2;
covModel.FixedEffectValues = initVals;

prob = fitproblem;
prob.Model = model;
prob.Data = data;  % groupedData with WT column
prob.ResponseMap = "DrugConc = Concentration";
prob.FitFunction = "sbiofitmixed";
prob.Estimated = covModel;
prob.ErrorModel = "proportional";
results = fit(prob);

When to use which:

• estimatedInfo — NLME without covariates (simpler, fewer parameters)
• CovariateModel — NLME with covariates (parameter-covariate relationships)

Expression rules: theta prefix for fixed effects, eta for random

effects. One random effect max per expression. Use verify(covModel) to

validate syntax before fitting.

Virtual Patient Generation

From assumed distributions (Scenarios)

Use SimBiology.Scenarios with makedist — avoids manual matrix construction:

sc = SimBiology.Scenarios;
add(sc, 'elementwise', 'ke', makedist('Lognormal', 'mu', log(0.1), 'sigma', 0.3), 'Number', 100);
add(sc, 'elementwise', 'ka', makedist('Lognormal', 'mu', log(0.5), 'sigma', 0.25), 'Number', 100);

simfun = createSimFunction(model, sc, {'Drug'}, []);
results = simfun(sc, 24);

From NLME results (sbiosampleparameters)

Use sbiosampleparameters to sample from fitted population parameters —

it respects the covariate model parameterization automatically:

% Extract from NLME results
covModel = covariateModel(nlmeResults);
thetas = nlmeResults.FixedEffects;
omega = nlmeResults.RandomEffectCovarianceMatrix;

% Sample 200 virtual patients
nVP = 200;
vpParams = sbiosampleparameters(covModel.Expression, thetas, omega, nVP);

% Simulate
simfun = createSimFunction(model, {'CL','Vd','ka'}, {'Cp'}, []);
vpSim = simfun(vpParams, 48);

Non-Compartmental Analysis (NCA)

From simulation output

Use explicit OutputTimes to ensure sufficient time-resolution for NCA

(the default solver output may have too few points near Cmax):

cs = getconfigset(m, 'active');
cs.SolverOptions.OutputTimes = linspace(0, 24, 200);

[t, x, names] = sbiosimulate(m);
drugIdx = find(strcmp(names, 'Drug'));
Vd = sbioselect(m, 'Type', 'parameter', 'Name', 'Vd');
conc = x(:, drugIdx) ./ Vd.Value;
evDose = NaN(size(t)); evDose(1) = 100;
data = table(t, conc, evDose, 'VariableNames', {'Time','Concentration','EVDose'});

opt = sbioncaoptions;
opt.concentrationColumnName = 'Concentration';
opt.timeColumnName = 'Time';
opt.EVDoseColumnName = 'EVDose';
opt.AdministrationRoute = 'ExtraVascular';
ncaResults = sbionca(data, opt);

Administration routes

| Route | Dose column | Extra config |

|-------|-------------|--------------|

| 'ExtraVascular' | opt.EVDoseColumnName | — |

| 'IVBolus' | opt.IVDoseColumnName | — |

| 'IVInfusion' | opt.IVDoseColumnName | opt.infusionRateColumnName |

Key NCA metrics

All metric names use underscores (e.g., C_max not Cmax):

| Metric | Description |

|--------|-------------|

| AUC_0_last | Area under curve (0 to last time) |

| AUC_infinity | AUC extrapolated to infinity |

| C_max | Maximum observed concentration |

| T_max | Time of Cmax |

| T_half | Terminal elimination half-life |

| CL | Clearance (dose / AUC) |

| V_z | Volume of distribution (terminal) |

| MRT | Mean residence time |

Multi-subject NCA

data.Properties.GroupVariableName = 'SubjectID';
opt.groupColumnName = 'SubjectID';
ncaResults = sbionca(data, opt);

Confidence Intervals and Profile Likelihood

Restriction: sbioparameterci only works with results from nonlinear

regression (sbiofit). It does NOT support NLME results (sbiofitmixed).

After fitting with sbiofit, compute confidence intervals:

Gaussian (asymptotic) CI — fast, default

ciResults = sbioparameterci(fitResults);
disp(ciResults.Results);  % table: Name (cell), Estimate, Bounds, ConfidenceInterval (Nx2 double), Status (categorical)
plot(ciResults);

Column types in .Results table:

• Name — cell array of char (Results.Name{i})
• ConfidenceInterval — Nx2 double matrix (Results.ConfidenceInterval(i,:))
• Status — categorical (Results.Status(i), NOT {i})

Profile likelihood CI — more accurate for nonlinear models

ciPL = sbioparameterci(fitResults, 'Type', 'ProfileLikelihood');
plot(ciPL);  % shows profile likelihood curves with CI bounds
% Custom confidence level: 'Alpha', 0.10 for 90% CI

| Type | Speed | Use when |

|------|-------|----------|

| 'Gaussian' (default) | Fast | Quick check, well-behaved problems |

| 'ProfileLikelihood' | Slower | Final results, parameter identifiability |

Conventions

• Species names must differ from compartment names. Use distinct names: compartment Depot with species DrugDepot (not species Depot inside compartment Depot).
• Units on compartment volumes: Always specify units on compartment volumes (e.g., 'liter'). Set DimensionalAnalysis = true and VariableUnits on groupedData. For pure amount-based models, set Value = 1 and omit units.
• Loading .sbproj: proj = sbioloadproject('f.sbproj'); model = proj.(fieldnames(proj){1});
• Extracting data: selectbyname(sbiosimulate(m), 'Drug') for specific variables; resample(sd, tSample, 'linear') for specific times
• Start with 'scattersearch' if unsure about parameter landscape
• Use 'log' transform for parameters spanning orders of magnitude (see Rule 5)
• Do NOT call sbioaccelerate(model) before fitting — no effect, wastes time
• Set cs.MaximumWallClock = 60 — stops hung simulations from bad guesses
• Set prob.ProgressPlot = true for long-running fits

Evaluating Fit Quality

plot(results);                      % observed vs predicted overlay
plotResiduals(results);             % residuals vs time
plotResidualDistribution(results);  % histogram — should be ~normal
plotActualVersusPredicted(results); % identity line check

results.LogLikelihood  % higher = better
results.AIC            % lower = better (penalizes complexity)
results.BIC            % lower = better (stronger penalty)
results.MSE            % mean squared error

Model comparison: Compare by BIC — deltaBIC < -10 = strong evidence

for complex model; deltaBIC > 0 = simpler model preferred.

When to escalate to NLME: Multiple subjects with different parameter

values, systematic subject-specific residual patterns, need to quantify

inter-individual variability, or covariates may explain differences.

• Pass model objects (not UUID strings) to fitting functions
• Do NOT call sbiofit/sbiofitmixed/sbionlmefit directly — use fitproblem

References

Load on demand for detailed guidance:

• references/nca-analysis-guidance.md — full NCA patterns, IV infusion, metrics interpretation

----

Copyright 2026 The MathWorks, Inc.

----