Моделирование фармакокинетики населения фенобарбитала в новорожденных

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

Данные были собраны по 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: [1×1 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® выбирает, журнал преобразовывают для всех предполагаемых параметров. Поэтому модель:

log(Vi)=log(ϕV,i)=θV+ηV,i

и

log(Cli)=log(ϕCl,i)=θCl+ηCl,i,

где θ, eta, и ϕ фиксированные эффекты, случайные эффекты и оцененные значения параметров соответственно, вычисленный для каждой группы 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"), меняющийся в зависимости от веса. Наша модель:

log(Vi)=log(ϕV,i)=θV+θV/weight*weighti+ηV,i

и

log(Cli)=log(ϕCl,i)=θCl+θCl/weight*weighti+ηCl,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.

Для просмотра документации необходимо авторизоваться на сайте