matlab-fit-simbiology-model

拟合 SimBiology 模型

使用 fitproblem 根据数据估计参数，使用 NLME 拟合群体模型， 生成虚拟患者并计算 NCA 指标。

适用场景

• “拟合”、“估计”或“校准”模型参数
• 根据实验数据／观测数据估计参数
• 群体 PK/PD、NLME、混合效应、个体间变异
• 虚拟患者、虚拟队列
• NCA、AUC、Cmax、Tmax、半衰期、清除率
• 关键词：“fit”、“estimate”、“calibrate”、“population”、“NCA”、“AUC”

不适用场景

• 构建模型或绘制模型图（使用 matlab-build-simbiology-model）
• 不进行拟合的模拟（使用 matlab-simulate-simbiology-model）
• 敏感性分析（使用 matlab-simulate-simbiology-model）

必须遵循的规则

1. 使用 fitproblem 进行参数估计

始终使用 fitproblem，不要直接调用 sbiofit 或 sbiofitmixed。 fitproblem 提供统一的声明式接口：

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

不要直接调用 sbiofit(model, data, ...) 或 sbiofitmixed(model, data, ...)——它们的位置参数签名容易出错。

2. 拟合需要使用 groupedData，而不是普通表格

始终对数据进行封装：

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

3. ResponseMap 将模型输出映射到数据列

格式始终为 "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";

除非同一物种名称存在于多个区室中（此时使用 Compartment.Species 限定），否则请使用未限定的物种名称。

4. 始终设置边界

防止出现非物理值（例如负速率）：

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

5. 对速率常数使用对数变换

跨越多个数量级的参数（清除率、速率常数）通常适合使用对数变换进行估计。创建后设置 .Transform：

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

或者，使用 'log(param)' 名称语法（结果等效）：

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

重要： InitialValue 和 Bounds 始终使用未变换（自然）域中的值。不要传入 log(value)。

可用变换：'log'、'logit'、'probit'

不要将 'Transform' 作为名称-值对传递给 estimatedInfo 构造函数，否则会报错。请始终在之后设置 .Transform 属性。

6. 群体拟合的误差模型

选择与噪声结构相匹配的误差模型：

• 'constant' — 绝对噪声均匀
• 'proportional' — 噪声随数值大小变化（大多数 PK 数据）
• 'combined' — 同时包含恒定和比例噪声
• 'exponential' — 对数正态残差

7. NCA 需要 sbioncaoptions 对象

不要使用名称-值对。列名采用 camelCase。 EVDose 列的非给药行使用 NaN。

决策表

场景	方法
单个受试者或合并拟合	使用 FitFunction="sbiofit" 的 fitproblem
按受试者分别拟合	使用 Pooled=false 的 fitproblem
群体 NLME（IIV、随机效应）	使用 FitFunction="sbiofitmixed" 的 fitproblem
与模型无关的 PK 指标	sbionca

fitproblem 工作流（推荐）

所有参数估计都使用 fitproblem。它提供统一的声明式接口，可替代直接调用 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);

关键 fitproblem 属性

属性	用途
Model	SimBiology 模型对象
Data	groupedData 表
Estimated	estimatedInfo 对象（不是 EstimatedParameters）
ResponseMap	将模型组分映射到数据列
Doses	剂量对象（不是 Dose）
FitFunction	"sbiofit"（默认）或 "sbiofitmixed"
FunctionName	算法：'scattersearch'、'nlinfit'、'fminsearch'、'lsqnonlin'、'particleswarm'
ProgressPlot	设为 true 可显示实时拟合进度
UseParallel	设为 true 可进行并行评估
Pooled	true/false/"auto"（仅限 sbiofit）
ErrorModel	"constant"、"proportional"、"combined"、"exponential"
Variants	拟合时应用的变体

常见属性名称错误： prob.Estimated（不是 EstimatedParameters）、 prob.Doses（不是 Dose）、prob.FunctionName（不是 Algorithm 或 Method）。

估计算法

方法	适用场景
'scattersearch'	内置全局搜索，无需额外工具箱——从这里开始
'nlinfit'	默认局部算法；适用于平滑问题
'lsqnonlin'	有界最小二乘（Optimization Toolbox）
'fminsearch'	无导数，适用于简单问题
'particleswarm'	全局搜索（Global Optimization Toolbox）

使用多受试者数据进行给药

当受试者接受不同剂量时，使用 createDoses 从数据中提取每位受试者的剂量对象。非给药行的剂量列必须为 NaN：

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

重要： 如果所有行都有剂量值（而不仅仅是给药时间）， createDoses 会将每一行都视为给药事件。非给药行请使用 NaN。

群体拟合（合并与个体）

data.Properties.GroupVariableName = 'SubjectID';

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

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

基于类别的合并（按组估计）

要按类别（例如剂量组）分别估计参数，请在 estimatedInfo 对象上使用 CategoryVariableName——不要在 fitproblem 或 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

