How do I fit PK models to multiple dose datasets using simbiology, specifically using the command line (sbiofit)?

I can fit to individual dose data using pooled fiting or mixed effects no problem. However, for some compounds I have multiple doses and I wish to fit to these simultaneously to obtain a single parameter set.

%% PREPARE MODEL
gData = groupedData(data);
gData.Properties.IndependentVariableName = 'Var1';
gData.Properties.GroupVariableName = 'Var3';
gData.Properties.VariableUnits = {'hour','nanogram/milliliter','',''};
gData.Properties;

% One-Compartment Extravascular
pkmd                    = PKModelDesign;
pkc1                    = addCompartment(pkmd,'Central');
pkc1.DosingType         = 'FirstOrder';
pkc1.EliminationType    = 'linear-clearance';
pkc1.HasResponseVariable = true;
model                   = construct(pkmd);
configset               = getconfigset(model);
configset.CompileOptions.UnitConversion = true;

% Single dose
dose                = sbiodose('dose');
dose.TargetName     = 'Dose_Central';
dose.StartTime      = 0;
dose.Amount         = 75;
dose.AmountUnits    = 'milligram';
dose.TimeUnits      = 'hour';
responseMap = {'Drug_Central = Var2'};

% Use NCA parameter estimates as initial parameter guess
NCA.Central = mean(ncaparameters.V_z,'omitnan');
NCA.CL = mean(ncaparameters.CL,'omitnan');
NCA.ka = 1;

paramsToEstimate    = {'log(Central)','log(Cl_Central)','log(ka_Central)'};
% estimatedParams     = estimatedInfo(paramsToEstimate,'InitialValue',[5e-4 1e-6 1],'Bounds',[1 5;0.5 2; 1e-3 10]);
estimatedParams     = estimatedInfo(paramsToEstimate,'InitialValue',[NCA.Central NCA.CL NCA.ka]);


% Fit [Individual, Pooled, Mixed Effects]
fitConst = sbiofit(model,gData,responseMap,estimatedParams,dose);
pooledFit = sbiofit(model,gData,responseMap,estimatedParams,dose,'Pooled',true);

Here is my code for 1 dose. The multiple dose data is in a table similar to the tutorial for phenobarbital, where the next dose and corresponding data follows the previous dose.

I wish to do this using PK models and the sbiofit or similar commands that utilise simbiology. Please can someone shed some light on this?

Cheers

Answer 1

[Jeremy Huard]

Hi Ross,

if the dosing information is included in your dataset (table), then you will be doing pretty much what is done in the Phenobarbital example, except that you will use sbiofit instead of sbiofitmixed.

The difference with the code snippet you shared lies in how the dosing schedule is defined:

• In your code snippet you define a dose object and specify its scheduling. Then, you pass it to sbiofit, which implies that all groups will be simulated with the same dosing schedule.
• For multiple doses, you can (a) define multiple dosing schedules and provide a column vector of doses to sbiofit (one row per group) if dosing is not defined in the dataset, or (b) extract the dosing information from the dataset with createDoses and pass the resulting array of doses to sbiofit.
• In the Phenobarbital example this is done with:

sampleDose = sbiodose('sample', 'TargetName', 'Drug_Central');

doses = createDoses(data, 'DOSE', '', sampleDose);

• Doses are created with information extracted from the DOSE column in the dataset and sampleDose is used as a dose template to define the dose target.

I can also highly recommend to set up the fit in the Model Analyzer App and then, use the View Program Code functionality to generate the code from this program. Then, you can copy/paste what you need instead of typing everyting from scratch.

Best regards,

Jérémy

[Ross Kelly]

Thanks Jeremy.

Two points:

1) When doing this 'from scratch', I get an error regarding dosing units. If I set the dosing units using the code in the initial comment, it does not like that it is assigned to only one.

Error using SimBiology.fit.internal.validateSimFunctionWithDoseInputs (line 53)

The dosing amount has invalid units.

Error in SimBiology.fit.internal.FitObject/fit (line 206)

obj.SimFunctionDoseInputs = SimBiology.fit.internal.validateSimFunctionWithDoseInputs(obj.SimFunction, obj.Phi0, obj.Dosing,

obj.OutputTimesCell, tfEmptyDose);

Error in sbiofit (line 298)

[varargout{1:nargout}] = fitObject.fit(modelObj, data, responseMap, estimInfo, varargin{:});

Error in SbioModelling_temp (line 75)

fitConst = sbiofit(model,gData,responseMap,estimatedParams,dose)

2) I have managed to fit to multiple doses using the Model Analyzer app, and I have viewed the code.

function args = runprogram(model, cs, data, variantsStruct)

% Initialize arguments.
args.input.model    = model;
args.input.cs       = cs;
args.input.data     = data;
args.input.variants = variantsStruct;

% Run fit.
args = runFit(args);

% -------------------------------------------------------------------------
function args = runFit(args)

% Extract the input arguments.
input    = args.input;
data     = input.data;
model    = input.model;
cs       = input.cs;
variants = input.variants.modelStep;

% Set the active configuration set.
originalConfigset = getconfigset(model, 'active');
setactiveconfigset(model, cs);

% Restore the original configset after the task has completed running.
cleanupConfigset = onCleanup(@() restoreActiveConfigset(model, originalConfigset));

% Classify Data Column Names.
GroupLabel = 'GROUP_NUMBER';
TimeLabel  = 'TIME_POINT_HOURS';
DoseLabel  = 'DOSE_TOTAL';
RateLabel  = '';

