Фон

Этот пример показывает, как подогнать данные профиля ПК индивидуума к однокамерной модели и оценить фармакокинетические параметры.

Предположим, что у вас есть данные о концентрации лекарственного средства в плазме от человека и вы хотите оценить объем центрального отделения и клиренс. Предположим, что концентрация лекарственного средства в сравнении с профилем времени следует за моноэкспоненциальным снижением $${}_{\mathrm{Ct}}{=}_{\mathrm{}}{}^{{}_{}}$$C0e-ket, $${\mathrm{где}}_{}$$Ct - концентрация лекарственного средства в момент t$${,}_{\mathrm{}}$$C0 - начальная концентрация, $${и}_{}$$ke - постоянная скорости удаления, которая зависит от зазора и объема центрального $${\mathrm{отделения}}_{}\mathrm{ke}=$$Cl/V.

Синтетические данные в этом примере были получены с использованием следующей модели, параметров и дозы:

Загрузка данных и визуализация

Данные хранятся в виде таблицы с переменными. 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

Преобразовать набор данных в groupedData объект, который является требуемым форматом данных для функции фитинга sbiofit для последующего использования. A groupedData объект также позволяет задавать независимые имена переменных и групп (если они существуют). Установка единиц измерения Time и Conc переменные. Установки являются необязательными и обязательными только для UnitConversion функция, которая автоматически преобразует соответствующие физические величины в одну согласованную систему единиц измерения.

gData = groupedData(data);
gData.Properties.VariableUnits = {'hour','milligram/liter'};
gData.Properties

ans = 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 переменная данных.

Построение однокамерной модели

Используйте встроенную библиотеку ПК для построения однокамерной модели с болюсным дозированием и устранением первого порядка, где скорость устранения зависит от зазора и объема центрального отсека. Используйте 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;

Дополнительные сведения о создании моделей компартментальных ПК с использованием встроенной библиотеки 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). В этом случае укажите логарифмическое преобразование для этих биологических параметров, поскольку они ограничены положительными. 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;
t

t=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       

Заключение

Этот пример показал, как оценить параметры ПК, а именно объем центрального отсека и параметр зазора индивидуума, путем подгонки данных профиля ПК к однокамерной модели. Вы сравнили информационные критерии каждой модели и оценочные значения параметров различных моделей ошибок, чтобы увидеть, какая модель лучше всего объяснила данные. Окончательные подогнанные результаты предлагали как постоянные, так и комбинированные модели ошибок, предоставляющие самые близкие оценки к значениям параметров, используемым для генерации данных. Тем не менее, модель постоянной ошибки является лучшей моделью, о чем свидетельствуют критерии источника информации, AIC и BIC.

Фон

Предположим, что у вас есть данные о концентрации лекарственного средства в плазме от 30 человек и вы хотите оценить фармакокинетические параметры, а именно объемы центрального и периферического отделения, клиренс и межкомпартментальный клиренс. Предположим, что концентрация лекарственного средства по сравнению с профилем времени следует за бикспоненциальным снижением $${}_{\mathrm{Ct}}={}^{}\mathrm{Ae-at}{}^{+}$$Be-bt, $${}_{}$$где Ct - концентрация лекарственного средства в момент t, $$$$а a $$$$и b - наклоны для соответствующих экспоненциальных понижений.

Загрузить данные

Эти синтетические данные содержат временной курс концентраций в плазме 30 особей после болюсной дозы (100 мг), измеренной в разное время как для центрального, так и для периферического отделения. Он также содержит категориальные переменные, а именно: Sex и Age.

clear
load('sd5_302RAgeSex.mat')

Преобразуйте в Формат groupedData

Преобразовать набор данных в groupedData объект, который является требуемым форматом данных для функции фитинга sbiofit. A groupedData объект также позволяет задавать независимые имена переменных и групп (если они существуют). Установка единиц измерения ID, Time, CentralConc, PeripheralConc, Age, и Sex переменные. Установки являются необязательными и обязательными только для UnitConversion функция, которая автоматически преобразует соответствующие физические величины в одну согласованную систему единиц измерения.

gData = groupedData(data);
gData.Properties.VariableUnits = {'','hour','milligram/liter','milligram/liter','',''};
gData.Properties

ans = struct with fields:
                Description: ''
                   UserData: []
             DimensionNames: {'Row'  'Variables'}
              VariableNames: {1x6 cell}
       VariableDescriptions: {}
              VariableUnits: {1x6 cell}
         VariableContinuity: []
                   RowNames: {}
           CustomProperties: [1x1 matlab.tabular.CustomProperties]
          GroupVariableName: 'ID'
    IndependentVariableName: 'Time'

IndependentVariableName и GroupVariableName свойства были автоматически установлены в значение Time и ID переменные данных.

Визуализация данных

Просмотрите данные ответа для каждого отдельного пользователя.

t = sbiotrellis(gData,'ID','Time',{'CentralConc','PeripheralConc'},...
                'Marker','+','LineStyle','none');
% Resize the figure.
t.hFig.Position(:) = [100 100 1280 800];

