Simulate and evaluate fitted SimBiology model
[ynew,parameterEstimates]=
predict(resultsObj)
[ynew,parameterEstimates]=
predict(resultsObj,data,dosing)
[
returns simulation results ynew
,parameterEstimates
]=
predict(resultsObj
)ynew
and parameter estimates
parameterEstimates
of a fitted SimBiology^{®} model.
[
returns simulation results ynew
,parameterEstimates
]=
predict(resultsObj
,data
,dosing
)ynew
and estimated parameter values
parameterEstimates
from evaluating the fitted SimBiology model
using the specified data
and dosing
information.
During simulations, predict
uses the parameter
values in the resultsObj.ParameterEstimates
property.
Use this method when you want to evaluate the fitted model and predict
responses using new data and/or dosing information.
resultsObj
— Estimation resultsOptimResults
object  NLINResults
objectEstimation results, specified as an OptimResults
object
or NLINResults object
,
which contains estimation results returned by sbiofit
.
It must be a scalar object.
data
— Output times or grouped datagroupedData
objectOutput times or grouped data, specified as a vector, or cell
array of vectors of output times, or groupedData
object
.
If it is a vector of time points, predict
simulates
the model with new time points using the parameter estimates from
the results object resultsObj
.
If it is a cell array of vectors of time points, predict
simulates
the model n times using the output times from each
time vector, where n is the length of data
.
If it is a groupedData
object, it must have
an independent variable such as Time. It must also have a group variable
if the training data used for fitting has such variable. You can use
a groupedData
object to query different combinations
of categories if the resultsObj
contains parameter
estimates for each category. predict
simulates
the model for each group with the specified categories. For instance,
suppose you have a set of parameter estimates for sex category (males
versus females), and age category (young versus old) in your training
data. You can use predict
to simulate the responses
of an old male (or any other combination) although such patient may
not exist in the training data.
If the resultsObj
is from estimating categoryspecific
parameters, data
must be a groupedData
object.
If UnitConversion
is turned on for the
underlying SimBiology model that was used for fitting and data
is
a groupedData
object, data
must
specify valid variable units via data.Properties.VariableUnits
property.
If it is a numeric vector or cell array of vectors of time points, predict
uses
the model’s TimeUnits
.
dosing
— Dosing information[]
 2D matrix of SimBiology dose objectsDosing information, specified as an empty array []
or 2D matrix of
SimBiology dose objects (ScheduleDose object
or RepeatDose object
). If dosing
is a matrix of dose
objects, the matrix must contain default doses or be consistent with the original dosing
data used with sbiofit
. That is, dose objects in
dosing
must have the same values for dose properties (such as
TargetName
) or must be parameterized in
the same way as the original dosing data. For instance, suppose that the original dosing
matrix has two columns of doses, where the doses in the first column target species
x and those in the second column target species y.
Then dosing
must have doses in the first column targeting species
x and those in the second column targeting species
y.
If empty, no doses are applied during simulation, even if the model has active doses.
If not empty, the matrix must have a single row or one row per group in the input data. If it has a single row, the same doses are applied to all groups during simulation. If it has multiple rows, each row is applied to a separate group, in the same order as the groups appear in the input data. If some groups (or time vectors) require more doses than others, then fill in the matrix with default (dummy) doses.
Multiple columns are allowed so that you can apply multiple dose objects to
each group or time vector. All doses in a column must be default doses or must
reference the same components in the model (for instance, the doses must have
the same dose target TargetName
) and must have consistent
parameterized properties as in the original dosing data used with
sbiofit
. For example, if the
Amount
property of a dose used in the original
sbiofit
call is parameterized to a modelscoped
parameter 'A'
, all doses for the corresponding group (column)
in dosing
must have the Amount
property parameterized to 'A'
.
A default dose has default values for all properties, except for the
Name
property. Create a default dose as
follows.
d1 = sbiodose('d1');
In addition to manually constructing dose objects using
sbiodose
, if the input data is a
groupedData
object and has dosing information, you can
use the createDoses
method to construct
doses from it.
The number of rows in the dosing
matrix and the number of groups or
output time vectors in data
determine the total number of simulation
results in the output ynew
. For details, see the table in the
ynew
argument description.
If UnitConversion
is turned on for the underlying SimBiology model that was used for fitting, dosing
must specify valid amount and time units.
ynew
— Simulation resultsSimData
objectsSimulation results, returned as a vector of SimData
objects. The states reported in ynew
are the states that were included in the responseMap
input argument of sbiofit
as well as any other states listed in the StatesToLog
property of the runtime options (RuntimeOptions
) of the SimBiology model.
The total number of simulation results in ynew
depends on the number of groups or output time vectors in data
and the number of rows in the dosing
matrix.
Number of groups or output time vectors in data  Number of rows in the dosing matrix  Simulation results 


 The total number of No doses are applied during simulation. 

 The total number of The given row of doses is applied during the simulation. 
 N  The total number of Each row of 
