Exponenta Banner

Симуляция PK/PD для руководства стратегией дозирования антибиотиков

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

Открыть скрипт

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

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

Фон

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

  • Разработайте модель PK/PD, чтобы описать антибактериальный эффект дорипенема против нескольких штаммов Pseudomonas aeruginosa

  • Используйте симуляции Монте-Карло, чтобы сравнить эффективность четырех распространенных режимов дозировки антибиотиков и определить наиболее эффективную стратегию дозирования

  • Исследуйте эффект функции почек на антибактериальную эффективность лечения

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

Ссылки

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

Модель PK/PD

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

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

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

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

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

Режимы дозировки

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

  • 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. Тип распределения и значения параметров распределения основывались на данных предыдущих клинических испытаний дорипенема, проведенных в Японии.

Примечание: В [1] в каждой группе дозировки моделировали 5000 виртуальных пациентов. В этом примере мы будем использовать 1000 пациентов в каждой группе. Чтобы симулировать другой размер населения, измените значение 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) вычисляли с помощью уравнения Кокрофта-Голта.

Входы в 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

Клиренс креатинина (с помощью уравнения Кокрофта-Голта)

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

Распределение фармакокинетических (ПК) параметров:

Параметры 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, были отобраны из lognormal распределения с типичным значением 5,59e-5 и 0,0297 1/час, соответственно, каждый с CV 20%. Kmax отбирали из логнормального распределения с типичным значением 3,5 1/час и 15,9% CV.

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 et al. принял, что значения k1, k2 и Kmax были независимы от обрабатываемого бактериального штамма. Значение Betaчистую константу скорости роста фиксировали на уровне 1,5 1/час.

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

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

где$\varepsilon$ - аддитивная остаточная ошибка, выбранная из нормального распределения со средним значением 0 и стандартным отклонением 1,06 мкг/мл. В симуляции 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

Настройка и проектирование симуляции

Создайте 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 true, чтобы включить параллельные вычисления.

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 ® Account. Для получения дополнительной информации смотрите Обнаружение кластеров и Использование профилей кластеров (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 et al. установите критерий для бактериальной элиминации в log10(CFU) < 0, где CFU - общее количество бактерий. Эффективность каждого режима дозы измеряли как долю населения, которая достигала критериев успеха в группе доз. Эта метрика эффективности, Pr{log10(CFU) < 0}, отслеживали как функцию времени для каждой дозированной группы.

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

  • Умеренная инфекция (Начальное количество бактерий = 1e4 КОЕ/мл)

  • Тяжелая инфекция (Начальное количество бактерий = 1e7 КОЕ/мл)

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

% 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 тид) является более эффективным, чем менее частое дозирование с более высоким количеством (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):

  • Группа клиренса креатинина 1: CrCL < 30

  • Группа клиренса креатинина 2:30 < = CrCL < 50

  • Группа клиренса креатинина 3:50 < = CrCL < 70

  • Группа клиренса креатинина 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

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

Поддержка