Моделирование PK/PD и Симуляция, чтобы Вести Стратегию Дозирования Антибиотиков

Этот пример использует:

В этом примере показано, как выполнить симуляцию Монте-Карло фармакокинетической/фармакодинамической модели (PK/PD) для антибактериального средства. Этот пример адаптируется от Katsube и др. [1], Этот пример также показывает, как использовать SimBiology® SimFunction объект выполнить сканы параметра параллельно.

Этот пример требует Statistics and Machine Learning Toolbox™. Производительность может улучшаться, если у вас есть программное обеспечение Parallel Computing Toolbox™.

Фон

Katsube и др. [1] использовал подход моделирования и симуляции PK/PD, чтобы определить самые эффективные режимы дозы для doripenem, антибиотика карбапенемы. Цели их исследования были:

Разработайте модель PK/PD, чтобы описать антибактериальный эффект doripenem против нескольких деформаций Pseudomonas aeruginosa
Используйте симуляции Монте-Карло, чтобы сравнить эффективность четырех общих антибиотических режимов дозы и определить самую эффективную стратегию дозирования
Исследуйте эффект почечной функции на антибактериальной эффективности обработок

В этом примере мы реализуем антибактериальную модель PK/PD, разработанную Katsube и др. [1] в SimBiology®, и реплицируем результаты симуляции Монте-Карло, описанной в их работе.

Ссылки

[1] Т. Кэтсьюб, И. Яно, T. Вадзима, И. Ямано и М. Такано. Фармакокинетическое/фармакодинамическое моделирование и симуляция, чтобы определить эффективные режимы дозы для doripenem. Журнал Фармацевтических Наук (2010) 99 (5), 2483-91.

Модель PK/PD

Katsube и др. принял модель вливания 2D отсека с линейным устранением из центрального отсека, чтобы описать фармакокинетику doripenem. Для бактериальной модели роста они приняли, что общее бактериальное население состоит из восприимчивых к препарату растущих ячеек и нечувствительных к препарату ячеек отдыха. Антибактериальный эффект препарата был включен в уровень уничтожения бактерий с помощью простой модели типа Emax:

$Killing Rate = \frac{Kmax*[Drug]*[Growing]}{KC50 + [Drug]}$

где [Drug] концентрация (ug/ml) препарата в центральном отсеке и [Growing] количество растущего бактериального населения в CFU/ml (CFU = Модули Формирования Колонии). Kmax максимальная константа скорости уничтожения (1/час) и KC50 константа скорости Михаэлиса-Ментен (ug/ml).

Графическое представление реализации SimBiology модели показывают ниже.

% Load model
sbioloadproject('AntibacterialPKPD.sbproj', 'm1') ;

Режимы дозы

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

250 мг два раза в день (b.i.d.)
250 мг три раза в день (t.i.d.)
500 мг два раза в день (b.i.d.)
500 мг три раза в день (t.i.d.)

Дозирование вливания использовалось во всех четырех режимах доз, и время вливания было установлено в 30 минут. В SimBiology были реализованы эти режимы дозы, когда доза возражает.

% Select dose objects in the model
doseNames = {'250 mg bid', '250 mg tid', '500 mg bid', '500 mg tid'};
for iDoseGrp = 1:length(doseNames)
    doseRegimens(iDoseGrp) = sbioselect(m1, 'Name', doseNames{iDoseGrp}) ;
end

Описание виртуального населения

Виртуальное население индивидуумов было сгенерировано на основе распределения демографических переменных и параметров PK/PD. Тип распределения и значения параметров распределения были основаны на данных из более ранних клинических испытаний doripenem, проводимого в Японии.

Примечание: В [1], 5 000 виртуальных пациентов были симулированы в каждой группе дозы. В этом примере мы будем использовать 1 000 пациентов в каждой группе. Чтобы симулировать различную численность населения, измените значение nPatients ниже.

% Setup
nPatients    = 1000         ; % Number of patients per dosage group
nDoseGrps    = 4            ; % Number of tested dosage regimens

Распределение демографических переменных:

