Products & Services Solutions Academia Support User Community Company

Learn more about SimBiology   

Fitting Pharmacokinetic Model Parameters at the Command Line

On this page…

Fitting Parameters

Specifying and Classifying the Data to Fit

Setting Initial Estimates

Specifying the Covariance Pattern of Random Effects and a Covariate Model

Performing Population Fitting Using sbionlmefit

Performing Individual Fitting Using sbionlinfit

About Simulation Settings and Specifying Alternate Values for Initial Estimates

Fitting Parameters

The following steps show one of the workflows you can use at the command line to fit a PK model and estimate parameters:

  1. Import data as shown in Importing Data at the Command Line.

  2. Create a PK model as shown in Creating PK Models at the Command Line. Alternatively, if you have a SimBiology model that you want to use in fitting, see Parameter Fitting Using Custom SimBiology Models.

  3. (Optional) Simulate the model containing the dosing information to visualize the effects of dosing on the model as shown in Simulating a Model Containing Dosing Information.

  4. Specify the data set to use in fitting, and classify the data columns containing group identifiers, the independent variable (for example, time), the dependent variable (for example, the measured response), the dose, the rate of infusion (if applicable), and the columns that represent the covariates. See Specifying and Classifying the Data to Fit.

  5. Specify the initial guesses for the parameters to be estimated as shown in Setting Initial Estimates.

  6. Perform individual or population fits:

  7. Run the task and visualize results as shown in Visualizing Parameter Fitting Results and Generating Diagnostic Plots.

Back to Top

Specifying and Classifying the Data to Fit

In order to use the imported data in fitting, you must identify required columns in the data set that was previously imported as shown in Importing Data at the Command Line.

Use the PKData object to specify the data set containing the observed data to use in fitting. The properties of the PKData object specify what each column in the data represents.

To create the PKData object:

  1. Create the PKData object for the data set data.

    pkDataObject = PKData(data);

    PKData assigns the data set (data) to the read-only DataSet property.

  2. Use the column headers in the data set to specify the following properties for the column in the data set.

    Column in Data Set RepresentsPKData Object Property
    Group identification labelsGroupLabel
    Independent variable (for example, time)IndependentVarLabel
    Dependent variable (for example, measured response)DependentVarLabel
    Amount of dose givenDoseLabel
    Rate of infusion (when applicable). Data must contain rate (amount/time) and not infusion time.RateLabel
    Covariates (such as age, gender, weight)CovariateLabels

    For example:

    pkDataObject.GroupLabel          = 'ID'; 
pkDataObject.IndependentVarLabel = 'TIME';
pkDataObject.DependentVarLabel   = 'Y';
pkDataObject.DoseLabel           = 'DOSE';
pkDataObject.CovariateLabels     = {'Age', 'Wt'};

    When you assign a column containing group identification labels to the GroupLabel property, PKData sets the following read-only properties as follows:

    • The GroupNames property is set to the unique names found in the group column.

    • The GroupID property is set to an integer corresponding to the unique names found in the group column.

Back to Top

Setting Initial Estimates

To set the initial estimates (or initial guesses ) for the parameters with fixed effects that are to be estimated, first identify the sequence of the parameters is in the model by querying the PKModelMap object, and then construct a vector (beta0) containing the initial conditions. For information about PKModelMap, see step 4 in Creating PK Models at the Command Line),

  1. Query the Estimated property of the PKModelMap object.

    PKModelMapObj.Estimated

    MATLAB returns the sequence of the parameters to be estimate. For example:

    ans = 

    'Central'
    'Cl_Central'
    'ka'

  2. Set the initial guesses for the parameters. For example:

    beta0 = [1.5, 0.0398, 1.5669];
For information on ...See ...
The parameters added to the model

Default units for the above parametersUnit Conversion for Imported Data and the Model

Back to Top

Specifying the Covariance Pattern of Random Effects and a Covariate Model

Suppose data for a nonlinear regression model falls into one of m distinct groups i = 1, ..., m. (Specifically, suppose that the groups are not nested.) To specify a general nonlinear mixed-effects model for this data:

  1. Define group-specific model parameters φi as linear combinations of fixed effects β and random effects bi.

  2. Define response values yi as a nonlinear function f of the parameters and group-specific covariate variables Xi.

The model is:

This formulation of the nonlinear mixed-effects model uses the following notation:

φiA vector of group-specific model parameters
βA vector of fixed effects, modeling population parameters
biA vector of multivariate normally distributed group-specific random effects
AiA group-specific design matrix for combining fixed effects
BiA group-specific design matrix for combining random effects
XiA data matrix of group-specific covariate values
yiA data vector of group-specific response values
fA general, real-valued function of φi and Xi
εiA vector of group-specific errors, assumed to be independent, identically, normally distributed, and independent of bi
ΨA covariance matrix for the random effects
σ2The error variance, assumed to be constant across observations

