Этот пример показывает, как выполнить симуляцию Монте-Карло фармакокинетической/фармакодинамической модели (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-гальта.
Входные параметры к функции 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-гальта)
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
