指导抗生素给药策略的PK/PD建模与仿真

这个例子演示了如何对抗菌剂的药代动力学/药效学(PK/PD)模型进行蒙特卡洛模拟。这个例子还展示了如何使用SimBiology®SimFunction对象以并行执行参数扫描。

这个例子需要Statistics和Machine Learning Toolbox™。如果您有并行计算工具箱™软件,则可以提高性能。

出身背景

Katsube等人[1]使用PK/PD建模和模拟方法确定碳青霉烯类抗生素多利培南的最有效剂量方案。他们的研究目标是:

  • 建立一个PK/PD模型来描述多利培南对几种铜绿假单胞菌菌株的抗菌作用

  • 使用蒙特卡洛模拟来比较四种常见抗生素剂量方案的效果,并确定最有效的剂量策略

  • 观察肾功能对抗菌治疗效果的影响

在本例中,我们将在SimBiology®中实现Katsube等人[1]开发的抗菌PK/PD模型,并复制他们工作中描述的蒙特卡罗模拟结果。

参考文献

[1] T.Katsube,Y.Yano,T.Wajima,Y.Yamano和M.Takano.确定多利培南有效剂量方案的药代动力学/药效学建模和模拟。制药科学杂志(2010) 99(5), 2483 - 91。

PK / PD模型

Katsube等人假设了一个由中央室线性消除的双室输注模型来描述多立培南的药代动力学。对于细菌生长模型,他们假设总细菌群由药物敏感的生长细胞和药物不敏感的静息细胞组成。通过简单的Emax模型,将药物的抗菌作用纳入对细菌的杀灭率:

$$ kill Rate = {Kmax*[Drug]*[Growing]}{KC50 + [Drug]} $$

哪里[药物]药物在中央腔室的浓度(ug/ml)和[成长]为CFU/ml中不断增长的细菌数量(CFU =菌落形成单位)。Kmax最大杀灭率常数(1/h)和KC50为Michaelis-Menten速率常数(ug/ml)。

模型的SimBiology实现的图形视图如下所示。

%负荷模型sbioloadproject (“AntibacterialPKPD.sbproj”“m1”);

剂量方案

Katsube等人使用四种常见的抗生素剂量策略模拟了该模型。

  • 每日两次250毫克(每日一次)

  • 每日三次,每次250毫克(t.i.d)

  • 500毫克,每天两次(b.i.d.)

  • 500毫克,每天三次(每日三次)

在所有四种剂量方案中都使用了输注剂量,输注时间设置为30分钟。在SimBiology中,这些剂量方案已作为剂量对象实施。

%在模型中选择剂量对象doseNames={“250毫克报价”‘每日三次每次250毫克’“500毫克出价”“500毫克tid”};iDoseGrp=1:length(doseNames)doseRegimens(iDoseGrp)=sbioselect(m1,“名字”,doseNames{iDoseGrp});结束

虚拟人口的描述

基于人口统计学变量和PK/PD参数的分布,生成了一个虚拟的个体群体。分布类型和分布参数的值基于在日本进行的多利培南早期临床试验的数据。

注:在[1]中,每个剂量组模拟5000例虚拟患者。在本例中,我们将在每组中使用1000名患者。要模拟不同的种群大小,请更改的值病人在下面

%设置nPatients = 1000;%每剂量组患者人数nDoseGrps = 4;%试验剂量方案数

人口统计变量分布:

体重(Wt)及年龄(年龄)为正态分布,均值分别为51.6 kg和71.8年,标准差分别为11.8 kg和11.9年。26%的人口被认为是女性。血清肌酐水平(可控硅)从对数正态分布中取样,典型值(几何平均值)为0.78 mg/dL,变异系数(CV)为32.8%。肌酐清除率(CrCL)的计算采用Cockcroft-Gault方程。

输入到lognrnd函数是平均值(μ)及标准偏差(σ)的相关正态分布。在这个例子中,μσ由所报告的对数正态分布的典型值和变异系数计算。您可以使用以下定义来计算它们。看到lognstat有关更多信息,请参阅文档。

= @(m,v) log(m²/√(v+m²))= @(m,v)√(log(v/m²+1))m = @(typicalValue) typicalValue;v = @(typicalValue,CV) typicalValue^2*CV^2;病人的人口百分比rng (“默认”)Wt=normrnd(51.6,11.8,N患者,nDoseGrps);%单位:千克年龄=正常(71.8,11.9,nPatients, nDoseGrps);%单位:年sc_mu = mu(m(0.78), v(0.78,0.328));sc_sigma = sigma(m(0.78), v(0.78,0.328));Scr = logrnd (Scr_mu, Scr_sigma, nPatients, nDoseGrps);%单位:毫升/分钟性别比例%id = 1:nPatients*nDoseGrps;idFemale = randsample(id, round(0.26*nDoseGrps*nPatients));% 26%为女性

肌酐清除率(使用科克罗夫特-高尔特方程)

CrCL =(140 -年龄).*Wt./(Scr*72);%单位:毫升/分钟CrCL(idFemale) = CrCL(idFemale)*0.85;%乘以0.85的雌性

