Моделирование фармакокинетики генеральной совокупности фенобарбитала в новорожденных

Этот пример использование:

Этот пример показывает, как создать простую нелинейную модель смешанных эффектов от клинических фармакокинетических данных.

Данные были собраны по 59 недоношенным детям, данным фенобарбитал для предотвращения занятости в течение первых 16 дней после рождения. Каждый человек получил начальную дозу, сопровождаемую одной или несколькими дозами поддержки внутривенным администрированием шарика. В общей сложности между 1 и 6 концентрациями измерения получались от каждого человека время от времени кроме времен дозы, как часть стандартного контроля, для в общей сложности 155 измерений. Младенческие веса и очки APGAR (мера новорожденного здоровья) были также зарегистрированы.

Этот пример использует данные, описанные в [1], исследование, финансируемое NIH/NIBIB, предоставляют P41-EB01975.

Этот пример требует Statistics and Machine Learning Toolbox™.

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

Эти данные были загружены с веб-сайта http://depts.washington.edu/rfpk/ (более не активный) из Средства Ресурса для Фармакокинетики Генеральной совокупности как текстовый файл и сохраненный как массив набора данных для простоты использования.

load pheno.mat ds

Визуализируйте данные в графике решетки

t = sbiotrellis(ds, 'ID', 'TIME', 'CONC', 'marker', 'o',...
       'markerfacecolor', [.7 .7 .7], 'markeredgecolor', 'r', ...
       'linestyle', 'none');

% Format the plot.
t.plottitle = 'States versus Time';
t.updateFigureForPrinting();

Опишите данные

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

data = groupedData(ds);
data.Properties

ans = struct with fields:
                Description: ''
                   UserData: []
             DimensionNames: {'Observations'  'Variables'}
              VariableNames: {'ID'  'TIME'  'DOSE'  'WEIGHT'  'APGAR'  'CONC'}
       VariableDescriptions: {}
              VariableUnits: {}
         VariableContinuity: []
                   RowNames: {}
           CustomProperties: [1x1 matlab.tabular.CustomProperties]
          GroupVariableName: 'ID'
    IndependentVariableName: 'TIME'

Задайте модель

Мы будем соответствовать простой модели с одним отсеком, с администрированием дозы шарика и линейным устранением разрешения, к данным.

pkmd = PKModelDesign;
pkmd.addCompartment('Central', 'DosingType', 'Bolus', 'EliminationType', ...
    'Linear-Clearance', 'HasResponseVariable', true);
model = pkmd.construct;

% The model species |Drug_Central| corresponds to the data variable |CONC|.
responseMap = 'Drug_Central = CONC';

Задайте первоначальные оценки для параметров

Параметры помещаются в эту модель, объем центрального отсека (Central) и уровень раскрываемости преступлений (Cl_Central). NLMEFIT вычисляет зафиксированные и случайные эффекты для каждого параметра. Базовый алгоритм принимает, что параметры нормально распределены. Это предположение не может содержать для биологических параметров, которые ограничиваются быть положительными, такими как объем и разрешение. Мы должны задать преобразование для предполагаемых параметров, таких, что преобразованные параметры действительно следуют за нормальным распределением. По умолчанию SimBiology® выбирает, журнал преобразовывают для всех предполагаемых параметров. Поэтому модель:

$l o g (V_{i}) = l o g (ϕ_{V, i}) = θ_{V} + η_{V, i}$

$l o g (C l_{i}) = l o g (ϕ_{C l, i}) = θ_{C l} + η_{C l, i},$

где $θ$ , $e t a$ , и $ϕ$ фиксированные эффекты, случайные эффекты и оцененные значения параметров соответственно, вычисленный для каждой группы $i$ . Мы обеспечиваем некоторые произвольные первоначальные оценки для V и Статья в отсутствие лучших эмпирических данных.

estimatedParams = estimatedInfo({'log(Central)', 'log(Cl_Central)'}, ...
    'InitialValue', [1 1]);

Извлеките информацию о дозировании от данных

Создайте демонстрационную дозу, которая предназначается для разновидностей Drug_Central, и используйте метод createDoses, чтобы сгенерировать дозы для каждого младенца на основе объема дозирования, перечисленного в переменной DOSE.

sampleDose = sbiodose('sample', 'TargetName', 'Drug_Central');
doses = createDoses(data, 'DOSE', '', sampleDose);