Настройка двухкамерной модели

Используйте встроенную библиотеку ПК для построения двухкамерной модели с инфузионным дозированием и устранением первого порядка, где скорость устранения зависит от зазора и объема центрального отсека. Используйте configset объект для включения преобразования единиц измерения.

pkmd                                    = PKModelDesign;
pkc1                                    = addCompartment(pkmd,'Central');
pkc1.DosingType                         = 'Bolus';
pkc1.EliminationType                    = 'linear-clearance';
pkc1.HasResponseVariable                = true;
pkc2                                    = addCompartment(pkmd,'Peripheral');
model                                   = construct(pkmd);
configset                               = getconfigset(model);
configset.CompileOptions.UnitConversion = true;

Дополнительные сведения о создании моделей компартментальных ПК с использованием встроенной библиотеки SimBiology ® см. в разделе Создание фармакокинетических моделей.

Определение дозирования

Предположим, что каждый человек получает болюсную дозу 100 мг в момент времени = 0. Для получения подробной информации о настройке различных стратегий дозирования см. Дозы в моделях SimBiology.

dose             = sbiodose('dose','TargetName','Drug_Central');
dose.StartTime   = 0;
dose.Amount      = 100;
dose.AmountUnits = 'milligram';
dose.TimeUnits   = 'hour';

Сопоставление данных ответа с соответствующими компонентами модели

Данные содержат измеренную концентрацию плазмы в центральном и периферийном отсеках. Сопоставьте эти переменные с соответствующими компонентами модели, которые являются Drug_Central и Drug_Peripheral.

responseMap = {'Drug_Central = CentralConc','Drug_Peripheral = PeripheralConc'};

Укажите параметры для оценки

Указать объемы центрального и периферийного отсеков Central и Peripheral, межотраслевой клиренс Q12и клиренс Cl_Central в качестве параметров для оценки. estimatedInfo объект позволяет дополнительно указать преобразования параметров, начальные значения и границы параметров. Поскольку оба Central и Peripheral должны быть положительными, укажите логарифмическое преобразование для каждого параметра.

paramsToEstimate    = {'log(Central)', 'log(Peripheral)', 'Q12', 'Cl_Central'};
estimatedParam      = estimatedInfo(paramsToEstimate,'InitialValue',[1 1 1 1]);

Оценка индивидуальных параметров

Оцените один набор параметров для каждого человека, установив 'Pooled' аргумент пары имя-значение для false.

unpooledFit =  sbiofit(model,gData,responseMap,estimatedParam,dose,'Pooled',false);

Показать результаты

Постройте график соответствия результатов исходным данным для каждого человека (группы).

t = plot(unpooledFit);
% Resize the figure.
t.hFig.Position(:) = [100 100 1280 800];

Для незаполненной посадки, sbiofit всегда возвращает один объект результатов для каждого отдельного объекта.

Анализ оценок параметров для зависимостей категорий

Изучите незаполненные оценки, чтобы выяснить, существуют ли какие-либо специфичные для категории параметры, то есть если некоторые параметры связаны с одной или несколькими категориями. При наличии каких-либо зависимостей категорий можно уменьшить число степеней свободы путем оценки только специфичных для категории значений для этих параметров.

Сначала извлеките значения идентификаторов и категорий для каждого идентификатора.

catParamValues = unique(gData(:,{'ID','Sex','Age'}));

Добавьте переменные в таблицу, содержащую оценку каждого параметра.

allParamValues            = vertcat(unpooledFit.ParameterEstimates);
catParamValues.Central    = allParamValues.Estimate(strcmp(allParamValues.Name, 'Central'));
catParamValues.Peripheral = allParamValues.Estimate(strcmp(allParamValues.Name, 'Peripheral'));
catParamValues.Q12        = allParamValues.Estimate(strcmp(allParamValues.Name, 'Q12'));
catParamValues.Cl_Central = allParamValues.Estimate(strcmp(allParamValues.Name, 'Cl_Central'));

Печать оценок каждого параметра для каждой категории. gscatter требуются Toolbox™ статистики и машинного обучения. Если он отсутствует, используйте другие альтернативные функции печати, такие как plot.

h           = figure;
ylabels     = {'Central','Peripheral','Cl\_Central','Q12'};
plotNumber  = 1;
for i = 1:4
    thisParam = estimatedParam(i).Name;
    
    % Plot for Sex category
    subplot(4,2,plotNumber);
    plotNumber  = plotNumber + 1;
    gscatter(double(catParamValues.Sex), catParamValues.(thisParam), catParamValues.Sex);
    ax          = gca;
    ax.XTick    = [];
    ylabel(ylabels(i));
    legend('Location','bestoutside')
    % Plot for Age category
    subplot(4,2,plotNumber);
    plotNumber  = plotNumber + 1;
    gscatter(double(catParamValues.Age), catParamValues.(thisParam), catParamValues.Age);
    ax          = gca;
    ax.XTick    = [];
    ylabel(ylabels(i));
    legend('Location','bestoutside')