警告： 使用 CategoryVariableName 时，不要设置 prob.Pooled。 设置 Pooled=false 会触发按受试者进行个体拟合，从而忽略 CategoryVariableName（MATLAB 会发出警告）。请勿设置 Pooled， 以便类别合并正常工作。

NLME 群体拟合

要估计个体间变异和随机效应， 请将 FitFunction 设为 "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

何时使用 NLME，何时使用 sbiofit

标准	FitFunction="sbiofit"	FitFunction="sbiofitmixed"
单个受试者	是
多个受试者，无 IIV	是（合并）
个体间变异		是
随机效应估计		是
协变量建模		是
小型数据集（< 5 名受试者）	是	可能无法收敛
参数边界	是（强制执行）	忽略——请改用合适的 InitialValue

使用协变量的 NLME（CovariateModel）

当协变量（例如体重、年龄）会影响参数时，请使用 CovariateModel，而不是 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);

如何选择：

• estimatedInfo — 不含协变量的 NLME（更简单，参数更少）
• CovariateModel — 含协变量的 NLME（参数与协变量之间的关系）

表达式规则： 固定效应使用 theta 前缀，随机效应使用 eta 前缀。 每个表达式最多包含一个随机效应。拟合前使用 verify(covModel) 验证语法。

虚拟患者生成

根据假定分布生成（Scenarios）

使用 SimBiology.Scenarios 和 makedist，避免手动构建矩阵：

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

根据 NLME 结果生成（sbiosampleparameters）

使用 sbiosampleparameters 从拟合的群体参数中抽样—— 它会自动遵循协变量模型的参数化方式：

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

非房室分析（NCA）

根据仿真输出进行分析

使用显式的 OutputTimes，确保 NCA 有足够的时间分辨率 （默认求解器输出在 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);

给药途径

路径	剂量列	额外配置
'ExtraVascular'	opt.EVDoseColumnName	—
'IVBolus'	opt.IVDoseColumnName	—
'IVInfusion'	opt.IVDoseColumnName	opt.infusionRateColumnName

关键 NCA 指标

所有指标名称均使用下划线（例如，C_max，而非 Cmax）：

指标	描述
AUC_0_last	曲线下面积（从 0 到最后一个时间点）
AUC_infinity	外推至无穷大的 AUC
C_max	观测到的最大浓度
T_max	Cmax 对应的时间
T_half	终末消除半衰期
CL	清除率（剂量 / AUC）
V_z	分布容积（终末）
MRT	平均滞留时间

多受试者 NCA

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

置信区间和剖面似然

限制：sbioparameterci 仅适用于非线性回归（sbiofit）的结果。不支持 NLME 结果（sbiofitmixed）。

使用 sbiofit 完成拟合后，计算置信区间：

高斯（渐近）置信区间——快速，默认

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

.Results 表中的列类型：

• Name — 字符数组的元胞数组（Results.Name{i}）
• ConfidenceInterval — Nx2 双精度矩阵（Results.ConfidenceInterval(i,:)）
• Status — 分类数据（Results.Status(i)，不是 {i}）

剖面似然置信区间——对非线性模型更准确

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

类型	速度	适用情形
'Gaussian'（默认）	快	快速检查、表现良好的问题
'ProfileLikelihood'	较慢	最终结果、参数可辨识性

约定

• **物种名称必须与区室名称不同。**使用不同的名称：区室 Depot 和物种 DrugDepot（而不是在区室 Depot 中使用物种 Depot）。
• **区室体积的单位：**始终为区室体积指定单位（例如，'liter'）。设置 DimensionalAnalysis = true，并在 groupedData 上设置 VariableUnits。对于纯基于量的模型，设置 Value = 1 并省略单位。
• 加载 .sbproj：proj = sbioloadproject('f.sbproj'); model = proj.(fieldnames(proj){1});
• **提取数据：**使用 selectbyname(sbiosimulate(m), 'Drug') 获取特定变量；使用 resample(sd, tSample, 'linear') 获取特定时间点
• 如果不确定参数空间的情况，请从 'scattersearch' 开始
• 对跨越多个数量级的参数使用 'log' 变换（参见规则 5）
• 不要在拟合前调用 sbioaccelerate(model)——这不起作用，只会浪费时间
• 设置 cs.MaximumWallClock = 60——避免错误猜测导致模拟挂起
• 对于长时间运行的拟合，设置 prob.ProgressPlot = true

评估拟合质量

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

**模型比较：**按 BIC 比较——deltaBIC < -10 表示复杂模型有强有力的证据支持；deltaBIC > 0 表示更倾向于简单模型。

**何时升级到 NLME：**多个受试者具有不同的参数值、存在系统性的受试者特异性残差模式、需要量化个体间变异，或协变量可能解释差异时。

• 向拟合函数传递模型对象（而不是 UUID 字符串）
• 不要直接调用 sbiofit/sbiofitmixed/sbionlmefit——请使用 fitproblem

参考资料

需要详细指导时按需加载：

• references/nca-analysis-guidance.md — 完整的 NCA 模式、IV 输注和指标解读

---

Copyright 2026 The MathWorks, Inc.

---