药代动力学(PK)参数分布:

PK参数,中央k12的,k21里面的样本为对数正态分布,典型值分别为7.64升、1.59 1/h和2.26 1/h,变异系数(CV)为20%。中央是中央隔间的分配量,和k12的k21里面传递速率常数是中央外围设备隔间。

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 = logrnd (Central_mu, Central_sigma, nPatients, nDoseGrps);%单位:升k12=lognrnd(k12μm,k12μsigma,N患者,NDOSEGRP);%单位:1 /小时k21 = logrnd (k21_mu, k21_sigma, nPatients, nDoseGrps);%单位:1 /小时

药物清除率,CL,通过以下方程假设其与肌酐清除率线性相关:

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

哪里$\varepsilon$为正态分布抽样的平均值为0 ml/min,标准差为22 ml/min的可加性残差。

CL=1.07*CrCL+45.6+正常值(0,22,N患者,NDOSEGRP);%单位:毫升/分钟

药效学(PD)参数分布:

生长-静止转换速率常数,k1k2,样本为对数正态分布,典型值分别为5.59e-5和0.0297 1/h, CV均为20%。Kmax为对数正态分布,典型值为3.5 1/h, 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 = logrnd (k1_mu, k1_sigma, nPatients, nDoseGrps);%单位:1 /小时k2 = logrnd (k2_mu, k2_sigma, nPatients, nDoseGrps);%单位:1 /小时Kmax = logrnd (Kmax_mu, Kmax_sigma, nPatients, nDoseGrps);%单位:1 /小时

Katsube等人假设了这个值k1k2Kmax与被处理的菌株无关。的价值β,净生长速率常数,固定在1.5 1/h。

基于对几种菌株的实验,作者得出结论KC50与菌株的最小抑制浓度(MIC)呈线性关系,公式如下:

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

哪里$\varepsilon$是从正态分布中采样的加性残差,平均值为0,标准偏差为1.06 ug/ml。在模拟中麦克风从离散分布中采样值,然后KC50为所选对象计算值麦克风用上面的方程。

% 71株铜绿假单胞菌MIC值的离散分布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];%使用randsample从离散分布中取样MIC值MIC=nan(n患者,nDOSEGRP);% preallocateiDoseGrp = 1:nDoseGrps MIC(:, iDoseGrp) = randsample(micValue, nPatients, true, micFreq);结束KC50=exp(-1.91+0.898*log(MIC)+1.06*randn(n患者,nDoseGrps));%单位:微克/毫升

模拟装置与设计

创建一个SimFunction对象,该对象允许您并行执行模型模拟和参数扫描。在本例中,您将更改8个参数,中央k12的k21里面CLk1k2Kmax,KC50.选择生长和静止细菌计数,增长的休息,作为响应,即在更改这些输入参数时要观察的模拟结果。

选择要更改的参数。

params = {“中央”“k12的”k21里面的“CL”“k1”“k2”“Kmax”“KC50”};

选择响应。

可见= {“细菌生长模型。日益增长的...“[细菌生长模型]。静息”};

设定一个模板剂量。

tempdose=sbiodose(“剂量”)目标=的中央。药物的;tempdose.AmountUnits=“毫克”;tempdose。TimeUnits =“小时”;tempdose。DurationParameterName =“TDose”

创建一个SimFunction对象。集UseParallel为true以启用并行计算。

simfunc = createSimFunction (m1, params,可见,tempdose,“UseParallel”,对);

创建一个输入矩阵φ每一剂量组。

φ=细胞(1、nDoseGrps);i = 1: nDoseGrpsφ{我}=[中央(:,i)、k12(:,我),k21(:,我)...CL(:,我),k1 (:, i), k2(:,我)...Kmax(:,我),KC50 (:, i)];结束

集群计算

这个示例使用预先配置到本地机器的本地集群配置文件。您还可以搜索其他在Amazon EC2®上运行的MATLAB®Parallel Server™集群。在首页选项卡中环境节,选择平行>发现集群.要访问这些集群,必须提供MathWorks®帐户登录信息。有关详细信息,请参见发现集群并使用集群配置文件(并行计算工具箱)。

如果不存在并行池,则创建并行池。

如果parpool isempty (gcp);结束
使用“local”配置文件启动并行池(parpool)…连接到并行池(工作人员数量:6)。

对于所有剂量场景,模拟模型至第一次给药时间t = 2周。在整个给药期间,每24小时(每天一次)采样一次细菌总数,即CFU。

tObs = 0:24:336;%小时nTPoints = length(tObs);%采样点数量

严重感染患者的蒙特卡罗模拟

通过不同的PK/PD指标可以测定药物的抗菌效果。Katsube等人在log10 (CFU) < 0,其中CFU是细菌总数。每个剂量方案的有效性是通过在该剂量组中达到成功标准的人口比例来衡量的。这个功效指标,Pr{log10(CFU)<0},作为各剂量组时间的函数进行跟踪。

在他们的模拟研究中,作者调查了两类患者的剂量方案的效果:

  • 中度感染(初始细菌计数= 1e4 CFU/ml)

  • 严重感染(初始细菌计数= 1e7 CFU/ml)

