Documentation 
The following steps show one of the workflows you can use at the command line to fit a PK model and estimate parameters:
Import data as shown in Importing Data.
Specify the structural model by creating a PK model as shown in Create a Pharmacokinetic Model Using the Command Line. Alternatively, if you have a SimBiology^{®} model that you want to use in fitting, see Prerequisites for Using Custom SimBiology Models in Data Fitting.
Classify the data set to use in fitting. See Specify and Classify the Data to Fit.
Specify the initial guesses for the parameters to be estimated, as shown in Set Initial Estimates.
Perform individual or population fits:
For individual fits:
(Optional) Specify an error model or weights. See Specify an Error Model.
(Optional) Set tolerances.
(Optional) Specify maximum iterations.
(Optional) Specify to pool the data, which simultaneously fits data from multiple dose levels using the same model parameters for each dose.
For population fits:
Specify the statistical model:
Specify the covariate model and the covariance matrix. See Specify a Covariate Model and Specify the Covariance Pattern of Random Effects.
(Optional) Specify the error model. See Specify an Error Model.
(Optional) Set tolerances.
(Optional) Specify maximum iterations.
Obtain and visualize results.
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.
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:
Create the PKData object for the data set data.
pkDataObject = PKData(data);
PKData assigns the data set data to the readonly DataSet property.
Use the column headers in the data set to specify the following properties for the column in the data set.
Column in Data Set Represents  PKData Object Property 

Group identification labels  GroupLabel 
Independent variable (For example, time)  IndependentVarLabel 
Dependent variable (For example, measured response)  DependentVarLabel 
Amount of dose given  DoseLabel 
Rate of infusion (when applicable). Data must contain rate (amount/time) and not infusion time.  RateLabel 
Covariates (For example, age, gender, weight)  CovariateLabels 
For example, for the tobramycin data set [1]:
pkDataObject.GroupLabel = 'ID'; pkDataObject.IndependentVarLabel = 'Time'; pkDataObject.DependentVarLabel = 'Response'; pkDataObject.DoseLabel = 'Dose'; pkDataObject.CovariateLabels = {'WT','HT','AGE', 'SEX', 'CLCR'};
Note: For the subset of data belonging to a single group (as defined by the column in your data set that represents group identification labels, which you map to the GroupLabel property), the software allows multiple observations made at the same time. If this is true for your data, be aware that:

Tip If dosing applies to more than one compartment in the model, specify the DoseLabel property as follows: pkDataObject.DoseLabel = {'Dose1', 'Dose2'}; Dose1 and Dose2 are names of columns containing dose information for compartments. A onetoone relationship must exist between the number and order of elements in the DoseLabel property and the Dosed property of the corresponding PKModelMap object. 
Tip If your model measures multiple responses, specify the DependentVarLabel property as follows: pkDataObject.DependentVarLabel = {'Response1', 'Response2'}; Response1 and Response2 are names of columns containing response measurements. A onetoone relationship must exist between the number and order of elements in the DependentVarLabel property and the Observed property of the corresponding PKModelMap object. 
When you assign a column containing group identification labels to the GroupLabel property, PKData sets these readonly 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.
If you specify a stochastic solver and options in the Configset object associated with your model, be aware that during fitting SimBiology temporarily changes:
SolverType property to the default solver of ode15s
SolverOptions property to the options last configured for a deterministic solver
Caution If your model includes active variants that specify alternate values for the parameters to estimate, the variants are ignored for those parameters during fitting. 
To set the initial estimates (or initial guesses) for the parameters with fixed effects to estimate, first identify the sequence of the parameters in the model by querying the PKModelMap object. Next, construct a vector, beta0, containing the initial conditions. For information about PKModelMap objects, see step 4 in Create a Pharmacokinetic Model Using the Command Line.
Query the Estimated property of the PKModelMap object:
PKModelMapObj.Estimated
MATLAB^{®} returns the sequence of the parameters to be estimated. For example:
ans = 'Central' 'Cl_Central'
Set the initial estimates for the parameters. For example:
beta0 = [10.0, 1.0];
For information on ...  See ... 