For example, consider a one compartment model with first-order dosing and linear clearance. The group-specific parameters (φ) in the model are clearance (Cl), compartment volume (V), and absorption rate constant (ka) . From the model:

In SimBiology, Bi is an identity matrix. That is, sbionlmefit assumes the random effects are not correlated, and does not support the specification of a different random-effects design matrix.

You can alter the design matrices, as necessary, to introduce weighting of individual effects.

The Statistics Toolbox function nlmefit fits the general nonlinear mixed-effects model to data, estimating the fixed and random effects. The function also estimates the covariance matrix Ψ for the random effects. Additional diagnostic outputs allow you to assess tradeoffs between the number of model parameters and the goodness of fit.

Specifying the Covariance Matrix Pattern of Random Effects

By default, the function used to perform population fits (nlmefit) assumes a diagonal covariance matrix (no covariance among the random effects). To specify a different covariance pattern of random effects, use the 'CovPattern' option. In the previous example, assuming that each of the parameters has random effects and that Cl and V exhibit covariance, the covariance pattern of random effects would be a logical array:

  1. Create an options struct with the specified covariance pattern.

    options.CovPattern = [1, 1, 0; 1, 1, 0; 0, 0, 1];

  2. Specify the arguments for sbionlmefit.

    [results, simdataI, simdataP] = sbionlmefit(modelObj,...
				 PKModelMapObj, pkDataObject, beta0, options)

For more information, see nlmefit in the Statistics Toolbox User's Guide.

Fitting the model and estimating the covariance matrix Ψ often leads to further refinements. A relatively small estimate for the variance of a random effect suggests that it can be removed from the model. Likewise, relatively small estimates for covariances among certain random effects suggests that a full covariance matrix is unnecessary. Since random effects are unobserved, Ψ must be estimated indirectly. Specifying a diagonal or block-diagonal covariance pattern for Ψ can improve convergence and efficiency of the fitting algorithm.

Specifying the Covariate Model

If the model in the above example assumes a group-dependent covariate such as weight (w) the model becomes:

Thus, the parameter for clearance (Cl) for an individual is Cli = βCli + bCli + βCl/w * wi.

To specify a covariate model, use the 'FEGroupDesign' option.

'FEGroupDesign' is a p-by-q-by-m array specifying a different p-by-q fixed-effects design matrix for each of the m groups. Using the previous example the array would resemble the following:

  1. Create the array as follows:

    %Number of parameters in the model (Phi)
num_params = 3;
%Number of covariates
num_cov = 1;
%Assuming number of groups in the data set is 7
num_groups = 7;
% Array of covariate values
covariates = [75; 52; 66; 55; 70; 58; 62 ];
A = repmat(eye(num_params, num_params+num_cov), [1,1,num_groups]);
A(1,num_params+1,1:num_groups) = covariates(:,1)

  2. Create a struct with the specified design matrix.

    options.FEGroupDesign = A;

  3. Specify the arguments for sbionlmefit as shown in Performing Population Fitting Using sbionlmefit.

For more information, see Nonlinear Regression Models in the Statistics Toolbox User's Guide.

Back to Top

Performing Population Fitting Using sbionlmefit

The function sbionlmefit lets you specify a SimBiology model that you want to use in fitting. sbionlmefit uses the nlmefit function from the Statistics Toolbox to fit data with both fixed and random sources of variation using nonlinear mixed-effects and returns the estimates. nlmefit fits the model by maximizing an approximation to the marginal likelihood with random effects integrated out assuming the following:

  1. (Optional) Set the tolerance and maximum iteration options. Use an options structure that is an argument for sbionlmefit.

    options.Options.TolX = 1.0E-4;
options.Options.TolFun = 1.0E-4;
options.Options.MaxIter = 200;

  2. Specify the model object, the PKModelMap object, the PKData object, a vector containing the initial estimates for the fixed effects, and the options.

    [results, simdataI, simdataP] = sbionlmefit(modelObj,...
				 PKModelMapObj, pkDataObject, beta0, options);

    sbionlmefit returns the following:

    • A results structure containing

      • beta — Estimates for the fixed effects

      • psi — Estimated covariance matrix of the random effects

      • stats — Structure with statistics such as AIC, and BIC. For a complete list see nlmefit in the Statistics Toolbox User's Guide

      • b — Estimated random effects for each group in observedData

    • simdataI contains the data from simulating the model using the estimated parameter values for individuals. That is, the fixed plus the random effects.

    • simdataP contains the data from simulating the model using the estimated parameter values for the population. That is the fixed effects only.

  3. Plot the data from the data set used in fitting. For example, data is the imported data set used for fitting, ID, TIME and Y are the column headers for the column containing group IDs, time, and the response variable respectively. 

    p = sbiotrellis(data, 'ID', 'TIME', 'Y')

  4. Plot the simulation results on the same plot. You can leave the second and third arguments for sbiotrellis empty for this example since the second argument requires a function handle which is not applicable here, and the third argument requires the column for the x-axis which is by default time and thus need not be specified.

    p.plot = (simdataP,'', '', PKModelMapObj.Observed)

