Modeling the Population Pharmacokinetics of Phenobarbital in Neonates

This example shows how to build 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™.

Load the Data

These data were downloaded from the website http://depts.washington.edu/rfpk/ (no longer active) 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.
t.plottitle = 'States versus Time';
t.updateFigureForPrinting();

Describe the Data

In order to perform nonlinear mixed-effects modeling on this dataset, we need to convert the data to a groupedData object, a container for holding tabular data that is divided into groups. The groupedData constructor automatically identifies commonly use variable names as the grouping variable or the 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
ans = 

                Description: ''
       VariableDescriptions: {}
              VariableUnits: {}
             DimensionNames: {'Observations'  'Variables'}
                   UserData: []
                   RowNames: {}
              VariableNames: {'ID'  'TIME'  'DOSE'  'WEIGHT'  'APGAR'  'CONC'}
          GroupVariableName: 'ID'
    IndependentVariableName: 'TIME'

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', 'DosingType', 'Bolus', 'EliminationType', ...
    'Linear-Clearance', 'HasResponseVariable', true);
model = pkmd.construct;

% The model species |Drug_Central| corresponds to the data variable |CONC|.
responseMap = 'Drug_Central = CONC';

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). NLMEFIT calculates fixed and random effects for each parameter. The underlying algorithm assumes parameters are normally distributed. This assumption may not hold for biological parameters that are constrained to be positive, such as volume and clearance. We need to specify a transform for the estimated parameters, such that the transformed parameters do follow a normal distribution. By default, SimBiology® chooses a log transform for all estimated parameters. Therefore, the model is:

$$log(V_i)=log(\phi_{V,i})=\theta_V+\eta_{V,i}$$

and

$$log(Cl_i)=log(\phi_{Cl,i})=\theta_{Cl}+\eta_{Cl,i},$$

where $\theta$, $eta$, and $\phi$ are the fixed effects, random effects, and estimated parameter values respectively, calculated for each group $i$. We provide some arbitrary initial estimates for V and Cl in the absence of better empirical data.

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

Extract the Dosing Information from the Data

Create a sample dose that targets species Drug_Central and use the createDoses method to generate doses for each infant based on the dosing amount listed in variable DOSE.

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

Fit the Model

Slightly loosen the tolerances to speed up the fit.

fitOptions.Options = statset('TolFun', 1e-3, 'TolX', 1e-3);
[nlmeResults, simI, simP] = sbiofitmixed(model, data, responseMap, ...
    estimatedParams, doses, 'nlmefit', fitOptions);

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, [], '', 'Drug_Central');
t.plot(simI, [], '', 'Drug_Central',...
    'legend',{'Observed', 'Fit-Pop.', 'Fit-Indiv.'});

Examine Fitted Parameters and Covariances

disp('Summary of initial results');
disp('Parameter Estimates Without Random Effects:');
disp(nlmeResults.PopulationParameterEstimates(1:2,:));
disp('Estimated Fixed Effects:');
disp(nlmeResults.FixedEffects);
disp('Estimated Covariance Matrix of Random Effects:');
disp(nlmeResults.RandomEffectCovarianceMatrix);
Summary of initial results
Parameter Estimates Without Random Effects:
    Group        Name        Estimate
    _____    ____________    ________

    1        'Central'         1.4086
    1        'Cl_Central'    0.006137

Estimated Fixed Effects:
      Name      Description     Estimate    StandardError
    ________    ____________    ________    _____________

    'theta1'    'Central'       0.34257       0.08627    
    'theta2'    'Cl_Central'    -5.0934     0.0004879    

Estimated Covariance Matrix of Random Effects:
             eta1       eta2  
            _______    _______

    eta1    0.20323          0
    eta2          0    0.19338

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 a box plot of the calculated random effects.
boxplot(nlmeResults);

Generate a Plot of the Residuals over Time

% 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 plot
% the residuals over time.
timeVec = data.(data.Properties.IndependentVariableName);
obsData = data.CONC;
indicesToKeep = ~isnan(obsData);
timeVec = timeVec(indicesToKeep);

% Get the residuals from the fitting results.
indRes = nlmeResults.stats.ires(indicesToKeep);
popRes = nlmeResults.stats.pres(indicesToKeep);