Вес (Wt) и возраст (Age) были произведены от нормального распределения со средним значением 51,6 кг и 71,8 года, соответственно, и стандартное отклонение 11,8 кг и 11,9 лет, соответственно. 26% населения были приняты, чтобы быть розеткой. Уровни креатинина сыворотки (Scr) были произведены от логарифмически нормального распределения с типичным значением (среднегеометрическим) из 0,78 мг/дл и коэффициентом изменчивости (CV) 32,8%. Уровень раскрываемости преступлений креатинина (CrCL) были вычислены с помощью уравнения Cockcroft-Gault.

Входные параметры к lognrnd функция является средним значением (mu) и стандартное отклонение (sigma) из связанного нормального распределения. Здесь и в примере, mu и sigma были вычислены от типичного значения, о котором сообщают, и коэффициента изменчивости логарифмически нормального распределения. Можно использовать следующие определения, чтобы вычислить их. Смотрите lognstat документация для получения дополнительной информации.

mu    = @(m,v) log(m^2/sqrt(v+m^2));
sigma = @(m,v) sqrt(log(v/m^2+1));
m     = @(typicalValue) typicalValue;
v     = @(typicalValue,CV) typicalValue^2*CV^2;

% Patient demographics
rng('default');
Wt        = normrnd(51.6, 11.8, nPatients , nDoseGrps )   ; % units: kg
Age       = normrnd(71.8, 11.9, nPatients , nDoseGrps )   ; % units: years
Scr_mu    =    mu(m(0.78), v(0.78,0.328));
Scr_sigma = sigma(m(0.78), v(0.78,0.328));
Scr       = lognrnd(Scr_mu, Scr_sigma, nPatients , nDoseGrps )   ; % units: ml/minute
% Gender ratio
id        = 1:nPatients*nDoseGrps                            ;
idFemale  = randsample(id , round(0.26*nDoseGrps*nPatients)) ; % 26% Female

Разрешение креатинина (использующий уравнение Cockcroft-Gault)

CrCL            = (140 - Age).*Wt./(Scr*72)                 ; % units: ml/minute
CrCL(idFemale)  = CrCL(idFemale)*0.85                       ; % multiply by 0.85 for females

Распределение Фармакокинетических (PK) параметры:

Параметры PK, Central, k12, и k21, были произведены от логарифмически нормального распределения с типичными значениями 7,64 литров, 1.59 1/час и 2.26 1/час, соответственно, и 20% коэффициента изменчивости (CV). Central объем распределения центрального отсека и k12 и k21 константы скорости передачи между Central и Peripheral отсеки.

Central_mu    =    mu(m(7.64), v(7.64,0.20));
Central_sigma = sigma(m(7.64), v(7.64,0.20));
k12_mu        =    mu(m(1.59), v(1.59,0.20));
k12_sigma     = sigma(m(1.59), v(1.59,0.20));
k21_mu        =    mu(m(2.26), v(2.26, 0.2));
k21_sigma     = sigma(m(2.26), v(2.26, 0.2));

Central = lognrnd(Central_mu , Central_sigma, nPatients , nDoseGrps); % units: liter
k12     = lognrnd(k12_mu, k12_sigma, nPatients , nDoseGrps)         ; % units: 1/hour
k21     = lognrnd(k21_mu, k21_sigma, nPatients , nDoseGrps)         ; % units: 1/hour

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

$CL = 1.07*CrCL + 45.6 + \varepsilon$

где $\varepsilon$ аддитивная остаточная ошибка производится от нормального распределения со средним значением 0 мл/минута и стандартным отклонением 22 мл/минута.

CL      = 1.07*CrCL + 45.6 + normrnd(0,22,nPatients,nDoseGrps); % units: ml/minute

Распределение Фармакодинамических (PD) параметры:

Константы скорости преобразования роста к отдыху, k1 и k2, были произведены от логарифмически нормального распределения с типичным значением 5.59e-5 и 0.0297 1/час, соответственно, каждый с CV 20%. Kmax был произведен от логарифмически нормального распределения с типичным значением 3,5 1/час и CV на 15,9%.