The parameters added to the model  
Default units for the above parameters  Unit Conversion for Imported Data 
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, mixedeffects (NLME) model for this data:
Define groupspecific model parameters φ_{i} as linear combinations of fixed effects β and random effects b_{i}.
Define response values y_{i} as a nonlinear function f of the parameters and groupspecific covariate variables X_{i}.
The model is:
$$\begin{array}{l}{\phi}_{i}={A}_{i}\theta +{B}_{i}{\eta}_{i},where\text{}{\eta}_{i}\sim N(0,\Psi )\\ {y}_{i}=f({\phi}_{i},{X}_{i})+{\epsilon}_{i}\\ Alternatively,\mathrm{log}{y}_{i}=\mathrm{log}f({\phi}_{i},{X}_{i})+{\epsilon}_{i}\end{array}$$
This formulation of the nonlinear, mixedeffects model uses the following notation:
φ_{i}  A vector of groupspecific model parameters 
θ  A vector of fixed effects, modeling population parameters 
η_{i}  A vector of multivariate, normally distributed, groupspecific, random effects 
A_{i}  A groupspecific design matrix for combining fixed effects 
B_{i}  A groupspecific design matrix for combining random effects 
X_{i}  A data matrix of groupspecific covariate values 
y_{i}  A data vector of groupspecific response values 
f  A general, realvalued function of φ_{i} and X_{i} 
ε_{i} 

