В этом примере показано, как выполнить симуляцию Монте-Карло фармакокинетической/фармакодинамической модели (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.
Katsube и др. принял модель вливания 2D отсека с линейным устранением из центрального отсека, чтобы описать фармакокинетику doripenem. Для бактериальной модели роста они приняли, что общее бактериальное население состоит из восприимчивых к препарату растущих ячеек и нечувствительных к препарату ячеек отдыха. Антибактериальный эффект препарата был включен в уровень уничтожения бактерий с помощью простой модели типа Emax:
где [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
, был принят, чтобы зависеть линейно от уровня раскрываемости преступлений креатинина через следующее уравнение:
где аддитивная остаточная ошибка производится от нормального распределения со средним значением 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) бактериального штамма через следующее уравнение.
где аддитивная остаточная ошибка производится от нормального распределения со средним значением 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
Создайте 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