在本例中,我们只复制严重感染病例的结果。注意,你可以很容易地模拟另一种情况,中度感染的患者,通过改变初始的细菌计数(the增长的种),在模型中为1e4 CFU/ml。

% Preallocatecfu=nan(nTPoints,NP患者);log10CFU=cell(1,nDoseGrps);i=1:nDoseGrps disp([“模拟组”num2str(我),'...'])%直接从现有剂量对象获取剂量表%给药方案。doseTable =可以获得的(doseRegimens (i));%模拟simdata = simfunc(φ{我},[],doseTable则);%每名患者生长和静止细菌计数之和j=1:n患者cfu(:,j)=和(simdata(j)。数据,2);结束%存储每个剂量组的日志转换计数。log10CFU{我}= log10 (cfu);结束%保存结果log10CFU_250bid = log10CFU{1};log10CFU_250tid = log10CFU{2};log10CFU_500bid = log10CFU{3};log10CFU_500tid = log10CFU{4};
模拟第一组…模拟第二组…模拟第三组…模拟第四组。。。

关闭并行池。

删除(gcp)(“nocreate”));

细菌计数的时程曲线

我们绘制了该区域的中位数(红色)和百分位(阴影)剖面图日志10(CFU)所有四种剂量方案的水平。观察所有四组的中位时间过程曲线显示,在治疗期结束前(336小时)细菌已完全根除。然而,从较高的百分位数中可以明显看出,治疗并不是对所有患者都成功。第95和90百分位分布也表明,较低剂量且较高频率(250 tid)的剂量比较低频率且较高剂量(500 bid)的剂量更有效。

hax1(1)=子标段(2,2,1)标绘量(tObs,log10CFUúu标,“a. 250剂量”) hax1(2) = subplot(2,2,2)“b.每天三次每次250剂”)hax1(3)=子标段(2,2,3)标绘量(tObs,log10CFU\U 500bid,“c.剂量500标”) hax1(4) = subplot(2,2,4)d.剂量为每天三次每次500粒%链接子地块轴linkaxes (hax1)
hax1 =轴与属性:XLim: [0 1] YLim: [0 1] XScale:“线性”YScale:“线性”GridLineStyle:“-”位置:[0.1300 0.5838 0.3347 0.3412)单位:“正常化”的使用可以显示所有属性hax1 = 1×2轴阵列:轴轴hax1 = 1×3轴阵列:轴轴轴hax1 = 1×4轴阵列:轴轴轴轴

肾功能对抗菌活性的影响

最后,作者比较了不同剂量方案对肾功能的影响。他们根据肌酐清除率将患者分为四个肾功能组(CrCL):

  • 肌酐清除组1:CrCL< 30

  • 肌酐清除组2:30<=CrCL< 50

  • 肌酐清除组3:50 <=CrCL< 70

  • 肌酐清除组4:CrCL> = 70

下图为肾功能(肌酐清除率)对四种给药方案抗菌效果的影响。观察正常肾功能组(CrCL>= 70),四种治疗策略的疗效表现有显著差异。在这种情况下,500毫克每分钟的剂量比其他方案更有效。与此相反,模拟涉及肾功能不全患者(CrCL<30和30<=CrCL< 50),我们没有看到治疗组之间的太大差异。这表明,对于肾功能不全的患者,低强度或低频率的给药策略与高频率或高剂量的给药策略几乎一样有效。

% PreallocateidCrCLGrp = false(nPatients, nDoseGrps);%线条样式ls = {的双相障碍:" b *:“理查德·道金斯:”的r *:};标题={“CL_c_r < 30”...'30 <= CL_c_r < 50'...'50 <= CL_c_r < 70'...“CL_c_r > 70年};f =图;f.Color =' w 'iCrCLGrp=1:4%肌酐清除组hax2(iCrCLGrp) = subplot(2,2, iCrCLGrp);title(titleStr{iCrCLGrp});ylabel (“概率(log10CFU < 0)”);包含(的时间(小时));结束设置坐标轴属性集(hax2,“XTick”, 0:48:336 ,...“XTickLabel”, 0:48:336 ,...“Ylim”, [0 1],...“Xlim”, [0 336] ,...“NextPlot”“添加”...“盒子”“上”);%按肾功能组绘制结果:iDoseGrp = 1: nDoseGrps%提取肾功能指标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;iCrCLGrp=1:4%肌酐清除组%计算概率Pr = sum((log10CFU{iDoseGrp}(:, idCrCLGrp(:,iCrCLGrp) ') < 0), 2)/sum(idCrCLGrp(:,iCrCLGrp));%的阴谋图(hax2(iCrCLGrp), tob, Pr, ls{iDoseGrp},“MarkerSize”7)结束结束图例(hax2(4){“250 b.i.d。”“250 t.i.d。”“500 b.i.d。”“500 t.i.d。”})传奇地方西北传说博科夫linkaxes (hax2)
f = Figure (1) with properties: Number: 1 Name: " Color: [1 1 1] Position: [680 678 560 420] Units: 'pixels'使用GET显示所有属性