Соответствуйте модели

Немного ослабьте допуски, чтобы ускорить подгонку.

fitOptions.Options = statset('TolFun', 1e-3, 'TolX', 1e-3);
[nlmeResults, simI, simP] = sbiofitmixed(model, data, responseMap, ...
    estimatedParams, doses, 'nlmefit', fitOptions);

Визуализируйте подобранную модель с данными

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

t.plot(simP, [], '', 'Drug_Central');
t.plot(simI, [], '', 'Drug_Central',...
    'legend',{'Observed', 'Fit-Pop.', 'Fit-Indiv.'});

Исследуйте подходящие параметры и ковариации

disp('Summary of initial results');

Summary of initial results

disp('Parameter Estimates Without Random Effects:');

Parameter Estimates Without Random Effects:

disp(nlmeResults.PopulationParameterEstimates(1:2,:));

    Group        Name        Estimate
    _____    ____________    ________

      1      'Central'         1.4086
      1      'Cl_Central'    0.006137

disp('Estimated Fixed Effects:');

Estimated Fixed Effects:

disp(nlmeResults.FixedEffects);

      Name      Description     Estimate    StandardError
    ________    ____________    ________    _____________

    'theta1'    'Central'       0.34257       0.061246   
    'theta2'    'Cl_Central'    -5.0934       0.079501

disp('Estimated Covariance Matrix of Random Effects:');

Estimated Covariance Matrix of Random Effects:

disp(nlmeResults.RandomEffectCovarianceMatrix);

             eta1       eta2  
            _______    _______

    eta1    0.20323          0
    eta2          0    0.19338

Сгенерируйте диаграмму предполагаемых параметров

Этот пример использует команды графического вывода MATLAB®, чтобы визуализировать результаты. Несколько обычно используемых графиков также доступны как встроенные опции при выполнении подгонок параметра через интерфейс рабочего стола SimBiology®.

% Create a box plot of the calculated random effects.
boxplot(nlmeResults);

Сгенерируйте график невязок в зависимости от времени

% The vector of observation data also includes NaN's at the time points at
% which doses were given. We need to remove these NaN's to be able to plot
% the residuals over time.
timeVec = data.(data.Properties.IndependentVariableName);
obsData = data.CONC;
indicesToKeep = ~isnan(obsData);
timeVec = timeVec(indicesToKeep);

% Get the residuals from the fitting results.
indRes = nlmeResults.stats.ires(indicesToKeep);
popRes = nlmeResults.stats.pres(indicesToKeep);

% Plot the residuals. Get a handle to the plot to be able to add more data
% to it later.
resplot = figure;
plot(timeVec,indRes,'d','MarkerFaceColor','b','markerEdgeColor','b');
hold on;
plot(timeVec,popRes,'d','MarkerFaceColor','w','markerEdgeColor','b');
hold off;

% Create a reference line representing a zero residual, and set its
% properties to omit this line from the plot legend.
refl = refline(0,0);
refl.Annotation.LegendInformation.IconDisplayStyle = 'off';

% Label the plot.
title('Residuals versus Time');
xlabel('Time');
ylabel('Residuals');
legend({'Individual','Population'});

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

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

covariateLabels = {'WEIGHT', 'APGAR'};
covariates = grpstats(ds, data.Properties.GroupVariableName, 'mean',...
    'DataVars', covariateLabels);

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

% Get the parameter values for each group (empirical Bayes estimates). 
paramValues = nlmeResults.IndividualParameterEstimates.Estimate;
isCentral = strcmp('Central', nlmeResults.IndividualParameterEstimates.Name);
isCl = strcmp('Cl_Central', nlmeResults.IndividualParameterEstimates.Name);

% Plot the parameter values vs. covariates for each group.
figure;
subplot(2,2,1);
plot(covariates.mean_WEIGHT,paramValues(isCentral), '.');
ylabel('Volume');

subplot(2,2,3);
plot(covariates.mean_WEIGHT,paramValues(isCl), '.');
ylabel('Clearance');
xlabel('weight');

subplot(2,2,2);
plot(covariates.mean_APGAR, paramValues(isCentral), '.');

subplot(2,2,4);
plot(covariates.mean_APGAR, paramValues(isCl), '.');
xlabel('APGAR');

Создайте CovariateModel, чтобы смоделировать ковариационные зависимости