Ψ  A covariance matrix for the random effects 
σ^{2}  The error variance, assumed to be constant across observations 
For example, consider a onecompartment model with firstorder dosing and linear clearance. The groupspecific parameters (φ) in the model are clearance (Cl), compartment volume (V), and absorption rate constant (k_{a}). From the model:
$$\left(\begin{array}{c}Cl\\ V\\ ka\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}\theta Cl\\ \theta V\\ {\theta}_{ka}\end{array}\right)+\text{\hspace{0.17em}}\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}\eta Cl\\ \eta V\\ {\eta}_{ka}\end{array}\right)$$
In SimBiology, B_{i} is an identity matrix. That is, sbionlmefit does not support the specification of a different randomeffects 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, mixedeffects 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. See MixedEffects Models in the Statistics Toolbox documentation for more information.
Construct a CovariateModel object to define the relationship between parameters and covariates. After constructing the object, modify the FixedEffectValues property of the object before using the object as an input argument to sbionlmefit or sbionlmefitsa, to estimate nonlinear mixed effects.
If the NLME model assumes a groupdependent covariate such as weight (w), the model becomes:
$$\left(\begin{array}{c}Cl\\ V\\ ka\end{array}\right)=\left(\begin{array}{cccc}1& 0& 0& wi\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right)\left(\begin{array}{c}\theta Cl\\ \theta V\\ \begin{array}{l}{\theta}_{ka}\\ {\theta}_{Cl/w}\end{array}\end{array}\right)+\text{\hspace{0.17em}}\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}\eta Cl\\ \eta V\\ {\eta}_{ka}\end{array}\right)$$
Thus, the parameter for clearance (Cl) for an individual is Cl_{i} = θ_{Cl} + θ_{Cl/w} * w_{i} + η_{Cli}.
Use the following procedure to specify a covariate model. If you are using the tobramycin data set, make sure you first complete the following procedures:
Use the CovariateModel constructor function to construct an empty CovariateModel object:
covModel = CovariateModel;
Set the Expression property of the object to define the relationship between parameters and covariates in the CovariateModel object, where Cl, v, and ka are parameters, weight is a covariate, theta1, theta2, theta3, and theta4 are fixed effects, and eta1, eta2, and eta3 are random effects.
covModel.Expression = {'Cl = exp(theta1 + theta4*weight + eta1)',... 'v = exp(theta2 + eta2)',... 'ka = exp(theta3 + eta3)'};
Display a list of the descriptions of the fixed effects (theta1 and theta2) in the CovariateModel object:
disp('Fixed Effects Descriptions:'); disp(covModel.FixedEffectDescription)
Your output appears as follows, where each string describes the role of a fixed effect in the expression equation:
Fixed Effects Descriptions: 'Cl' 'v' 'ka' 'Cl/weight'
Use the constructDefaultFixedEffectValues method of the CovariateModel object to create a structure containing the initial estimates for the fixed effects in the object. The initial estimates in this structure are set to a default of zero:
initialEstimates = covModel.constructDefaultFixedEffectValues
Your output appears as:
initialEstimates = theta1: 0 theta2: 0 theta3: 0 theta4: 0
Edit the initialEstimates structure to set the initial estimates of the fixed effects:
initialEstimates.theta1 = 1.408; initialEstimates.theta2 = 0.061; initialEstimates.theta3 = 0.31;
Use the modified initialEstimates structure to update the FixedEffectValues property of the CovariateModel object:
covModel.FixedEffectValues = initialEstimates;
Now covModel, the CovariateModel object, is ready to submit as an input argument to sbionlmefit or sbionlmefitsa.
By default, the function you use to perform population fits (nlmefit or nlmefitsa) 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:
$$\left(\begin{array}{ccc}1& 1& 0\\ 1& 1& 0\\ 0& 0& 1\end{array}\right)$$
Create an options struct with the specified covariance pattern:
options.CovPattern = [1, 1, 0; 1, 1, 0; 0, 0, 1];
Specify the arguments for sbionlmefit or sbionlmefitsa:
[results, simdataI, simdataP] = sbionlmefit(modelObj,... PKModelMapObj, pkDataObject, beta0, options)
If you are using the tobramycin data set [1], do the following:
Create an options struct with the specified covariance pattern:
options.CovPattern = [1, 0; 0, 1];
Specify the arguments for sbionlmefit:
[results, simdataI, simdataP] = sbionlmefit(modelObj,... PKModelMapObj, pkDataObject, beta0, options)
results = NLMEResults handle Properties: FixedEffects: [2x3 dataset] RandomEffects: [97x2 dataset] IndividualParameterEstimates: [97x2 dataset] PopulationParameterEstimates: [97x2 dataset] RandomEffectCovarianceMatrix: [2x2 dataset] EstimatedParameterNames: {2x1 cell} CovariateNames: {'WT' 'HT' 'AGE' 'SEX' 'CLCR'} FixedEffectsStruct: [1x1 struct] stats: [1x1 struct]
results.FixedEffects
ans = Description Estimate StandardError 'Central' 3.0478 0.064369 'Cl_Central' 1.3054 0.061095
For more information, see nlmefit or nlmefitsa in the Statistics Toolbox documentation.
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. Similarly, relatively small estimates for covariances among certain random effects suggest that a full covariance matrix is unnecessary. Since random effects are unobserved, Ψ must be estimated indirectly. Specifying a diagonal or blockdiagonal covariance pattern for Ψ can improve convergence and efficiency of the fitting algorithm.
You can specify any of the following error models when using the sbionlinfit, sbionlmefit, or sbionlmefitsa function. Each model defines the error using a standard normal (Gaussian) variable e, the function value f, and one or two parameters, a and b. The default error model is 'constant'.
'constant': y = f + a*e (default)
'proportional': y = f + b*abs(f)*e
'combined': y = f + (a+b*abs(f))*e
'exponential': y = f*exp(a*e), or equivalently log(y) = log(f) + a*e
To define an error model:
Create an optionStruct input argument and set the ErrorModel field to specify one of the above error models. For example:
optionStruct.ErrorModel = 'proportional';
Specify the optionStruct input argument for sbionlinfit, sbionlmefit, or sbionlmefitsa, as shown in Perform Individual Fitting or Perform Population Fitting.
See also or nlinfit, nlmefit, or nlmefitsa in the Statistics Toolbox documentation.
To specify parameter transformations, use the ParamTransform option in sbionlinfit, sbionlmefit and sbionlmefitsa. The ParamTransform option lets you specify either no transformation, or the log, probit, or logit transformation.
Note: Do not use the ParamTransform option to specify parameter transformations when providing a CovariateModel object to a fitting function. The CovariateModel object provides the parameter transformation. 
The underlying algorithm in nlmefit assumes that parameters follow a normal distribution. This assumption may not hold for biological parameters that are constrained to be positive, such as volume and clearance. You may specify a transformation function for the estimated parameters, so that the transformed parameters follow a normal distribution.
By default, the SimBiology fitting functions choose a log transform for all estimated parameters. Parameters that are constrained between the values 0 and 1, like absorption fraction, can be transformed by the probit or logit transformations described below.
The probit function is the inverse cumulative distribution function (CDF) associated with the standard normal distribution. To apply the probit transform to a variable x in MATLAB, use the Statistics Toolbox function norminv: t = norminv(x). To untransform a variable t, use the function normcdf: x = normcdf(t).
The logit function is the inverse of the sigmoid function. To apply the logit transform to a variable x in MATLAB, use the following expression: t = log(x)  log(1x). To untransform the variable t, use x = 1/(1+exp(t)).
For the ParamTransform option, specify a vector of values equal to the number of parameters to be estimated. The values must be one of the integer codes listed in nlmefitsa or nlmefit specifying the transformation for the corresponding value of the parameters to be estimated. For example
options.ParamTransform = [0 1 2];
Specify the arguments for sbionlmefit or sbionlmefitsa, as shown in Perform Population Fitting.
For individual fitting, see Perform Individual Fitting.
The sbionlmefit and sbionlmefitsa functions let you specify a SimBiology model that you want to use in fitting. These functions use the nlmefit and nlmefitsa functions from the Statistics Toolbox to fit data with both fixed and random sources of variation using nonlinear mixedeffects and return the estimates. Both nlmefitand nlmefitsa fit the model by maximizing an approximation to the marginal likelihood with random effects integrated out assuming the following:
Random effects are multivariate, normally distributed, and independent between groups.
Observation errors are independent, identically normally distributed, and independent of random effects.
(Optional) Set the tolerance or maximum iteration options. Use an options structure that is an input argument for sbionlmefit or sbionlmefitsa:
optionStruct.Options.TolX = 1.0E4; optionStruct.Options.TolFun = 1.0E4; optionStruct.Options.MaxIter = 200;
Specify the model object, the PKModelMap object, the PKData object, the PKCovariateModel object, a vector containing the initial estimates for the fixed effects, and the options:
[results, simdataI, simdataP] = sbionlmefit(modelObj,... PKModelMapObj, pkDataObject, CovariateModelObject, beta0, optionStruct);
Note: If your population fit uses multiple doses, make sure each element in the Dosed property of the PKModelMap object is unique. 
Note: In your PKData object, for each subset of data belonging to a single group (as defined in the data column specified by the GroupLabel property), the software allows multiple observations made at the same time. If this is true for your data, be aware that:

