Этот пример показывает аппроксимацию данных профиля PK индивидуума к модели с одним отделением и оценку фармакокинетических параметров.
Предположим, что у вас есть данные о концентрации лекарственного средства в плазме от индивидуума и вы хотите оценить объем центрального отделения и зазор. Предположим, что концентрация препарата в зависимости от временного профиля следует за моноэкспоненциальным снижением , где - концентрация препарата в момент t, является начальной концентрацией, и - константа скорости устранения, которая зависит от зазора и объема центрального отсека .
Синтетические данные в этом примере были сгенерированы с использованием следующей модели, параметров и дозы:
Однокамерная модель с болюсным дозированием и устранением первого порядка
Объем центрального отсека (Central) = 1,70 л
Параметр зазора (Cl_Central) = 0,55 л/час
Модель постоянной ошибки
Болусовая доза 10 мг
Данные хранятся как таблица с переменными Time и Conc которые представляют временное течение концентрации индивидуума в плазме после внутривенного болюсного введения, измеренное в 13 различных временных точках. Переменные модули для Time и Conc час и миллиграмм/литр, соответственно.
clear all load('data15.mat')
plot(data.Time,data.Conc,'b+') xlabel('Time (hour)'); ylabel('Drug Concentration (milligram/liter)');
Преобразуйте набор данных в groupedData объект, который является необходимым форматом данных для функции аппроксимации sbiofit для дальнейшего использования. A groupedData Объект также позволяет задать независимые имена переменных и групп (если они существуют). Установите модули Time и Conc переменные. Модули являются необязательными и требуются только для UnitConversion функция, которая автоматически преобразует совпадающие физические величины в одну последовательную единичную систему.
gData = groupedData(data);
gData.Properties.VariableUnits = {'hour','milligram/liter'};
gData.Propertiesans = struct with fields:
Description: ''
UserData: []
DimensionNames: {'Row' 'Variables'}
VariableNames: {'Time' 'Conc'}
VariableDescriptions: {}
VariableUnits: {'hour' 'milligram/liter'}
VariableContinuity: []
RowNames: {}
CustomProperties: [1x1 matlab.tabular.CustomProperties]
GroupVariableName: ''
IndependentVariableName: 'Time'
groupedData автоматическая установка имени IndependentVariableName свойство для Time переменная данных.
Используйте встроенную библиотеку PK для создания однокамерной модели с болюсным дозированием и устранением первого порядка, где скорость устранения зависит от зазора и объема центрального отделения. Используйте configset объект для включения модуля.
pkmd = PKModelDesign; pkc1 = addCompartment(pkmd,'Central'); pkc1.DosingType = 'Bolus'; pkc1.EliminationType = 'linear-clearance'; pkc1.HasResponseVariable = true; model = construct(pkmd); configset = getconfigset(model); configset.CompileOptions.UnitConversion = true;
Для получения дополнительной информации о создании компартментных моделей PK с помощью встроенной библиотеки SimBiology ®, смотрите, Создают Фармакокинетические модели.
Задайте одну дозу болюса 10 миллиграммов, заданную во время = 0. Для получения дополнительной информации о настройке различных графиков дозирования, смотрите Дозы в SimBiology Моделях.
dose = sbiodose('dose'); dose.TargetName = 'Drug_Central'; dose.StartTime = 0; dose.Amount = 10; dose.AmountUnits = 'milligram'; dose.TimeUnits = 'hour';
Данные содержат данные о концентрации препарата, хранящиеся в Conc переменная. Эти данные соответствуют Drug_Central виды в модели. Поэтому сопоставьте данные с Drug_Central следующим образом.
responseMap = {'Drug_Central = Conc'};Параметрами, подходящими в этой модели, являются объем центрального отсека (Central) и скорость зазора (Cl_Central). В этом случае задайте логарифмическое преобразование для этих биологических параметров, поскольку они ограничены, чтобы быть положительными. The estimatedInfo Объект позволяет при необходимости задать преобразования параметров, начальные значения и ограничения параметров.
paramsToEstimate = {'log(Central)','log(Cl_Central)'};
estimatedParams = estimatedInfo(paramsToEstimate,'InitialValue',[1 1],'Bounds',[1 5;0.5 2]);Теперь, когда вы определили модель с одним отделением, данные для подгонки, сопоставили данные отклика, параметры для оценки и дозирования, используйте sbiofit для оценки параметров. Функция оценки по умолчанию, которая sbiofit использование будет изменяться в зависимости от того, какие тулбоксы доступны. Чтобы увидеть, какая функция использовалась во время подбора кривой, проверьте EstimationFunction свойство соответствующего объекта результатов.
fitConst = sbiofit(model,gData,responseMap,estimatedParams,dose);
Обратите внимание, что оценки параметров были недалеко от истинных значений (1,70 и 0,55), которые использовались для генерации данных. Можно также попробовать различные модели ошибок, чтобы увидеть, могут ли они еще больше улучшить оценки параметров.
fitConst.ParameterEstimates
ans=2×4 table
Name Estimate StandardError Bounds
______________ ________ _____________ __________
{'Central' } 1.6993 0.034821 1 5
{'Cl_Central'} 0.53358 0.01968 0.5 2
s.Labels.XLabel = 'Time (hour)'; s.Labels.YLabel = 'Concentration (milligram/liter)'; plot(fitConst,'AxesStyle',s);
Попробуйте три другие поддерживаемые модели ошибок (пропорциональные, комбинация постоянных и пропорциональных моделей ошибок и экспоненциальные).
fitProp = sbiofit(model,gData,responseMap,estimatedParams,dose,... 'ErrorModel','proportional'); fitExp = sbiofit(model,gData,responseMap,estimatedParams,dose,... 'ErrorModel','exponential'); fitComb = sbiofit(model,gData,responseMap,estimatedParams,dose,... 'ErrorModel','combined');
Можно задать веса как числовую матрицу, где количество столбцов соответствует количеству откликов. Установка всех весов равной 1 эквивалентна модели постоянной ошибки.
weightsNumeric = ones(size(gData.Conc));
fitWeightsNumeric = sbiofit(model,gData,responseMap,estimatedParams,dose,'Weights',weightsNumeric);Также можно использовать указатель на функцию, который принимает вектор предсказанных значений отклика и возвращает вектор весов. В этом примере используйте указатель на функцию, который эквивалентен пропорциональной модели ошибки.
weightsFunction = @(y) 1./y.^2;
fitWeightsFunction = sbiofit(model,gData,responseMap,estimatedParams,dose,'Weights',weightsFunction);Сравните значения логарифмической правдоподобности, AIC и BIC каждой модели, чтобы увидеть, какая модель ошибки лучше всего подходит для данных. Большее значение вероятности указывает, что соответствующая модель лучше подходит модели. Для AIC и BIC меньшие значения лучше.
allResults = [fitConst,fitWeightsNumeric,fitWeightsFunction,fitProp,fitExp,fitComb];
errorModelNames = {'constant error model','equal weights','proportional weights', ...
'proportional error model','exponential error model',...
'combined error model'};
LogLikelihood = [allResults.LogLikelihood]';
AIC = [allResults.AIC]';
BIC = [allResults.BIC]';
t = table(LogLikelihood,AIC,BIC);
t.Properties.RowNames = errorModelNames;
tt=6×3 table
LogLikelihood AIC BIC
_____________ _______ _______
constant error model 3.9866 -3.9732 -2.8433
equal weights 3.9866 -3.9732 -2.8433
proportional weights -3.8472 11.694 12.824
proportional error model -3.8257 11.651 12.781
exponential error model 1.1984 1.6032 2.7331
combined error model 3.9163 -3.8326 -2.7027
Основываясь на информационных критериях, постоянная модель ошибки (или равные веса) лучше всего подходит для данных, поскольку она имеет самое большое значение логарифмической правдоподобности и самые маленькие AIC и BIC.
Покажите оцененные значения параметров каждой модели.
Estimated_Central = zeros(6,1); Estimated_Cl_Central = zeros(6,1); t2 = table(Estimated_Central,Estimated_Cl_Central); t2.Properties.RowNames = errorModelNames; for i = 1:height(t2) t2{i,1} = allResults(i).ParameterEstimates.Estimate(1); t2{i,2} = allResults(i).ParameterEstimates.Estimate(2); end t2
t2=6×2 table
Estimated_Central Estimated_Cl_Central
_________________ ____________________
constant error model 1.6993 0.53358
equal weights 1.6993 0.53358
proportional weights 1.9045 0.51734
proportional error model 1.8777 0.51147
exponential error model 1.7872 0.51701
combined error model 1.7008 0.53271
Этот пример показал, как оценить параметры PK, а именно объем центрального отсека и параметр зазора индивидуума, путем подгонки данных профиля PK к модели с одним отсеком. Вы сравнили информационные критерии каждой модели и оцененные значения параметров различных моделей ошибок, чтобы увидеть, какая модель лучше всего объяснила данные. Окончательные подгоняемые результаты позволили предположить, что как постоянная, так и комбинированная модели ошибок предоставили самые близкие оценки значениям параметров, используемым для генерации данных. Однако модель постоянной ошибки является лучшей моделью, на которую указывает логарифмическую правдоподобность, AIC и информационные критерии BIC.