На основе графиков, кажется, существует корреляция между объемом и весом, разрешением и весом, и возможно счетом APGAR и объемом. Мы принимаем решение фокусироваться на эффекте веса путем моделирования двух из этих ковариационных зависимостей: объем ("Центральное") меняться в зависимости от веса и разрешения ("Cl_Central"), меняющийся в зависимости от веса. Наша модель:

$l o g (V_{i}) = l o g (ϕ_{V, i}) = θ_{V} + θ_{V / w e i g h t} * w e i g h t_{i} + η_{V, i}$

$l o g (C l_{i}) = l o g (ϕ_{C l, i}) = θ_{C l} + θ_{C l / w e i g h t} * w e i g h t_{i} + η_{C l, i}$

% Define the covariate model.
covmodel            = CovariateModel;
covmodel.Expression = ({'Central = exp(theta1 + theta2*WEIGHT + eta1)','Cl_Central = exp(theta3 + theta4*WEIGHT + eta2)'});

% Use constructDefaultInitialEstimate to create a initialEstimates struct.
initialEstimates = covmodel.constructDefaultFixedEffectValues;
disp('Fixed Effects Description:');

Fixed Effects Description:

disp(covmodel.FixedEffectDescription);

    'Central'
    'Cl_Central'
    'Central/WEIGHT'
    'Cl_Central/WEIGHT'

Обновите первоначальные оценки с помощью значений, оцененных от подбора кривой базовой модели.

initialEstimates.theta1 = nlmeResults.FixedEffects.Estimate(1);
initialEstimates.theta3 = nlmeResults.FixedEffects.Estimate(2);
covmodel.FixedEffectValues = initialEstimates;

Соответствуйте модели

[nlmeResults_cov, simI_cov, simP_cov] = sbiofitmixed(model, data, responseMap, ...
    covmodel, doses, 'nlmefit', fitOptions);

Исследуйте подходящие параметры и ковариации

disp('Summary of results when modeling covariate dependencies');

Summary of results when modeling covariate dependencies

disp('Parameter Estimates Without Random Effects:');

Parameter Estimates Without Random Effects:

