Borrar filtros
Borrar filtros

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

9 visualizaciones (últimos 30 días)
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

Respuesta aceptada

Jeremy Huard
Jeremy Huard el 22 de Jul. de 2022
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
  3 comentarios
Jeremy Huard
Jeremy Huard el 2 de Ag. de 2022
Editada: Jeremy Huard el 2 de Ag. de 2022
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.

Iniciar sesión para comentar.

Más respuestas (0)

Comunidades de usuarios

Más respuestas en  SimBiology Community

Categorías

Más información sobre Nonlinear Regression en Help Center y File Exchange.

Productos


Versión

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by