k1_mu      =    mu(m(5.59e-5), v(5.59e-5, 0.2));
k1_sigma   = sigma(m(5.59e-5), v(5.59e-5, 0.2));
k2_mu      =    mu(m(0.0297) , v(0.0297, 0.2));
k2_sigma   = sigma(m(0.0297) , v(0.0297, 0.2));
Kmax_mu    =    mu(m(3.50)   , v(3.50, 0.159));
Kmax_sigma = sigma(m(3.50)   , v(3.50, 0.159));

k1   = lognrnd(k1_mu, k1_sigma, nPatients , nDoseGrps)     ; % units: 1/hour
k2   = lognrnd(k2_mu, k2_sigma, nPatients , nDoseGrps)     ; % units: 1/hour
Kmax = lognrnd(Kmax_mu, Kmax_sigma, nPatients , nDoseGrps) ; % units: 1/hour

Katsube и др. принял это значения k1, k2 и Kmax были независимы от обрабатываемого бактериального штамма. Значение Beta, сетевой постоянный темп роста, был зафиксирован в 1,5 1/час.

На основе экспериментов с несколькими деформациями авторы пришли к заключению что значение KC50 было линейно зависимо на минимальной концентрации ингибирования (MIC) бактериального штамма через следующее уравнение.

$ln(KC50) = -1.91 + 0.898*ln(MIC) + \varepsilon$

где $\varepsilon$ аддитивная остаточная ошибка производится от нормального распределения со средним значением 0 и стандартным отклонением 1.06 ug/ml. В симуляции, MIC значения были произведены от дискретного распределения и KC50 значение было вычислено для выбранного MIC использование вышеупомянутого уравнения.

% Discrete distribution of MIC values based on 71 P. aeruginosa strains
micValue  = [0.0625, 0.125, 0.25, 0.5 , 1 , 2 , 4 , 8 , 16 , 32 ] ;
micFreq   = [   5  ,    8 ,  9  , 14  , 7 , 8 , 9 , 5 , 2  , 4  ] ;

% Sample MIC values from a discrete distribution using randsample
MIC = nan(nPatients, nDoseGrps) ; % preallocate
for iDoseGrp = 1:nDoseGrps
    MIC(:, iDoseGrp) = randsample(micValue , nPatients, true , micFreq);
end

KC50 = exp(-1.91 + 0.898*log(MIC) + 1.06*randn(nPatients , nDoseGrps)) ; % units: microgram/milliliter

Симуляция Setup & Design

Создайте SimFunction объект, который позволяет вам выполнить симуляции модели и сканы параметра параллельно. В этом примере вы будете варьироваться 8 параметров, Central, k12, k21, CL, k1, k2, Kmax, и KC50. Выберите рост и отдых бактериальных количеств, Growing и Resting, как ответы, то есть, результаты симуляции, которые вы хотите наблюдать при варьировании тех входных параметров.

Выберите параметры, чтобы варьироваться.

params      = {'Central', 'k12', 'k21', 'CL', 'k1', 'k2', 'Kmax', 'KC50'};

Выберите ответы.

observables = {'[Bacterial Growth Model].Growing',...
               '[Bacterial Growth Model].Resting'};

Настройте дозу шаблона.

tempdose = sbiodose('dose');
tempdose.Target = 'Central.Drug';
tempdose.AmountUnits = 'milligram';
tempdose.TimeUnits = 'hour';
tempdose.DurationParameterName = 'TDose';

Создайте SimFunction объект. Установите UseParallel к истине, чтобы включить параллельные вычисления.

simfunc  = createSimFunction(m1,params,observables,tempdose,'UseParallel',true);

Создайте входную матрицу phi для каждой группы дозы.

phi = cell(1,nDoseGrps);
for i = 1:nDoseGrps
    phi{i} = [Central(:,i),k12(:,i),k21(:,i), ...
              CL(:,i), k1(:,i), k2(:,i), ...
              Kmax(:,i), KC50(:,i)];
end

Кластерные вычисления

Этот пример использует локальный кластерный профиль, который предварительно сконфигурирован к вашей локальной машине. Можно также искать другие кластеры MATLAB® Parallel Server™, которые работают на Amazon EC2®. На вкладке Home в разделе Environment выберите Parallel> Discover Clusters. Чтобы получить доступ к этим кластерам, необходимо предоставить информацию о входе в систему Учетной записи MathWorks®. Для получения дополнительной информации смотрите, Обнаруживают Кластеры и Профили Кластера Использования (Parallel Computing Toolbox).

