%% Fit a One-Compartment PK Model to the Phenobarbital Data
% This example uses data collected on 59 preterm infants given phenobarbital during the first 16 days after birth [1]. Each infant received an initial dose followed by one or more sustaining doses by intravenous bolus administration. A total of between 1 and 6 concentration measurements were obtained from each infant at times other than dose times, for a total of 155 measurements. Infant weights and APGAR scores (a measure of newborn health) were also recorded.
% Load the data.

load pheno.mat ds

%
% Convert the dataset to a groupedData object, a container for holding tabular data that is divided into groups. It can automatically identify commonly used variable names as the grouping variable or independent (time) variable. Display the properties of the data and confirm that GroupVariableName and IndependentVariableName are correctly identified as 'ID' and 'TIME', respectively.

data = groupedData(ds);
data.Properties

%
% Create a simple one-compartment PK model with bolus dosing and linear clearance to fit such data. Use the PKModelDesign object to construct the model. Each compartment is defined by a name, dosing type, a clearance type, and whether or not the dosing requires a lag parameter. After constructing the model, you can also get a PKModelMap object map that lists the names of species and parameters in the model that are most relevant for fitting.

pkmd = PKModelDesign;
addCompartment(pkmd,'Central','DosingType','Bolus',...
                    'EliminationType','linear-clearance',...
                    'HasResponseVariable',true,'HasLag',false);
[onecomp, map] = pkmd.construct;

%
% Describe the experimentally measured response by mapping the appropriate model component to the response variable. In other words, indicate which species in the model corresponds to which response variable in the data. The PKModelMap property Observed indicates that the relevant species in the model is Drug_Central, which represents the drug concentration in the system. The relevant data variable is CONC, which you visualized previously.

map.Observed

%
% Map the Drug_Central species to the CONC variable.

responseMap = 'Drug_Central = CONC';

%
% The parameters to estimate in this model are the volume of the central compartment Central and the clearance rate Cl_Central. The PKModelMap property Estimated lists these relevant parameters. The underlying algorithm of sbiofit assumes parameters are normally distributed, but this assumption may not be true for biological parameters that are constrained to be positive, such as volume and clearance. Specify a log transform for the estimated parameters so that the transformed parameters follow a normal distribution. Use an estimatedInfo object to define such transforms and initial values (optional).

map.Estimated

%
% Define such estimated parameters, appropriate transformations, and initial values.

estimatedParams = estimatedInfo({'log(Central)','log(Cl_Central)'},'InitialValue',[1 1]);

%
% Each infant received a different schedule of dosing. The amount of drug is listed in the data variable DOSE. To specify these dosing during fitting, create dose objects from the data. These objects use the property TargetName to specify which species in the model receives the dose. In this example, the target species is Drug_Central, as listed by the PKModelMap property Dosed.

map.Dosed

%
% Create a sample dose with this target name and then use the createDoses method of groupedData object data to generate doses for each infant based on the dosing data DOSE.

sampleDose = sbiodose('sample','TargetName','Drug_Central');
doses = createDoses(data,'DOSE','',sampleDose);

%
% Fit the model.

[nlmeResults,simI,simP] = sbiofitmixed(onecomp,data,responseMap,estimatedParams,doses,'nlmefit');

%
% Visualize the fitted results using individual-specific parameter estimates.

plot(nlmeResults,'ParameterType','individual');

%
% Visualize the fitted results using population parameter estimates.

plot(nlmeResults,'ParameterType','population');

%
% Display the variation of estimated parameters using boxplot.

boxplot(nlmeResults)

%
% Compare the model predictions to the actual data.

plotActualVersusPredicted(nlmeResults)

%
% Plot the distribution of residuals.

plotResidualDistribution(nlmeResults)

%
% Plot residuals for each response using the model predictions on x-axis.

plotResiduals(nlmeResults,'Predictions')

% Copyright 2012 The MathWorks, Inc.