sbionlmefit and sbionlmefitsa return the following:
results, an object containing estimated values and other statistics. For more information, see the sbionlmefit and sbionlmefitsa reference pages.
simdataI, a SimData object containing the data from simulating the model using the estimated parameter values for individuals, which includes both the fixed and random effects.
simdataP, SimData object containing the data from simulating the model using the estimated parameter values for the population, which includes only the fixed effects.
Plot the data from the data set. For example, in the imported data set used for fitting, ds, ID, Time, and Response are the column headers for the columns containing group IDs, time, and the response variable, respectively.
p = sbiotrellis(ds, 'ID', 'Time', 'Response')
Use the plot method on the trellis plot object p, returned by sbiotrellis to overlay data, using default values for the second and third input arguments.
p.plot(simdataP, [], '', PKModelMapObj.Observed);
For a description of the results, see sbionlmefit in the SimBiology documentation.
For more information, see the following topics in the Statistics Toolbox documentation:
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 loglikelihood. This functionality is useful for large and complex models when you expect the time to return the results to be longer than a few minutes. Use the options structure that is an argument for the sbionlmefit function:
% 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:
The fitting function tries to maximize the loglikelihood. When the plot begins to display a flat line, this might indicate that maximization is complete. Try setting the maximum iterations to a lower number to reduce the number of iterations you need and improve performance. For information on how to set iteration options, see Perform Population Fitting.
Plots for the fixed effects (β) and the estimates for the diagonal elements of the covariance matrix for the random effects (Ψ), should show convergence. If you see oscillations, or jumps without accompanying improvements in the loglikelihood, the model may be overparameterized. Try the following:
Reduce the number of fixed effects.
Reduce the number of random effects.
Simplify the covariance matrix pattern of random effects.
See also sbiofitstatusplot in the SimBiology documentation.
When performing population fitting using nonlinear regression, you can simultaneously fit data from multiple dose levels by either:
Using sbionlmefit with a CovariateModel object input argument and omitting the random effect (eta) from the expressions in the CovariateModel object.
Using sbionlmefit with an InitEstimates input argument and setting the REParamsSelect field or namevalue pair input argument to a 1byn logical vector, with all entries set to false, where n equals the number of fixed effects.
The sbionlinfit function lets you specify a SimBiology model to fit using the nlinfit function (individual fit). The nlinfit function uses nonlinear least squares and returns parameter estimates, residuals, and the estimated coefficient covariance matrix.
(Optional) Specify an error model, set the tolerance, set the maximum iteration, or set the data pooling option, which lets you simultaneously fit data from multiple dose levels using the same model parameters for each dose. Use an options structure that is an input argument for sbionlinfit:
optionStruct.ErrorModel = 'proportional'; optionStruct.TolX = 1.0E8; optionStruct.TolFun = 1.0E8; optionStruct.MaxIter = 100; optionStruct.Pooled = true;
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, PKDataObj, beta0, optionStruct);
Note: If your individual fit uses multiple doses, ensure each element in the Dosed property of the PKModelMap object is unique. 
Note: In your PKData object, for each subset of data belonging to a single group (as defined in the data column specified by the GroupLabel property), the software allows multiple observations made at the same time. If this is true for your data, be aware that:

sbionlinfit returns the following:
A results array of objects, with each object containing the following for one group:
ParameterEstimates — A dataset array containing fitted coefficients and their standard errors.
CovarianceMatrix — Estimated covariance matrix for the fitted coefficients.
beta — Vector of scalars specifying the fitted coefficients in transformed space.
R — Residuals.
J — Jacobian of modelObject.
COVB — Estimated covariance matrix for the transformed coefficients.
mse — Scalar specifying the estimate of the error of the variance term.
errorparam — Estimated parameters of the error model or an empty array if you specified weights using the 'Weights' namevalue pair argument.
simdataI, a SimData object containing the data from simulating the model using the estimated parameter values, for individuals.
Plot the data from the data set. For example, in the imported data set (ds), ID, Time and Response are the column headers for the columns containing group IDs, time, and the response variable respectively.
p = sbiotrellis(ds, 'ID', 'Time', 'Response')
Use the plot method on the trellis plot object p, returned by sbiotrellis to overlay data, using default values for the second and third input arguments.
p.plot(simdataI, [], '', PKModelMapObj.Observed);
For more information, see Nonlinear Regression and nlinfit in the Statistics Toolbox documentation.