N 
 The total number of No doses are applied during simulation. 
N 
 The total number of The same row of doses is applied to each simulation. 
N  N  The total number of Each row of 
M  N  The function throws an error when M ≠ N. 
parameterEstimates
— Estimated parameter valuesEstimated parameter values, returned as a table. This is identical
to resultsObj.ParameterEstimates
property. The predict
method
uses these parameter values during simulation.
This example uses the yeast heterotrimeric G protein model and experimental data reported by [1]. For details about the model, see the Background section in Parameter Scanning, Parameter Estimation, and Sensitivity Analysis in the Yeast Heterotrimeric G Protein Cycle.
Load the G protein model.
sbioloadproject gprotein
Enter the experimental data containing the time course for the fraction of active G protein, as reported in the reference paper [1].
time = [0;10;30;60;110;210;300;450;600]; GaFracExpt = [0;0.35;0.4;0.36;0.39;0.33;0.24;0.17;0.2];
Create a groupedData
object based on
the experimental data.
tbl = table(time,GaFracExpt); grpData = groupedData(tbl);
Map the appropriate model component to the experimental
data. In other words, indicate which species in the model corresponds
to which response variable in the data. In this example, map the model
parameter GaFrac
to the experimental data variable GaFracExpt
from grpData
.
responseMap = 'GaFrac = GaFracExpt';
Use an estimatedInfo
object to define
the model parameter kGd
as a parameter to be estimated.
estimatedParam = estimatedInfo('kGd');
Perform the parameter estimation.
fitResult = sbiofit(m1,grpData,responseMap,estimatedParam);
View the estimated parameter value of kGd
.
fitResult.ParameterEstimates
ans = Name Estimate StandardError _____ ________ _____________ 'kGd' 0.11 3.683e05
Suppose you want to simulate the fitted model using different
output times than those in the training data. You can use the predict
method
to do so.
Create a new variable T
with different
output times.
T = [0;10;50;80;100;150;300;350;400;450;500;550];
Use the predict
method to simulate
the fitted model on the new time points. No dosing was specified when
you first ran sbiofit
. Hence, you cannot use any
dosing information with the predict
method, and
an empty array must be specified as the third input argument.
ynew = predict(fitResult,T,[]);
Plot the simulated data with the new output times.
sbioplot(ynew)
Expand Run1 and only select GaFrac to display its simulated data.
This example shows how to estimate categoryspecific (such as young versus old, male versus female) PK parameters using the profile data from multiple individuals using a twocompartment model. The parameters to estimate are the volumes of central and peripheral compartment, the clearance, and intercompartmental clearance.
The synthetic data used in this example contains the time
course of plasma concentrations of multiple individuals after a bolus
dose (100 mg) measured at different times for both central and peripheral
compartments. It also contains categorical variables, namely Sex
and Age
.
clear load(fullfile(matlabroot,'examples','simbio','sd5_302RAgeSex.mat'));
Convert the data set to a groupedData
object,
which is the required data format for the fitting function sbiofit
.
A groupedData
object allows you to set independent
variable and group variable names (if they exist). Set the units of
the ID
, Time
, CentralConc
, PeripheralConc
, Age
,
and Sex
variables. The units are optional and only
required for the UnitConversion
feature, which
automatically converts matching physical quantities to one consistent
unit system.
gData = groupedData(data); gData.Properties.VariableUnits = {'','hour','milligram/liter','milligram/liter','',''};
The IndependentVariableName
and GroupVariableName
properties
have been automatically set to the Time
and ID
variables
of the data.
gData.Properties
ans = Description: '' VariableDescriptions: {} VariableUnits: {'' 'hour' 'milligram/liter' 'milligram/liter' '' ''} DimensionNames: {'Row' 'Variable'} UserData: [] RowNames: {} VariableNames: {'ID' 'Time' 'CentralConc' 'PeripheralConc' 'Sex' 'Age'} GroupVariableName: 'ID' IndependentVariableName: 'Time'
For illustration purposes, use the first five individual data for training and the 6th individual data for testing.
trainData = gData([gData.ID < 6],:); testData = gData([gData.ID == 6],:);
Display the response data for each individual in the training set.
sbiotrellis(trainData,'ID','Time',{'CentralConc','PeripheralConc'});
Use the builtin PK library to construct a twocompartment
model with infusion dosing and firstorder elimination where the elimination
rate depends on the clearance and volume of the central compartment.
Use the configset
object to turn on unit conversion.
pkmd = PKModelDesign; pkc1 = addCompartment(pkmd,'Central'); pkc1.DosingType = 'Bolus'; pkc1.EliminationType = 'linearclearance'; pkc1.HasResponseVariable = true; pkc2 = addCompartment(pkmd,'Peripheral'); model = construct(pkmd); configset = getconfigset(model); configset.CompileOptions.UnitConversion = true;
Assume every individual receives a bolus dose of 100 mg at time = 0.
dose = sbiodose('dose','TargetName','Drug_Central'); dose.StartTime = 0; dose.Amount = 100; dose.AmountUnits = 'milligram'; dose.TimeUnits = 'hour';
The data contains measured plasma concentration in the
central and peripheral compartments. Map these variables to the appropriate
model components, which are Drug_Central
and Drug_Peripheral
.
responseMap = {'Drug_Central = CentralConc','Drug_Peripheral = PeripheralConc'};
Specify the volumes of central and peripheral compartments Central
and Peripheral
,
intercompartmental clearance Q12
, and clearance Cl_Central
as
parameters to estimate. The estimatedInfo
object
lets you optionally specify parameter transforms, initial values,
and parameter bounds. Since both Central and Peripheral are constrained
to be positive, specify a logtransform for each parameter.
paramsToEstimate = {'log(Central)', 'log(Peripheral)', 'Q12', 'Cl_Central'}; estimatedParam = estimatedInfo(paramsToEstimate,'InitialValue',[1 1 1 1]);
Use the 'CategoryVariableName'
property
of the estimatedInfo
object to specify which category
to use during fitting. Use 'Sex'
as the group to
fit for the clearance Cl_Central
and Q12
parameters.
Use 'Age'
as the group to fit for the Central
and Peripheral
parameters.
estimatedParam(1).CategoryVariableName = 'Age'; estimatedParam(2).CategoryVariableName = 'Age'; estimatedParam(3).CategoryVariableName = 'Sex'; estimatedParam(4).CategoryVariableName = 'Sex'; categoryFit = sbiofit(model,trainData,responseMap,estimatedParam,dose)
When fitting by category (or group), sbiofit
always
returns one results object, not one for each category level. This
is because both male and female individuals are considered to be part
of the same optimization using the same error model and error parameters,
similarly for the young and old individuals.
categoryFit = OptimResults with properties: ExitFlag: 3 Output: [1x1 struct] GroupName: [] Beta: [8x5 table] ParameterEstimates: [20x6 table] J: [40x8x2 double] COVB: [8x8 double] CovarianceMatrix: [8x8 double] R: [40x2 double] MSE: 0.1240 SSE: 8.9269 Weights: [] EstimatedParameterNames: {'Central' 'Peripheral' 'Q12' 'Cl_Central'} ErrorModelInfo: [1x3 table] EstimationFunction: 'lsqnonlin' ErrorModel: 'constant' ErrorParameters: [1x1 table]
Plot the categoryspecific estimated results. For the Cl_Central
and Q12
parameters,
all males had the same estimates, and similarly for the females. For
the Central
and Peripheral
parameters,
all young individuals had the same estimates, and similarly for the
old individuals.
plot(categoryFit);
As for testing purposes, simulate the responses of the
6th individual who is an old male. Since you have estimated one set
of parameters for the Age
category (young versus
old), and another set for Sex
category (male versus
female), you can simulate the responses of an old male even though
there is no such individual in the training data.
Use the predict
method to simulate the responses.
ynew
contains simulation data and
paramestim
contains parameter estimates used during
simulation.
[ynew,paramestim] = predict(categoryFit,testData,dose);
Plot the simulated responses of the old male.
sbioplot(ynew);
The paramestim
variable contains the
estimated parameters used by the predict
method.
The parameter estimates for corresponding categories were obtained
from the categoryFit.ParameterEstimates
property.
Specifically, Central
and Peripheral
parameter
estimates are obtained from the Old
group, and Q12
and Cl_Central
parameter
estimates are obtained from the Male
group.
paramestim
paramestim = Name Estimate StandardError Group CategoryVariableName CategoryValue ____________ ________ _____________ _____ ____________________ _____________ 'Central' 1.1993 0.0050591 6 'Age' Old 'Peripheral' 0.5627 0.02908 6 'Age' Old 'Q12' 1.3119 0.062733 6 'Sex' Male 'Cl_Central' 0.56631 0.0079784 6 'Sex' Male
Overlay the experimental results on the simulated data.
figure; plot(testData.Time,testData.CentralConc,'LineStyle','none','Marker','d','MarkerEdgeColor','b'); hold on plot(testData.Time,testData.PeripheralConc,'LineStyle','none','Marker','d','MarkerEdgeColor','r'); plot(ynew.Time,ynew.Data(:,1),'b'); plot(ynew.Time,ynew.Data(:,2),'r'); hold off legend({'OBS1(CentralConc)','OBS2(PeripheralConc)',... 'PRED1(Central.Drug\_Central)','PRED2(Peripheral.Drug\_Peripheral)'});
[1] Yi, TM., Kitano, H., and Simon, M. (2003). A quantitative characterization of the yeast heterotrimeric G protein cycle. PNAS. 100, 10764–10769.
You clicked a link that corresponds to this MATLAB command:
Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
Select web siteYou can also select a web site from the following list:
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.