SimBiology 3.1

Modeling the Population Pharmacokinetics of Phenobarbital in Neonates

This example builds a simple nonlinear mixed-effects model from clinical pharmacokinetic data.

Data were collected on 59 pre-term infants given phenobarbital for prevention of seizures during the first 16 days after birth. Each individual 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 individual at times other than dose times, as part of routine monitoring, for a total of 155 measurements. Infant weights and APGAR scores (a measure of newborn health) were also recorded.

This example uses data described in [1], a study funded by NIH/NIBIB grant P41-EB01975.

This example requires Statistics Toolbox™.

Contents

Load the Data

These data were downloaded from the website of the Resource Facility for Population Pharmacokinetics as a text file, and saved as a dataset array for ease of use.

load pheno.mat ds

Visualize the Data in a Trellis Plot

t = sbiotrellis(ds, 'ID', 'TIME', 'CONC', 'marker', 'o',...
       'markerfacecolor', [.7 .7 .7], 'markeredgecolor', 'r', ...
       'linestyle', 'none');

% Format the plot.
set(t.labels,'backgroundColor',[1 1 1]);
set(t.titleLabel, 'backgroundColor', [1 1 1], 'String', 'States versus Time'
);
set(t.hFig, 'color', [1 1 1]);

Describe the Data

In order to perform nonlinear mixed-effects modeling on this dataset, we need to add additional metadata to describe the role played in the modeling by particular variables.

pkd = PKData(ds);
pkd.GroupLabel          = 'ID';
pkd.IndependentVarLabel = 'TIME';
pkd.DependentVarLabel   = 'CONC';
pkd.DoseLabel           = 'DOSE';
pkd.CovariateLabels     = {'WEIGHT','APGAR'};

Define the Model

We will fit a simple one-compartment model, with bolus dose administration and linear clearance elimination, to the data.

pkmd = PKModelDesign;
pkmd.addCompartment('Central', 'Bolus', 'linear-clearance', true);
[model, modelMap] = pkmd.construct;

Specify Initial Estimates for the Parameters

The parameters fit in this model are the volume of the central compartment (Central) and the clearance rate (Cl_Central). We provide some arbitrary initial conditions in the absence of better empirical data.

beta0 = [1, 1];

Fit the Model

As nonlinear mixed-effects modeling can take up to several minutes, the results of the sbionlmefit command are provided for convenience in a MAT file. The command would be:

[nlmeResults simI simP] = sbionlmefit(model, modelMap, pkd, beta0);

load phenoresults.mat nlmeResults simI simP

Visualise the Fitted Model with the Data

Overlay the fitted simulation results on top of the observed data already displayed on the trellis plot. For the population results, simulations are run using the estimated fixed effects as the parameter values. For the individual results, simulations are run using the sum of the fixed and random effects as the parameter values.

t.plot(simP, [], '', modelMap.Observed);
t.plot(simI, [], '', modelMap.Observed,...
    'legend',{'Observed', 'Fit-Pop.', 'Fit-Indiv.'});

Examine Fitted Parameters and Covariances

The results are log-transformed.

beta = exp(nlmeResults.beta)
logbeta = nlmeResults.beta
sebeta = nlmeResults.stats.sebeta
psi = nlmeResults.psi

beta =

    1.4086
    0.0061


logbeta =

    0.3426
   -5.0934


sebeta =

    0.0612    0.0796


psi =

    0.2029         0
         0    0.1941

Generate a Box Plot of the Estimated Parameters

This example uses MATLAB® plotting commands to visualize the results. Several commonly used plots are also available as built-in options when performing parameter fits through the SimBiology® desktop interface.

% Create an appropriate name for each parameter.
names = {};
for i = 1:length(modelMap.Estimated)
    names{i} = ['log(' modelMap.Estimated{i} ')'];
end