% 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);
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 mean value of each covariate for each group.

covariateLabels = {'WEIGHT', 'APGAR'};
covariates = grpstats(ds, data.Properties.GroupVariableName, 'mean',...
    'DataVars', covariateLabels);

Examine NLME Parameter Fit Results for Possible Covariate Dependencies

% Get the parameter values for each group (empirical Bayes estimates).
paramValues = nlmeResults.IndividualParameterEstimates.Estimate;
isCentral = strcmp('Central', nlmeResults.IndividualParameterEstimates.Name);
isCl = strcmp('Cl_Central', nlmeResults.IndividualParameterEstimates.Name);

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

subplot(2,2,3);
plot(covariates.mean_WEIGHT,paramValues(isCl), '.');
ylabel('Clearance');
xlabel('weight');

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

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

Create a CovariateModel 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 ("Central") varying with weight and clearance ("Cl_Central") varying with weight. Our model is:

$$log(V_i)=log(\phi_{V,i})=\theta_V+\theta_{V/weight}*weight_i+\eta_{V,i}$$

and

$$log(Cl_i)=log(\phi_{Cl,i})=\theta_{Cl}+\theta_{Cl/weight}*weight_i+\eta_{Cl,i}$$

% Define the covariate model.
covmodel            = CovariateModel;
covmodel.Expression = ({'Central = exp(theta1 + theta2*WEIGHT + eta1)','Cl_Central = exp(theta3 + theta4*WEIGHT + eta2)'});

% Use constructDefaultInitialEstimate to create a initialEstimates struct.
initialEstimates = covmodel.constructDefaultFixedEffectValues;
disp('Fixed Effects Description:');
disp(covmodel.FixedEffectDescription);
Fixed Effects Description:
    'Central'
    'Cl_Central'
    'Central/WEIGHT'
    'Cl_Central/WEIGHT'

Update the initial estimates using the values estimated from fitting the base model.

initialEstimates.theta1 = nlmeResults.FixedEffects.Estimate(1);
initialEstimates.theta3 = nlmeResults.FixedEffects.Estimate(2);
covmodel.FixedEffectValues = initialEstimates;

Fit the Model

[nlmeResults_cov, simI_cov, simP_cov] = sbiofitmixed(model, data, responseMap, ...
    covmodel, doses, 'nlmefit', fitOptions);

Examine Fitted Parameters and Covariances