end
% Resize the figure.
h.Position(:) = [100 100 1280 800];

Исходя из сюжета, кажется, что молодые люди, как правило, имеют более высокие объемы центрального и периферического отделений (Central, Peripheral), чем старые люди (то есть объемы, по-видимому, зависят от возраста). Кроме того, мужчины, как правило, имеют более высокие показатели клиренса (Cl_Central), чем женщины, но наоборот для параметра Q12 (то есть клиренс и Q12, по-видимому, зависят от пола).

Параметры, специфичные для категории оценки

Используйте 'CategoryVariableName' имущества estimatedInfo для указания категории, используемой во время фитинга. Использовать 'Sex' как группа, подходящая для зазора Cl_Central и Q12 параметры. Использовать 'Age' как группа, подходящая для Central и Peripheral параметры.

estimatedParam(1).CategoryVariableName = 'Age';
estimatedParam(2).CategoryVariableName = 'Age';
estimatedParam(3).CategoryVariableName = 'Sex';
estimatedParam(4).CategoryVariableName = 'Sex';
categoryFit = sbiofit(model,gData,responseMap,estimatedParam,dose)

categoryFit = 
  OptimResults with properties:

                   ExitFlag: 3
                     Output: [1x1 struct]
                  GroupName: []
                       Beta: [8x5 table]
         ParameterEstimates: [120x6 table]
                          J: [240x8x2 double]
                       COVB: [8x8 double]
           CovarianceMatrix: [8x8 double]
                          R: [240x2 double]
                        MSE: 0.4362
                        SSE: 205.8690
                    Weights: []
              LogLikelihood: -477.9195
                        AIC: 971.8390
                        BIC: 1.0052e+03
                        DFE: 472
             DependentFiles: {1x3 cell}
    EstimatedParameterNames: {'Central'  'Peripheral'  'Q12'  'Cl_Central'}
             ErrorModelInfo: [1x3 table]
         EstimationFunction: 'lsqnonlin'

При подборе по категориям (или группам), sbiofit всегда возвращает один объект результатов, а не один для каждого уровня категории. Это связано с тем, что как мужчины, так и женщины считаются частью одной и той же оптимизации, используя одну и ту же модель ошибок и параметры ошибок, аналогично для молодых и старых людей.

Результаты графика

Постройте график оценочных результатов для конкретной категории.

t = plot(categoryFit);
% Resize the figure.
t.hFig.Position(:) = [100 100 1280 800];

Для Cl_Central и Q12 параметры, все мужчины имели одинаковые оценки, и аналогично для женщин. Для Central и Peripheral параметры, все молодые люди имели одинаковые оценки и аналогично для старых людей.

Оценка параметров для всего населения

Чтобы лучше сравнить результаты, подберите модель для всех данных, объединенных вместе, то есть оцените один набор параметров для всех пользователей, установив 'Pooled' аргумент пары имя-значение для true. Предупреждающее сообщение указывает на то, что эта опция игнорирует любую информацию, относящуюся к конкретной категории (если она существует).

pooledFit = sbiofit(model,gData,responseMap,estimatedParam,dose,'Pooled',true);

Warning: CategoryVariableName property of the estimatedInfo object is ignored when using the 'Pooled' option.

Результаты графика

Постройте график соответствия результатов исходным данным. Хотя для каждого индивидуума был сформирован отдельный график, данные подгонялись с использованием одного и того же набора параметров (то есть все индивидуумы имели одинаковую подгоняемую линию).

t = plot(pooledFit);
% Resize the figure.
t.hFig.Position(:) = [100 100 1280 800];

Сравнить остатки

Сравнить остатки CentralConc и PeripheralConc отклики для каждой посадки.

t = gData.Time;
allResid(:,:,1) = pooledFit.R;
allResid(:,:,2) = categoryFit.R;
allResid(:,:,3) = vertcat(unpooledFit.R);

h = figure;
responseList = {'CentralConc', 'PeripheralConc'};
for i = 1:2
    subplot(2,1,i);
    oneResid = squeeze(allResid(:,i,:));
    plot(t,oneResid,'o');
    refline(0,0); % A reference line representing a zero residual
    title(sprintf('Residuals (%s)', responseList{i}));
    xlabel('Time');
    ylabel('Residuals');
    legend({'Pooled','Category-Specific','Unpooled'});
end
% Resize the figure.
h.Position(:) = [100 100 1280 800];

Как показано на графике, незаполненная посадка обеспечивает наилучшее соответствие данным, поскольку она соответствует данным каждого отдельного пользователя. Это ожидалось, поскольку в нем использовалось наибольшее количество степеней свободы. Подгонка по категориям уменьшила число степеней свободы за счет подгонки данных к двум категориям (пол и возраст). В результате остатки были больше, чем незаполненная посадка, но все же меньше, чем популяция, которая оценивала только один набор параметров для всех индивидуумов. Соответствие категориям может быть хорошим компромиссом между незаполненным и объединенным фитингом при условии, что в данных существует любая иерархическая модель.