Этот пример показывает, как выполнить моделирование Монте-Карло фармакокинетической/фармакодинамической (PK/PD) модели для антибактериального агента. Этот пример адаптирован из Katsube et al. [1] В этом примере также показано, как использовать SimBiology ®SimFunction объект для параллельного сканирования параметров.
В этом примере требуется Toolbox™ статистики и машинного обучения. Производительность можно повысить при наличии программного обеспечения Parallel Computing Toolbox™.
Кацубе и др. [1] использовали метод моделирования и моделирования PK/PD для определения наиболее эффективных режимов дозирования для дорипенема, антибиотика карбапенема. Их исследование преследовало следующие цели:
Разработать модель PK/PD для описания антибактериального эффекта дорипенема против нескольких штаммов Pseudomonas aeruginosa
Используйте моделирование Монте-Карло, чтобы сравнить эффективность четырех общих режимов дозирования антибиотиков и определить наиболее эффективную стратегию дозирования
Исследовать влияние функции почек на антибактериальную эффективность лечения
В этом примере мы осуществим антибактериальную модель PK/PD, разработанную Katsube et al. [1] в SimBiology ® и воспроизвести результаты моделирования Монте-Карло, описанные в их работе.
[1] Т. Кацубе, Я. Яно, Т. Вадзима, Я. Ямано и М. Такано. Фармакокинетическое/фармакодинамическое моделирование и моделирование для определения эффективных режимов дозирования дорипенема. Журнал фармацевтических наук (2010) 99 (5), 2483-91.
Кацубе и др. предполагали двухкамерную инфузионную модель с линейным выведением из центрального отделения для описания фармакокинетики дорипенема. Для модели роста бактерий они предположили, что общая популяция бактерий состоит из чувствительных к лекарствам растущих клеток и нечувствительных к лекарствам покоящихся клеток. Антибактериальный эффект препарата был включен в скорость уничтожения бактерий с помощью простой модели типа Emax:
![$$ Killing Rate = \frac{Kmax*[Drug]*[Growing]}{KC50 + [Drug]} $$](../../examples/simbio/win64/antibioticdemo_eq06411219986050949670.png)
где [Drug] - концентрация (мкг/мл) лекарственного средства в центральном отделении, и [Growing] - количество растущей бактериальной популяции в КОЕ/мл (КОЕ = колониеобразующие единицы). Kmax - максимальная константа скорости уничтожения (1/час) и KC50 является константой скорости Михаэлиса-Ментена (мкг/мл).
Графическое представление реализации модели SimBiology показано ниже.
% Load model sbioloadproject('AntibacterialPKPD.sbproj', 'm1') ;
Кацубе и др. смоделировала модель с использованием четырех общих стратегий дозировки антибиотиков.
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Предполагалось, что линейно зависит от скорости клиренса креатинина через следующее уравнение:
где
- остаточная ошибка добавки, отобранная из нормального распределения со средним значением 0 мл/мин и стандартным отклонением 22 мл/мин.
CL = 1.07*CrCL + 45.6 + normrnd(0,22,nPatients,nDoseGrps); % units: ml/minute
Распределение фармакодинамических (PD) параметров:
Константы скорости преобразования с ростом до покоя, k1 и k2, были отобраны из логнормального распределения с типичным значением 5,59-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
Кацубе и др. предполагали, что значения k1, k2 и Kmax были независимы от обрабатываемого бактериального штамма. Значение Beta, постоянная чистая скорость роста была зафиксирована на уровне 1,5 1/час.
Основываясь на экспериментах с несколькими штаммами, авторы пришли к выводу, что ценность KC50 линейно зависела от минимальной концентрации ингибирования (MIC) бактериального штамма посредством следующего уравнения.
где - остаточная ошибка добавки, отобранная из нормального распределения со средним значением 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 ®. На вкладке Главная в разделе Среда выберите Параллельный > Обнаружить кластеры. Чтобы получить доступ к этим кластерам, необходимо предоставить данные для входа в учетную запись MathWorks ®. Дополнительные сведения см. в разделах Обнаружение кластеров и использование профилей кластеров (панель инструментов параллельных вычислений).
Создайте параллельный пул, если он не существует.
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. Кацубе и др. установить критерий элиминации бактерий при 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