Создайте параллельный пул, если ни один не существует.

if isempty(gcp)
    parpool;
end

Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 6).

Для всех сценариев дозы модель была симулирована до t = 2 недели со времени первой дозы. Общее бактериальное количество, CFU, было произведено каждые 24 часа (один раз в день) на целое время режима дозы.

tObs      = 0:24:336       ; % hour
nTPoints  = length(tObs)   ; % Number of sampling points

Симуляция Монте-Карло пациентов с серьезным заражением

Антибактериальная эффективность препарата может быть измерена с помощью различных индексов PK/PD. Katsube и др. устанавливают критерий бактериального устранения в log10(CFU) < 0, где CFU является общим бактериальным количеством. Эффективность каждого режима дозы была измерена как часть населения, которое достигло критериев успеха в группе дозы. Эта метрика эффективности, Pr{log10(CFU) < 0}, был прослежен в зависимости от времени для каждой группы дозы.

В их исследованиях симуляции авторы исследовали эффективность режимов дозы на двух классах пациентов:

Умеренное заражение (Начальное бактериальное количество = 1e4 CFU/ml)
Серьезное заражение (Начальное бактериальное количество = 1e7 CFU/ml)

В этом примере мы реплицируем результаты для серьезного случая заражения только. Обратите внимание на то, что можно легко симулировать другой сценарий, пациентов с умеренным заражением, путем изменения начальной суммы бактериального количества (Growing разновидности), в модели к 1e4 CFU/ml.

% Preallocate
cfu          = nan(nTPoints,nPatients);
log10CFU     = cell(1,nDoseGrps) ;

for i = 1:nDoseGrps
    disp(['Simulating group ', num2str(i),'...'])

    % Get the dose table directly from an existing dose object for each
    % dosing regimen.
    doseTable = getTable(doseRegimens(i));

    % Simulate
    simdata = simfunc(phi{i},[],doseTable,tObs);

    % Sum of growing and resting bacterial counts for each patient
    for j = 1:nPatients
       cfu(:,j) = sum(simdata(j).Data,2);
    end
    % Store log-transformed counts for each dose group.
    log10CFU{i} = log10(cfu);
end

% Save results
log10CFU_250bid = log10CFU{1} ;
log10CFU_250tid = log10CFU{2} ;
log10CFU_500bid = log10CFU{3} ;
log10CFU_500tid = log10CFU{4} ;

Simulating group 1...
Simulating group 2...
Simulating group 3...
Simulating group 4...

Закройте параллельный пул.

delete(gcp('nocreate'));

Профили курса времени бактериальных количеств

Мы строим медиану (в красном), и процентиль (заштриховала) профили log10(CFU) уровни для всех четырех режимов дозы. Заметьте, что во всех четырех группах, средний профиль курса времени показывает, что бактериальное уничтожение завершено перед концом периода лечения (336 часов). Однако очевидно из более высоких профилей процентили, что обработки не успешны для всех пациентов. 95-е и 90-е профили процентили также указывают, что дозирование более низкой суммы с более высокой частотой (250 tid) является более эффективным, чем менее частое дозирование с более высокой суммой (500 предложений).

hax1(1) = subplot(2,2,1)
plotCFUCount(tObs, log10CFU_250bid, 'a. Dose 250 bid' )
hax1(2) = subplot(2,2,2)
plotCFUCount(tObs, log10CFU_250tid, 'b. Dose 250 tid' )
hax1(3) = subplot(2,2,3)
plotCFUCount(tObs, log10CFU_500bid, 'c. Dose 500 bid' )
hax1(4) = subplot(2,2,4)
plotCFUCount(tObs, log10CFU_500tid, 'd. Dose 500 tid' )

% Link subplot axes
linkaxes(hax1)

hax1 = 

  Axes with properties:

             XLim: [0 1]
             YLim: [0 1]
           XScale: 'linear'
           YScale: 'linear'
    GridLineStyle: '-'
         Position: [0.1300 0.5838 0.3347 0.3412]
            Units: 'normalized'

  Use GET to show all properties


