Simbiology Fit 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 simbiology-build-model)
• Simulation without fitting (use simbiology-simulate-model)
• Sensitivity analysis (use simbiology-simulate-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

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

Common property name mistakes: prob.Estimated (not EstimatedParameters), prob.Doses (not Dose), prob.FunctionName (not Algorithm or Method).

Estimation algorithms

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

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

Key NCA metrics

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

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

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.