For a description of the results, see sbionlmefit in the SimBiology Reference.

For more information, see the following in the Statistics Toolbox User's Guide:

Obtaining the Status of Fitting

The sbiofitstatusplot function dynamically plots the progress of the fitting task. During the task, the function plots the fixed effects (β), the estimates for the diagonal elements of the covariance matrix for the random effects (Ψ), and the log-likelihood. This functionality is useful for large and complex models when time taken to return the results is expected to be longer than a few minutes. Use the options structure that is an argument for sbionlmefit.

% Create options structure with 'OutputFcn'.      
options.Options.OutputFcn = @sbiofitstatusplot;
% Pass options structure with OutputFcn to sbionlmefit function.
results = sbionlmefit(..., options);

The following figure shows the type of plots obtained.

Tips for interpreting status plots:

See also sbiofitstatusplot in the SimBiology Reference documentation.

Back to Top

Performing Individual Fitting Using sbionlinfit

The function sbionlinfit lets you specify a SimBiology model to fit using the nlinfit function (Individual fit). The nlinfit function fits using nonlinear least squares and returns parameter estimates, residuals, and the estimated coefficient covariance matrix.

  1. (Optional) Set the tolerance and maximum iteration options.

    options.TolX = 1.0E-8;
options.TolFun = 1.0E-8;
options.MaxIter = 100;

  2. Specify the model object, the PKModelMap object, the PKData object, a vector containing the initial estimates for the fixed effects, and the options.

    [results, simdataI] = sbionlinfit(modelobj,...
      PKModelMapObj, pkData, beta0, options);

    sbionlinfit returns the following:

    • A results array of structures, each containing the following for one group:

      • beta — Fitted coefficients

      • R — Residuals

      • J — Jacobian of modelObject

      • COVB — Estimated covariance matrix for the fitted coefficients

      • mse — Estimate of the error of the variance term

    • simdataI contains the data from simulating the model using the estimated parameter values, for individuals.

  3. Plot the data from the data set used in fitting. For example, data is the imported data set used for fitting, ID, TIME and Y are the column headers for the column containing group IDs, time, and the response variable respectively. 

    p = sbiotrellis(data, 'ID', 'TIME', 'Y')

  4. Plot the simulation results on the same plot. You can leave the second and third arguments for sbiotrellis empty for this example since the second argument requires a function handle which is not applicable here, and the third argument requires the column for the x-axis which is by default time and thus need not be specified.

    p.plot = (simdataI,'', '', PKModelMapObj.Observed)

For more information, see Nonlinear Regression Models and nlinfit in the Statistics Toolbox User's Guide.

Back to Top

About Simulation Settings and Specifying Alternate Values for Initial Estimates

Use the Variant object to store and apply alternate values for model components during a simulation. For more information about variants and how to add variants see Storing and Applying Alternate Model Values Using Variants.

Use the Configset object to change settings for simulations. For more information about performing simulations, see sbiosimulate. The model object created using the PKModelDesign object's construct method contains a default configuration set that uses the sundials solver.

If you change the solver type in the configuration set being used, during fitting there is a temporary change to sundials to support events in the model. The change is reversed after returning the results. If you want to change tolerance options for simulations select a deterministic solver and use the tolerance options provided in the deterministic solver. The fitting functions (sbionlmefit or sbionlinfit) will use the tolerance options specified in the deterministic solver.

Back to Top

  

Parameter Fitting in Pharmacokinetic ModelsParameter Fitting in Pharmacokinetic Models Fitting Pharmacokinetic Model Parameters in the SimBiology DesktopFitting Pharmacokinetic Model Parameters in the SimBiology Desktop

Recommended Products

Get MATLAB trial software

Includes the most popular MATLAB recorded presentations with Q&A sessions led by MATLAB experts.

Get the Interactive Kit

 © 1984-2009- The MathWorks, Inc.    -   Site Help   -   Patents   -   Trademarks   -   Privacy Policy   -   Preventing Piracy   -   RSS