hax1 = 

  1×2 Axes array:

    Axes    Axes


hax1 = 

  1×3 Axes array:

    Axes    Axes    Axes


hax1 = 

  1×4 Axes array:

    Axes    Axes    Axes    Axes

Эффект почечной функции на антибактериальном действии

Наконец, авторы сравнили профили эффективности режимов доз в зависимости от почечной функции. Они классифицировали пациентов в четыре почечных функциональных группы на основе уровня раскрываемости преступлений креатинина (CrCL):

Creatinine Clearance Group 1: CrCL < 30
Creatinine Clearance Group 2: 30 <= CrCL < 50
Creatinine Clearance Group 3: 50 <= CrCL < 70
Creatinine Clearance Group 4: CrCL >= 70

Следующий рисунок показывает эффект почечной функции (уровень раскрываемости преступлений креатинина) на антибактериальной эффективности четырех режимов дозы. Заметьте это в нормальной почечной функциональной группе (CrCL > = 70), профили эффективности этих четырех стратегий лечения существенно отличаются друг от друга. В этом случае 500 мг t.i.d. доза являются намного более эффективными, чем другие режимы. В отличие от этого симуляции включающие пациенты с почечной дисфункцией (CrCL <30 и 30 <= CrCL <50), мы не видим большого различия между контрольными группами. Это указывает, что для пациентов с почечной дисфункцией, менее интенсивная или менее частая стратегия дозирования работала бы почти, а также стратегия дозирования с более высокой частотой или дозированием суммы.

% Preallocate
idCrCLGrp   = false(nPatients, nDoseGrps) ;

% Line Style
ls = {'bd:', 'b*:', 'rd:', 'r*:'} ;

titleStr = {'CL_c_r < 30'           , ...
            '30 <= CL_c_r < 50'     , ...
            '50 <= CL_c_r < 70'     , ...
            'CL_c_r > 70'            };

f = figure;
f.Color = 'w'

for iCrCLGrp = 1:4 % Creatinine Clearance Groups

    hax2(iCrCLGrp) = subplot(2,2, iCrCLGrp) ;
    title( titleStr{iCrCLGrp}   ) ;
    ylabel('Prob(log10CFU < 0)' ) ;
    xlabel('Time (hours)'       ) ;

end

% Set axes properties
set(hax2,   'XTick'        , 0:48:336      , ...
            'XTickLabel'   , 0:48:336      , ...
            'Ylim'         , [0 1]         , ...
            'Xlim'         , [0 336]       , ...
            'NextPlot'     , 'add'         , ...
            'Box'          , 'on'          );

% Plot results by renal function group:
for iDoseGrp = 1:nDoseGrps

    % Extract indices for renal function
    idCrCLGrp(:, 1) = CrCL(:,iDoseGrp) < 30                                    ;
    idCrCLGrp(:, 2) = CrCL(:,iDoseGrp) >= 30 & CrCL(:,iDoseGrp) < 50           ;
    idCrCLGrp(:, 3) = CrCL(:,iDoseGrp) >= 50 & CrCL(:,iDoseGrp) < 70           ;
    idCrCLGrp(:, 4) = CrCL(:,iDoseGrp) >= 70                                   ;

    for iCrCLGrp = 1:4 % Creatinine Clearance Groups

        % Calculate probability
        Pr = sum( ( log10CFU{iDoseGrp}(:, idCrCLGrp(:, iCrCLGrp)') < 0) , 2 )/sum(idCrCLGrp(:,iCrCLGrp)) ;

        % Plot
        plot(hax2(iCrCLGrp), tObs, Pr , ls{iDoseGrp}, 'MarkerSize', 7)
    end

end

legend(hax2(4), {'250 b.i.d.', '250 t.i.d.', '500 b.i.d.', '500 t.i.d.'} )
legend location NorthWest
legend boxoff

linkaxes(hax2)

f = 

  Figure (1) with properties:

      Number: 1
        Name: ''
       Color: [1 1 1]
    Position: [680 678 560 420]
       Units: 'pixels'

  Use GET to show all properties

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