Фон

Этот пример показывает, как соответствовать данным о профиле 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(fullfile(matlabroot,'examples','simbio','data15.mat'))

plot(data.Time,data.Conc,'b+')
xlabel('Time (hour)');
ylabel('Drug Concentration (milligram/liter)');

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

Преобразуйте набор данных в объект groupedData, который является необходимым форматом данных для подходящего функционального sbiofit для дальнейшего использования. Объект 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. Для получения дополнительной информации при подготовке различных расписаний дозирования, смотрите Дозы.

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'};

Задайте параметры, чтобы оценить

Параметры, чтобы поместиться в эту модель являются объемом центрального (Центрального) отсека и уровень раскрываемости преступлений (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);

Сравните информационные критерии образцового выбора

Сравните loglikelihood, 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

На основе информационных критериев постоянная ошибочная модель (или равные веса) соответствует данным лучше всего, поскольку это имеет самое большое loglikelihood значение и самый маленький 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 к модели с одним отсеком. Вы сравнили информационные критерии каждой модели и оценили, что значения параметров различных ошибочных моделей видели, какая модель лучше всего объяснила данные. Финал соответствовал результатам, предложенным обоих, которых постоянные и объединенные ошибочные модели предоставили, самые близкие оценки к значениям параметров раньше генерировали данные. Однако постоянная ошибочная модель является лучшей моделью, как обозначено loglikelihood, AIC и критериями информации о BIC.

Фон

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

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

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

clear
load(fullfile(matlabroot,'examples','simbio','sd5_302RAgeSex.mat'))

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

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

Настройте модель 2D отсека

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

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, каждый заканчивается объект для каждого человека.

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

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

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

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' в качестве группы к пригодному для разрешения параметры Q12 и Cl_Central. Используйте 'Age' в качестве группы к пригодному для параметров Peripheral и Central.

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];

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