% Define a description of the data.
groupedDataObj                                    = groupedData(data);
groupedDataObj.Properties.GroupVariableName       = GroupLabel;
groupedDataObj.Properties.IndependentVariableName = TimeLabel;

% Define objects being estimated and their initial estimates.
estimatedInfoObj = estimatedInfo({'log(Central)', 'log(Cl_Central)', 'log(ka_Central)'});
estimatedInfoObj(1).InitialValue = 0.1;
estimatedInfoObj(2).InitialValue = 0.05;
estimatedInfoObj(3).InitialValue = 1;

% Define description of doses.
dose1                  = sbiodose('DOSE_TOTAL');
dose1.TargetName       = 'Central.Dose_Central';
dose1.LagParameterName = '';

% Define doses.
doseTemplate = dose1;
dosing       = createDoses(groupedDataObj, DoseLabel, RateLabel, doseTemplate);

% Define response information.
responseMap = {'Central.Drug_Central = CONCENTRATION'};

% Define fitting options.
fittingOptions = {'ProgressPlot', true, 'ErrorModel', 'constant'};

% Define Algorithm options.
options                     = optimoptions('lsqnonlin');
options.StepTolerance       = 1e-08;
options.FunctionTolerance   = 1e-08;
options.OptimalityTolerance = 1e-06;
options.MaxIterations       = 400;

% Estimate parameter values.
[results, simdataI] = sbiofit(model, groupedDataObj, responseMap, estimatedInfoObj, dosing, 'lsqnonlin', options, variants, fittingOptions{:}, 'Pooled', true);

% Configure the names of the resulting simulation data.
for i = 1:length(simdataI)
    simdataI(i).Name = 'SimData Individual';
end

% Assign output arguments.
args.output.results  = results;
args.output.simdataI = simdataI;


% -------------------------------------------------------------------------
function restoreActiveConfigset(model, cs)

% Restore active configset.
setactiveconfigset(model, cs);

I can't seem to fathom what input arguments I need to use to run this function?

Ultimately, I have data for hundreds of compounds that I wish to use this functionality for. It is important therefore to automate this process.

[Jeremy Huard]

Hi Ross,

1) Do you mean that you used createDoses to create the dose objects?

If yes, the units should be read from the groupedDataObj, provided units are defined in the dataset. They should be defined in data.Properties.VariableUnits .

Alternatively, you could define AmountUnits and TimeUnits of the template sampleDose:

sampleDose.AmountUnits = "milligram";

sampleDose.TimeUnits = "hour";

Or you could set those properties in the final doses objects:

set(doses,'AmountUnits','milligram','TimeUnits','hour');

But if units for the dose column are not defined in the dataset, I suspect units for the dependent and independent variables might not be defined either.

2) If you use R2021a or newwer, you can right-click on the program as you did to view the program code and select 'Export arguments for program code'. This will export a cell array args to the MATLAB workspace.

You can then run the program code with:

runprogram(args{:})

If you use an older version, you will need to generate the input arguments. But if you plan to run this code on different datasets, I recommend to trim down the runprogram code.

For example, if your fit program does not use any baseline variant, you could have a function such as:

function [results, simdataI] = fitData(data)

% One-Compartment Extravascular

pkmd = PKModelDesign;

pkc1 = addCompartment(pkmd,'Central');

pkc1.DosingType = 'FirstOrder';

pkc1.EliminationType = 'linear-clearance';

pkc1.HasResponseVariable = true;

model = construct(pkmd);

configset = getconfigset(model);

configset.CompileOptions.UnitConversion = true;

% Classify Data Column Names.

GroupLabel = 'GROUP_NUMBER';

TimeLabel = 'TIME_POINT_HOURS';

DoseLabel = 'DOSE_TOTAL';

RateLabel = '';

% Define a description of the data.

groupedDataObj = groupedData(data);

groupedDataObj.Properties.GroupVariableName = GroupLabel;

groupedDataObj.Properties.IndependentVariableName = TimeLabel;

% Define objects being estimated and their initial estimates.

estimatedInfoObj = estimatedInfo({'log(Central)', 'log(Cl_Central)', 'log(ka_Central)'});

estimatedInfoObj(1).InitialValue = 0.1;

estimatedInfoObj(2).InitialValue = 0.05;

estimatedInfoObj(3).InitialValue = 1;

% Define description of doses.

dose1 = sbiodose('DOSE_TOTAL');

dose1.TargetName = 'Central.Dose_Central';

dose1.LagParameterName = '';

% Define doses.

doseTemplate = dose1;

dosing = createDoses(groupedDataObj, DoseLabel, RateLabel, doseTemplate);

% Define response information.

responseMap = {'Central.Drug_Central = CONCENTRATION'};

% Define fitting options.

fittingOptions = {'ProgressPlot', true, 'ErrorModel', 'constant'};

% Define Algorithm options.

options = optimoptions('lsqnonlin');

options.StepTolerance = 1e-08;

options.FunctionTolerance = 1e-08;

options.OptimalityTolerance = 1e-06;

options.MaxIterations = 400;

% Estimate parameter values.

[results, simdataI] = sbiofit(model, groupedDataObj, responseMap, estimatedInfoObj, ...

dosing, 'lsqnonlin', options, [], fittingOptions{:}, 'Pooled', true);

end

Then, you can have a script that loads the datasets (for example with readtable) and run this function.

[Ross Kelly]

Thank you, Jeremy!