disp('Summary of results when modeling covariate dependencies');
disp('Parameter Estimates Without Random Effects:');
disp(nlmeResults_cov.PopulationParameterEstimates);
disp('Estimated Fixed Effects:');
disp(nlmeResults_cov.FixedEffects);
disp('Estimated Covariance Matrix:');
disp(nlmeResults_cov.RandomEffectCovarianceMatrix);
Summary of results when modeling covariate dependencies
Parameter Estimates Without Random Effects:
    Group        Name        Estimate 
    _____    ____________    _________

    1        'Central'          1.2973
    1        'Cl_Central'    0.0061576
    2        'Central'          1.3682
    2        'Cl_Central'    0.0065512
    3        'Central'          1.3682
    3        'Cl_Central'    0.0065512
    4        'Central'         0.99431
    4        'Cl_Central'    0.0045173
    5        'Central'          1.2973
    5        'Cl_Central'    0.0061576
    6        'Central'          1.1664
    6        'Cl_Central'      0.00544
    7        'Central'          1.0486
    7        'Cl_Central'     0.004806
    8        'Central'          1.1664
    8        'Cl_Central'      0.00544
    9        'Central'          1.2973
    9        'Cl_Central'    0.0061576
    10       'Central'          1.2973
    10       'Cl_Central'    0.0061576
    11       'Central'          1.1664
    11       'Cl_Central'      0.00544
    12       'Central'          1.2301
    12       'Cl_Central'    0.0057877
    13       'Central'          1.1059
    13       'Cl_Central'    0.0051132
    14       'Central'          1.1059
    14       'Cl_Central'    0.0051132
    15       'Central'          1.2301
    15       'Cl_Central'    0.0057877
    16       'Central'          1.1664
    16       'Cl_Central'      0.00544
    17       'Central'          1.1059
    17       'Cl_Central'    0.0051132
    18       'Central'          1.0486
    18       'Cl_Central'     0.004806
    19       'Central'          1.0486
    19       'Cl_Central'     0.004806
    20       'Central'          1.1664
    20       'Cl_Central'      0.00544
    21       'Central'           1.605
    21       'Cl_Central'    0.0078894
    22       'Central'          1.3682
    22       'Cl_Central'    0.0065512
    23       'Central'          3.2052
    23       'Cl_Central'     0.017654
    24       'Central'          3.3803
    24       'Cl_Central'     0.018782
    25       'Central'         0.89394
    25       'Cl_Central'    0.0039908
    26       'Central'          3.9653
    26       'Cl_Central'     0.022619
    27       'Central'          1.6927
    27       'Cl_Central'    0.0083936
    28       'Central'          3.3803
    28       'Cl_Central'     0.018782
    29       'Central'          1.0486
    29       'Cl_Central'     0.004806
    30       'Central'           1.605
    30       'Cl_Central'    0.0078894
    31       'Central'          1.2973
    31       'Cl_Central'    0.0061576
    32       'Central'          4.1819
    32       'Cl_Central'     0.024064
    33       'Central'          1.5218
    33       'Cl_Central'    0.0074154
    34       'Central'          1.5218
    34       'Cl_Central'    0.0074154
    35       'Central'          2.3292
    35       'Cl_Central'     0.012173
    36       'Central'          1.3682
    36       'Cl_Central'    0.0065512
    37       'Central'          1.1664
    37       'Cl_Central'      0.00544
    38       'Central'          1.2301
    38       'Cl_Central'    0.0057877
    39       'Central'          1.6927
    39       'Cl_Central'    0.0083936
    40       'Central'          1.1059
    40       'Cl_Central'    0.0051132
    41       'Central'          1.5218
    41       'Cl_Central'    0.0074154
    42       'Central'          2.7323
    42       'Cl_Central'     0.014659
    43       'Central'         0.99431
    43       'Cl_Central'    0.0045173
    44       'Central'          1.2973
    44       'Cl_Central'    0.0061576
    45       'Central'         0.94279
    45       'Cl_Central'    0.0042459
    46       'Central'          1.1059
    46       'Cl_Central'    0.0051132
    47       'Central'          2.4565
    47       'Cl_Central'     0.012951
    48       'Central'         0.89394
    48       'Cl_Central'    0.0039908
    49       'Central'          1.2301
    49       'Cl_Central'    0.0057877
    50       'Central'          1.1059
    50       'Cl_Central'    0.0051132
    51       'Central'         0.99431
    51       'Cl_Central'    0.0045173
    52       'Central'         0.99431
    52       'Cl_Central'    0.0045173
    53       'Central'          1.5218
    53       'Cl_Central'    0.0074154
    54       'Central'           1.605
    54       'Cl_Central'    0.0078894
    55       'Central'          1.1059
    55       'Cl_Central'    0.0051132
    56       'Central'         0.84763
    56       'Cl_Central'     0.003751
    57       'Central'          1.8827
    57       'Cl_Central'    0.0095009
    58       'Central'          1.2973
    58       'Cl_Central'    0.0061576
    59       'Central'          1.1059
    59       'Cl_Central'    0.0051132

Estimated Fixed Effects:
      Name          Description        Estimate    StandardError
    ________    ___________________    ________    _____________

    'theta1'    'Central'              -0.48453      0.041857   
    'theta3'    'Cl_Central'            -5.9575    0.00031552   
    'theta2'    'Central/WEIGHT'        0.53203      0.040788   
    'theta4'    'Cl_Central/WEIGHT'     0.61957      0.074264   

Estimated Covariance Matrix:
              eta1       eta2  
            ________    _______

    eta1    0.029793          0
    eta2           0    0.04644

Visualise the Fitted Model with the Data

t.plot(simP_cov, [], '', 'Drug_Central');
t.plot(simI_cov, [], '', 'Drug_Central',...
    '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 shows that the population residuals are smaller in the covariate model fit compared to the original fit.

% Get the residuals from the fitting results.
indRes = nlmeResults_cov.stats.ires(indicesToKeep);
popRes = nlmeResults_cov.stats.pres(indicesToKeep);

% 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

Was this topic helpful?