% Create a box plot of the calculated random effects.
boxplot(nlmeResults.b',names);

% Label the plot.
xlabel('Random Effects');
ylabel('Values');
title('Estimates Box Plot');

Generate a Plot of the Residuals over Time

% Get the simulation data.
indSimData = selectbyname(simI,modelMap.Observed);
popSimData = selectbyname(simP,modelMap.Observed);

% Concatenate results from each group into a single vector.
indData = vertcat(indSimData.Data);
popData = vertcat(popSimData.Data);

% The simulations produce results only for the time points at which
% observations were made.  The vector of observation data also includes
% NaN's at the time points at which doses were given.  We need to remove
% these NaN's to be able to compare the vectors and calculate residuals.
ds = pkd.DataSet;
timeVec = ds.(pkd.IndependentVarLabel);
obsData = ds.(pkd.DependentVarLabel);
indicesToKeep = ~isnan(obsData);
timeVec = timeVec(indicesToKeep);
obsData = obsData(indicesToKeep);

% Calculate the residuals.
indRes = obsData - indData;
popRes = obsData - popData;

% Plot the residuals.  Get a handle to the plot to be able to add more data
% to it later.
resplot = figure;
plot(timeVec,indRes,'d','MarkerFaceColor','b','markerEdgeColor','b');
hold on;
plot(timeVec,popRes,'d','MarkerFaceColor','w','markerEdgeColor','b');
hold off;

% Create a reference line representing a zero residual, and set its
% properties to omit this line from the plot legend.
refl = refline(0,0);
set(get(get(refl,'Annotation'),'LegendInformation'),'IconDisplayStyle','off'
);

% Label the plot.
title('Residuals versus Time');
xlabel('Time');
ylabel('Residuals');
legend({'Individual','Population'});

Extract Group-dependent Covariate Values from the Dataset

% Get the number of groups and covariates.
nGroups     = numel(pkd.GroupNames);
nCovariates = numel(pkd.CovariateLabels);

% Get the mean value of each covariate for each group.
covariates = grpstats(pkd.DataSet, pkd.GroupLabel, 'mean',...
    'DataVars', pkd.CovariateLabels);

Examine NLME Parameter Fit Results for Possible Covariate Dependencies

% Calculate parameter values for each group by adding each group's random
% effects to the estimated fixed effects.
paramValues = bsxfun(@plus, nlmeResults.b, nlmeResults.beta)';

% Plot the parameter values vs. covariates for each group.
figure;subplot(2,2,1);
plot(covariates.mean_WEIGHT,paramValues(:,1), '.');
ylabel('log(V)');

subplot(2,2,3);
plot(covariates.mean_WEIGHT,paramValues(:,2), '.');
ylabel('log(Cl)');
xlabel('weight');

subplot(2,2,2);
plot(covariates.mean_APGAR,paramValues(:,1), '.');

subplot(2,2,4);
plot(covariates.mean_APGAR,paramValues(:,2), '.');
xlabel('APGAR');

Create a Group Design Matrix to Model the Covariate Dependencies

Based on the plots, there appears to be a correlation between volume and weight, clearance and weight, and possibly volume and APGAR score. We choose to focus on the effect of weight by modeling two of these covariate dependencies: volume varying with weight and clearance varying with weight. Our model is:

and

nParameters = 2;
nCovariateDependencies = 2;

% We need to set up the fixed effect group design matrices.  First, create
% a skeleton design matrix that is an nParameters-by-nParameters (2x2)
% identity matrix (representing the fixed effect for each individual
% parameter) concatenated with a placeholder column for each of the
% covariate dependencies.
A = eye(nParameters,nParameters+nCovariateDependencies);

% Next, copy this design matrix for each group.
FEGroupDesignMatrices = repmat(A,[1,1,nGroups]);

% Next, fill the placeholder columns with the covariate values for each
% group.  The row to fill for each column is based on the affected
% parameter.

% We need to find the order of estimated parameters to specify the correct
% parameter row indices.
modelMap.Estimated

% The first placeholder column models volume depending on weight.
dependency_col_index = 1;
parameter_row_index = 1;
FEGroupDesignMatrices(parameter_row_index, nParameters+dependency_col_index,
 1:nGroups) = ...
    covariates.mean_WEIGHT;

% The second placeholder column models clearance depending on weight.
dependency_col_index = 2;
parameter_row_index = 2;
FEGroupDesignMatrices(parameter_row_index, nParameters+dependency_col_index,
 1:nGroups) = ...
    covariates.mean_WEIGHT;

ans =

    'Central'
    'Cl_Central'

Update Initial Estimates for the Parameters and Covariate Dependency Factors

beta_covariatefactors = [0 0];
beta0 = [exp(nlmeResults.beta)' beta_covariatefactors];

Fit the Model

As nonlinear mixed-effects modeling can take up to several minutes, the results of the sbionlmefit command are provided for convenience in a MAT file. The commands would be:

options.FEGroupDesign = FEGroupDesignMatrices;
[nlmeResults_cov simI_cov simP_cov] = sbionlmefit(model, modelMap, pkd,...
   beta0, options);

load covariates_phenoresults.mat nlmeResults_cov simI_cov simP_cov

Examine Fitted Parameters and Covariances

The results are log-transformed.

beta = exp(nlmeResults_cov.beta)
logbeta = nlmeResults_cov.beta
sebeta = nlmeResults_cov.stats.sebeta
psi = nlmeResults_cov.psi

beta =

    0.6160
    0.0026
    1.7024
    1.8583


logbeta =

   -0.4845
   -5.9578
    0.5321
    0.6197


sebeta =

    0.0681    0.1219    0.0409    0.0742


psi =

    0.0300         0
         0    0.0466

Visualise the Fitted Model with the Data

t.plot(simP_cov, [], '', modelMap.Observed);
t.plot(simI_cov, [], '', modelMap.Observed,...
    'legend',{'Observed', 'Fit-Pop..', 'Fit-Indiv.', 'Cov. Fit-Pop.', '
Cov. Fit-Indiv.'});

Compare the Residuals to Those from a Model Without Covariate Dependencies

The following plot demonstrates that the population residuals are smaller in the covariate model fit compared to the original fit.

% Get the simulation data.
indSimData = selectbyname(simI_cov,modelMap.Observed);
popSimData = selectbyname(simP_cov,modelMap.Observed);

% Concatenate results from each group into a single vector.
indData = vertcat(indSimData.Data);
popData = vertcat(popSimData.Data);

% Calculate the residuals.
indRes = obsData - indData;
popRes = obsData - popData;

% Return to the original residual plot figure and add the new data.
figure(resplot);
hold on;
plot(timeVec,indRes,'d','MarkerFaceColor','r','markerEdgeColor','r');
plot(timeVec,popRes,'d','MarkerFaceColor','w','markerEdgeColor','r');
hold off;

% Update the legend.
legend('off');
legend({'Individual','Population','Individual(Cov.)','Population(Cov.)'
});

References

[1] Grasela TH Jr, Donn SM. Neonatal population pharmacokinetics of phenobarbital derived from routine clinical data. Dev Pharmacol Ther 1985:8(6). 374-83. PubMed Abstract