www.gusucode.com > simbio 案例源码程序 matlab代码 > simbio/FitAOneCompartmentModelToAnIndividualPKProfileExample.m
%% Fit a One-Compartment Model to an Individual's PK Profile
%% Background
% This example shows how to fit an individual's PK profile data to
% one-compartment model and estimate pharmacokinetic parameters.
%
% Suppose you have drug plasma concentration data from an individual and
% want to estimate the volume of the central compartment and the clearance.
% Assume the drug concentration versus the time profile follows the
% monoexponential decline ${C_t} = {C_0}{e^{ - {k_e}t}}$, where $C_t$ is
% the drug concentration at time t, $C_0$ is the initial concentration, and
% $k_e$ is the elimination rate constant that depends on the clearance and
% volume of the central compartment ${k_e} = {Cl/V}$.
%
% The synthetic data in this example was generated using the following
% model and parameters:
%
% * One-compartment model with bolus dosing and first-order elimination
% * Volume of the central compartment (|Central|) = 1.70 liter
% * Clearance parameter (|Cl_Central|) = 0.55 liter/hour
% * Constant error model
%
%% Load Data and Visualize
% The data is stored as a table with variables |Time| and |Conc| that
% represent the time course of the plasma concentration of an individual
% after an intravenous bolus administration measured at 13 different time
% points. The variable units for |Time| and |Conc| are hour and
% milligram/liter, respectively.
clear all
load(fullfile(matlabroot,'examples','simbio','data15.mat'))
%%
plot(data.Time,data.Conc,'b+')
xlabel('Time (hour)');
ylabel('Drug Concentration (milligram/liter)');
%% Convert to groupedData Format
% Convert the data set to a |groupedData| object, which is the required
% data format for the fitting function |sbiofit| for later use. A
% |groupedData| object also lets you set independent variable and group
% variable names (if they exist). Set the units of the |Time| and |Conc|
% variables. The units are optional and only required for the
% <docid:simbio_ref.bqfnt1l-1> feature, which automatically converts
% matching physical quantities to one consistent unit system.
gData = groupedData(data);
gData.Properties.VariableUnits = {'hour','milligram/liter'};
gData.Properties
%%%
% |groupedData| automatically set the name of the |IndependentVariableName|
% property to the |Time| variable of the data.
%% Construct a One-Compartment Model
% Use the built-in PK library to construct a one-compartment model with
% bolus dosing and first-order elimination where the elimination rate
% depends on the clearance and volume of the central compartment. Use the
% |configset| object to turn on unit conversion.
pkmd = PKModelDesign;
pkc1 = addCompartment(pkmd,'Central');
pkc1.DosingType = 'Bolus';
pkc1.EliminationType = 'linear-clearance';
pkc1.HasResponseVariable = true;
model = construct(pkmd);
configset = getconfigset(model);
configset.CompileOptions.UnitConversion = true;
%%%
% For details on creating compartmental PK models using the SimBiology(R)
% built-in library, see <docid:simbio_ug.bsehazy-47>.
%% Define Dosing
% Define a single bolus dose of 10 milligram given at time = 0. For details
% on setting up different dosing schedules, see
% <docid:simbio_ug.bsrr3oq-1>.
dose = sbiodose('dose');
dose.TargetName = 'Drug_Central';
dose.StartTime = 0;
dose.Amount = 10;
dose.AmountUnits = 'milligram';
dose.TimeUnits = 'hour';
%% Map Response Data to the Corresponding Model Component
% The data contains drug concentration data stored in the |Conc|
% variable. This data corresponds to the |Drug_Central| species in the model.
% Therefore, map the data to |Drug_Central| as follows.
responseMap = {'Drug_Central = Conc'};
%% Specify Parameters to Estimate
% The parameters to fit in this model are the volume of the central
% compartment (Central) and the clearance rate (Cl_Central). In this case,
% specify log-transformation for these biological parameters since they are
% constrained to be positive. The |estimatedInfo| object lets you specify
% parameter transforms, initial values, and parameter bounds if needed.
paramsToEstimate = {'log(Central)','log(Cl_Central)'};
estimatedParams = estimatedInfo(paramsToEstimate,'InitialValue',[1 1],'Bounds',[1 5;0.5 2]);
%% Estimate Parameters
% Now that you have defined one-compartment model, data to fit, mapped
% response data, parameters to estimate, and dosing, use |sbiofit| to
% estimate parameters. The default estimation function that |sbiofit| uses
% will change depending on which toolboxes are available. To see which
% function was used during fitting, check the |EstimationFunction| property
% of the corresponding results object.
fitConst = sbiofit(model,gData,responseMap,estimatedParams,dose);
%% Display Estimated Parameters and Plot Results
% Notice the parameter estimates were not far off from the true values
% (1.70 and 0.55) that were used to generate the data. You may also try
% different error models to see if they could further improve the parameter
% estimates.
fitConst.ParameterEstimates
s.Labels.XLabel = 'Time (hour)';
s.Labels.YLabel = 'Concentration (milligram/liter)';
plot(fitConst,'AxesStyle',s);
%% Use Different Error Models
% Try three other supported error models (proportional, combination of
% constant and proportional error models, and exponential).
fitProp = sbiofit(model,gData,responseMap,estimatedParams,dose,...
'ErrorModel','proportional');
fitExp = sbiofit(model,gData,responseMap,estimatedParams,dose,...
'ErrorModel','exponential');
fitComb = sbiofit(model,gData,responseMap,estimatedParams,dose,...
'ErrorModel','combined');
%% Compare Information Criteria for Model Selection
% Compare the loglikelihood, AIC, and BIC values of each model to see which
% eror model best fits the data. A larger likelihood value indicates the
% corresponding model fits the model better. For AIC and BIC, the smaller
% values are better.
allResults = [fitConst,fitProp,fitExp,fitComb];
errorModelNames = {'constant', 'proportional', 'exponential', 'combined'};
LogLikelihood = [allResults.LogLikelihood]';
AIC = [allResults.AIC]';
BIC = [allResults.BIC]';
t = table(LogLikelihood,AIC,BIC);
t.Properties.RowNames = errorModelNames;
t
%%
% Based on the information criteria, the constant error model fits the
% data best since it has the largest loglikelihood value and the smallest
% AIC and BIC.
%% Display Estimated Parameter Values
% Show the estimated parameter values of each error model.
allResults = [fitConst,fitProp,fitExp,fitComb];
Estimated_Central = zeros(4,1);
Estimated_Cl_Central = zeros(4,1);
t = table(Estimated_Central,Estimated_Cl_Central);
t.Properties.RowNames = errorModelNames;
for i = 1:height(t)
t{i,1} = allResults(i).ParameterEstimates.Estimate(1);
t{i,2} = allResults(i).ParameterEstimates.Estimate(2);
end
t
%% Conclusion
% This example showed how to estimate PK parameters, namely the volume of
% the central compartment and clearance parameter of an individual, by
% fitting the PK profile data to one-compartment model. You compared the
% information criteria of each model and estimated parameter values of
% different error models to see which model best explained the data. Final
% fitted results suggested both the constant and combined error models
% provided the closest estimates to the parameter values used to generate
% the data. However, the constant error model is a better model as
% indicated by the loglikelihood, AIC, and BIC information criteria.