Model 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 and Machine Learning 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.hFig.Color = [1 1 1];
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 = struct with fields:
Description: ''
UserData: []
DimensionNames: {'Observations' 'Variables'}
VariableNames: {'ID' 'TIME' 'DOSE' 'WEIGHT' 'APGAR' 'CONC'}
VariableTypes: ["double" "double" "double" "double" "double" "double"]
VariableDescriptions: {}
VariableUnits: {}
VariableContinuity: []
RowNames: {}
CustomProperties: [1x1 matlab.tabular.CustomProperties]
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:
and
where , , and are the fixed effects, random effects, and estimated parameter values respectively, calculated for each group . 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);
Visualize 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');
Summary of initial results
disp('Parameter Estimates Without Random Effects:');
Parameter Estimates Without Random Effects:
disp(nlmeResults.PopulationParameterEstimates(1:2,:));
Group Name Estimate _____ ______________ ________ 1 {'Central' } 1.4086 1 {'Cl_Central'} 0.006137
disp('Estimated Fixed Effects:');
Estimated Fixed Effects:
disp(nlmeResults.FixedEffects);
Name Description Estimate StandardError __________ ______________ ________ _____________ {'theta1'} {'Central' } 0.34257 0.061246 {'theta2'} {'Cl_Central'} -5.0934 0.079501
disp('Estimated Covariance Matrix of Random Effects:');
Estimated Covariance Matrix of Random Effects:
disp(nlmeResults.RandomEffectCovarianceMatrix);
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:
and
% 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 an initialEstimates struct. initialEstimates = covmodel.constructDefaultFixedEffectValues; disp('Fixed Effects Description:');
Fixed Effects Description:
disp(covmodel.FixedEffectDescription);
{'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');
Summary of results when modeling covariate dependencies
disp('Parameter Estimates Without Random Effects:');
Parameter Estimates Without Random Effects:
disp(nlmeResults_cov.PopulationParameterEstimates);
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
disp('Estimated Fixed Effects:');
Estimated Fixed Effects:
disp(nlmeResults_cov.FixedEffects);
Name Description Estimate StandardError __________ _____________________ ________ _____________ {'theta1'} {'Central' } -0.48453 0.067952 {'theta3'} {'Cl_Central' } -5.9575 0.12199 {'theta2'} {'Central/WEIGHT' } 0.53203 0.040788 {'theta4'} {'Cl_Central/WEIGHT'} 0.61957 0.074264
disp('Estimated Covariance Matrix:');
Estimated Covariance Matrix:
disp(nlmeResults_cov.RandomEffectCovarianceMatrix);
eta1 eta2 ________ _______ eta1 0.029793 0 eta2 0 0.04644
Visualize 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.