disp(nlmeResults_cov.PopulationParameterEstimates);

    Group        Name        Estimate 
    _____    ____________    _________

     1       'Central'          1.2973
     1       'Cl_Central'    0.0061576
     2       'Central'          1.3682
     2       'Cl_Central'    0.0065512
     3       'Central'          1.3682
     3       'Cl_Central'    0.0065512
     4       'Central'         0.99431
     4       'Cl_Central'    0.0045173
     5       'Central'          1.2973
     5       'Cl_Central'    0.0061576
     6       'Central'          1.1664
     6       'Cl_Central'      0.00544
     7       'Central'          1.0486
     7       'Cl_Central'     0.004806
     8       'Central'          1.1664
     8       'Cl_Central'      0.00544
     9       'Central'          1.2973
     9       'Cl_Central'    0.0061576
     10      'Central'          1.2973
     10      'Cl_Central'    0.0061576
     11      'Central'          1.1664
     11      'Cl_Central'      0.00544
     12      'Central'          1.2301
     12      'Cl_Central'    0.0057877
     13      'Central'          1.1059
     13      'Cl_Central'    0.0051132
     14      'Central'          1.1059
     14      'Cl_Central'    0.0051132
     15      'Central'          1.2301
     15      'Cl_Central'    0.0057877
     16      'Central'          1.1664
     16      'Cl_Central'      0.00544
     17      'Central'          1.1059
     17      'Cl_Central'    0.0051132
     18      'Central'          1.0486
     18      'Cl_Central'     0.004806
     19      'Central'          1.0486
     19      'Cl_Central'     0.004806
     20      'Central'          1.1664
     20      'Cl_Central'      0.00544
     21      'Central'           1.605
     21      'Cl_Central'    0.0078894
     22      'Central'          1.3682
     22      'Cl_Central'    0.0065512
     23      'Central'          3.2052
     23      'Cl_Central'     0.017654
     24      'Central'          3.3803
     24      'Cl_Central'     0.018782
     25      'Central'         0.89394
     25      'Cl_Central'    0.0039908
     26      'Central'          3.9653
     26      'Cl_Central'     0.022619
     27      'Central'          1.6927
     27      'Cl_Central'    0.0083936
     28      'Central'          3.3803
     28      'Cl_Central'     0.018782
     29      'Central'          1.0486
     29      'Cl_Central'     0.004806
     30      'Central'           1.605
     30      'Cl_Central'    0.0078894
     31      'Central'          1.2973
     31      'Cl_Central'    0.0061576
     32      'Central'          4.1819
     32      'Cl_Central'     0.024064
     33      'Central'          1.5218
     33      'Cl_Central'    0.0074154
     34      'Central'          1.5218
     34      'Cl_Central'    0.0074154
     35      'Central'          2.3292
     35      'Cl_Central'     0.012173
     36      'Central'          1.3682
     36      'Cl_Central'    0.0065512
     37      'Central'          1.1664
     37      'Cl_Central'      0.00544
     38      'Central'          1.2301
     38      'Cl_Central'    0.0057877
     39      'Central'          1.6927
     39      'Cl_Central'    0.0083936
     40      'Central'          1.1059
     40      'Cl_Central'    0.0051132
     41      'Central'          1.5218
     41      'Cl_Central'    0.0074154
     42      'Central'          2.7323
     42      'Cl_Central'     0.014659
     43      'Central'         0.99431
     43      'Cl_Central'    0.0045173
     44      'Central'          1.2973
     44      'Cl_Central'    0.0061576
     45      'Central'         0.94279
     45      'Cl_Central'    0.0042459
     46      'Central'          1.1059
     46      'Cl_Central'    0.0051132
     47      'Central'          2.4565
     47      'Cl_Central'     0.012951
     48      'Central'         0.89394
     48      'Cl_Central'    0.0039908
     49      'Central'          1.2301
     49      'Cl_Central'    0.0057877
     50      'Central'          1.1059
     50      'Cl_Central'    0.0051132
     51      'Central'         0.99431
     51      'Cl_Central'    0.0045173
     52      'Central'         0.99431
     52      'Cl_Central'    0.0045173
     53      'Central'          1.5218
     53      'Cl_Central'    0.0074154
     54      'Central'           1.605
     54      'Cl_Central'    0.0078894
     55      'Central'          1.1059
     55      'Cl_Central'    0.0051132
     56      'Central'         0.84763
     56      'Cl_Central'     0.003751
     57      'Central'          1.8827
     57      'Cl_Central'    0.0095009
     58      'Central'          1.2973
     58      'Cl_Central'    0.0061576
     59      'Central'          1.1059
     59      'Cl_Central'    0.0051132

disp('Estimated Fixed Effects:');

Estimated Fixed Effects:

disp(nlmeResults_cov.FixedEffects);

      Name          Description        Estimate    StandardError
    ________    ___________________    ________    _____________

    'theta1'    'Central'              -0.48453      0.067952   
    'theta3'    'Cl_Central'            -5.9575       0.12199   
    'theta2'    'Central/WEIGHT'        0.53203      0.040788   
    'theta4'    'Cl_Central/WEIGHT'     0.61957      0.074264

disp('Estimated Covariance Matrix:');

Estimated Covariance Matrix:

disp(nlmeResults_cov.RandomEffectCovarianceMatrix);

              eta1       eta2  
            ________    _______

    eta1    0.029793          0
    eta2           0    0.04644

Визуализируйте подобранную модель с данными

t.plot(simP_cov, [], '', 'Drug_Central');
t.plot(simI_cov, [], '', 'Drug_Central',...
    'legend',{'Observed', 'Fit-Pop..', 'Fit-Indiv.', 'Cov. Fit-Pop.', 'Cov. Fit-Indiv.'});

Сравните невязки с теми из модели без ковариационных зависимостей

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

% Get the residuals from the fitting results.
indRes = nlmeResults_cov.stats.ires(indicesToKeep);
popRes = nlmeResults_cov.stats.pres(indicesToKeep);

% Return to the original residual plot figure and add the new data.
figure(resplot);
hold on;
plot(timeVec,indRes,'d','MarkerFaceColor','r','markerEdgeColor','r');
plot(timeVec,popRes,'d','MarkerFaceColor','w','markerEdgeColor','r');
hold off;

% Update the legend.
legend('off');
legend({'Individual','Population','Individual(Cov.)','Population(Cov.)'});

Ссылки

[1] Грэзела Т младший, Донн СМ. Неонатальная фармакокинетика генеральной совокупности фенобарбитала выведена от стандартных клинических данных. Фармакол Dev Там 1985:8 (6). 374-83. Краткий обзор PubMed

Документация