Фон

Этот пример показывает аппроксимацию данных профиля PK индивидуума к модели с одним отделением и оценку фармакокинетических параметров.

Предположим, что у вас есть данные о концентрации лекарственного средства в плазме от индивидуума и вы хотите оценить объем центрального отделения и зазор. Предположим, что концентрация препарата в зависимости от временного профиля следует за моноэкспоненциальным снижением $$C_t=C_0{e}^{-k_{e}t}$$, где $$C_t$$ - концентрация препарата в момент t, $$C_0$$ является начальной концентрацией, и $$k_e$$ - константа скорости устранения, которая зависит от зазора и объема центрального отсека $$k_e=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 объект, который является необходимым форматом данных для функции аппроксимации 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 переменная данных.

Создайте модель с одним отсеком

Используйте встроенную библиотеку 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;
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       

Заключение

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

Фон

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

Загрузка данных

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

clear
load('sd5_302RAgeSex.mat')

Преобразование в формат сгруппированных данных

Преобразуйте набор данных в 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'

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

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

Отображение данных отклика для каждого индивидуума.

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

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

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

Для получения дополнительной информации о создании компартментных моделей PK с помощью встроенной библиотеки 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 как параметры для оценки. The 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 требует Statistics and Machine Learning 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];

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