PK/PD建模与仿真指导抗生素给药策略
本例演示如何对抗菌剂的药代动力学/药效学(PK/PD)模型进行蒙特卡罗模拟。这个例子改编自Katsube等人。[1]这个例子还展示了如何使用SimBiology®SimFunction对象并行执行参数扫描。
这个例子需要统计和机器学习工具箱™。如果您拥有并行计算工具箱™软件，则可以提高性能。
背景
Katsube等人[1]使用PK/PD建模和模拟方法确定碳青霉烯类抗生素多利培南(doripenem)的最有效剂量方案。他们的研究目标是:
建立PK/PD模型来描述多oripenem对几种铜绿假单胞菌的抗菌作用
使用蒙特卡罗模拟比较四种常见抗生素剂量制度的疗效，并确定最有效的剂量策略
观察肾功能对药物抗菌效果的影响
在本例中，我们将在SimBiology®中实现Katsube等人开发的抗菌药物PK/PD模型，并复制他们工作中描述的蒙特卡罗模拟的结果。
参考文献
Katsube, Y.矢野，T.和岛，Y.山野和M.高野。药代动力学/药效学建模和模拟，以确定多利培南的有效剂量方案。药学杂志(2010) 99(5)， 2483-91。
PK / PD模型
Katsube等人假设了一个双室输注模型，从中央室线性消除，以描述多oripenem的药代动力学。对于细菌生长模型，他们假设细菌总数由药物敏感的生长细胞和药物不敏感的静息细胞组成。通过简单的Emax型模型，将药物的抗菌作用包括在细菌的杀灭率中:
![$$杀伤率= \frac{Kmax*[药物]*[生长]}{KC50 +[药物]}$$](http://www.tatmou.com/au/au/help/examples/simbio/win64/antibioticdemo_eq06411219986050949670.png)
在哪里(药物)药物的浓度(ug/ml)在中央室，和(增长)为不断增长的细菌数量，单位为CFU/ml (CFU =菌落形成单位)。Kmax最大杀伤速率是常数(1/小时)和KC50为米克利斯-门腾速率常数(ug/ml)。
模型的SimBiology实现的图形化视图如下所示。
%负载模型sbioloadproject (“AntibacterialPKPD.sbproj”，“m1”);
剂量方案
Katsube等人使用四种常见的抗生素剂量策略模拟了该模型。
250毫克，每日2次(每日口服)
250毫克，每日三次(定时注射)
500毫克，每日2次(口服)
500毫克，每日三次(定时注射)
所有四种给药方案均采用输注给药，输注时间设置为30分钟。在SimBiology中，这些剂量方案已经作为剂量对象实现。
选择模型中的剂量对象doseNames = {“250毫克竞标”，“250毫克tid”，“500毫克bid”，“500毫克tid”};为iDoseGrp = 1:length(doseNames) doseRegimens(iDoseGrp) = sbioselect(m1，“名字”， doseNames{iDoseGrp});结束
虚拟人口描述
基于人口统计变量的分布和PK/PD参数，生成个体的虚拟种群。分布类型和分布参数的值基于日本进行的多立培南早期临床试验的数据。
注:[1]中，每个给药组模拟虚拟患者5000人。在本例中，我们将在每组中使用1000名患者。要模拟不同的总体大小，请更改的值nPatients在下面。
%设置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^2/√(v+m^2));= @(m,v)√(log(v/m^2+1));m = @(typicalValue);v = @(typicalValue,CV) typicalValue²*CV²;患者人口统计rng (“默认”）;Wt = normrnd(51.6, 11.8, nPatients, nDoseGrps);%单位:kg年龄= normrnd(71.8, 11.9, nPatients, nDoseGrps);%单位:年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);%单位:ml/分钟性别比率id = 1:nPatients*nDoseGrps;idFemale = randsample(id, round(0.26*nDoseGrps*nPatients));% 26%女性
肌酐清除(使用Cockcroft-Gault方程)
CrCL =(140 -年龄).*Wt./(Scr*72);%单位:ml/分钟CrCL(idFemale) = CrCL(idFemale)*0.85;%乘以0.85为女性
药代动力学(PK)参数分布:
PK参数,中央，k12的,k21里面，均采对数正态分布，典型值分别为7.64升、1.59 1升/小时和2.26 1升/小时，变异系数(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 = lognrnd(Central_mu, Central_sigma, nPatients, nDoseGrps);%单位:升k12 = lognrnd(k12_mu, k12_sigma, nPatients, nDoseGrps);%单位:1/小时k21 = lognrnd(k21_mu, k21_sigma, nPatients, nDoseGrps);%单位:1/小时
药物清除率，CL，由下式假设与肌酐清除率线性相关:
在哪里
为从正态分布中抽样的加性残差，均值为0 ml/min，标准差为22 ml/min。
CL = 1.07*CrCL + 45.6 + normrnd(0,22,nPatients,nDoseGrps);%单位:ml/分钟
药效学参数分布:
生长到静止的转换速率常数，k1而且k2，均采对数正态分布，典型值分别为5.59e-5和0.0297 1/hour, 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);%单位:1/小时k2 = lognrnd(k2_mu, k2_sigma, nPatients, nDoseGrps);%单位:1/小时Kmax = lognrnd(Kmax_mu, Kmax_sigma, nPatients, nDoseGrps);%单位:1/小时
Katsube等人假设k1，k2而且Kmax与被治疗的细菌菌株无关。的价值β，净生长速率常数，固定为1.5 1/h。
通过对几种菌株的实验，作者得出结论KC50与菌株的最低抑菌浓度(MIC)呈线性相关，由下式可知。
在哪里是均值为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];%来自离散分布的抽样MIC值，使用randsampleMIC = nan(nPatients, nDoseGrps);% preallocate为iDoseGrp = 1:nDoseGrps MIC(:， iDoseGrp) = randsample(micValue, nPatients, true, micFreq);结束KC50 = exp(-1.91 + 0.898*log(MIC) + 1.06*randn(nPatients, nDoseGrps));%单位:微克/毫升
仿真设置与设计
创建一个SimFunction对象，使您可以并行执行模型模拟和参数扫描。在这个例子中，你将改变8个参数，中央，k12的，k21里面，CL，k1，k2，Kmax,KC50．选择生长和静止的细菌计数，日益增长的而且休息，作为响应，即您希望在改变这些输入参数时观察到的模拟结果。
选择需要变化的参数。
参数= {“中央”，“k12的”，k21里面的，“氯”，“k1”，“k2”，“Kmax”，“KC50”};
选择回答。
可观测值= {[细菌生长模型]。日益增长的，.．.[细菌生长模型]。休息的};
设定一个模板剂量。
Tempdose = sbiodose(“剂量”）;tempdose。目标=的中央。药物的；tempdose。AmountUnits =毫克的；tempdose。TimeUnits =“小时”；tempdose。DurationParameterName =“TDose”；
创建一个SimFunction对象。集UseParallel为true以启用并行计算。
simfunc = createSimFunction(m1,params,observables,tempdose，“UseParallel”,真正的);
创建一个输入矩阵φ对于每个剂量组。
phi = cell(1,nDoseGrps);为i = 1: nDoseGrpsφ{我}=[中央(:,i)、k12(:,我),k21(:,我).．.CL(:，i) k1(:，i) k2(:，i)，.．.Kmax(:,我),KC50 (:, i)];结束
集群计算
本示例使用预先配置到本地计算机的本地集群概要文件。您还可以搜索在Amazon EC2®上运行的其他MATLAB®Parallel Server™集群。在首页页中的环境部分中,选择平行>发现集群．若要访问这些集群，您必须提供MathWorks®帐户登录信息。详细信息请参见发现集群并使用集群概要文件(并行计算工具箱)．
如果不存在并行池，请创建并行池。
如果parpool isempty (gcp);结束
使用“本地”配置文件启动并行池(parpool)…连接到并行池(工人数:6)。
对于所有剂量方案，模拟模型直到t =第一次剂量后2周。在整个给药过程中，每24小时(每天1次)采样一次总细菌计数(CFU)。
tObs = 0:24:336;%小时nTPoints =长度(tObs);%抽样点数
重症感染患者的蒙特卡罗模拟
一种药物的抗菌效果可以通过不同的PK/PD指数来衡量。Katsube等人将细菌消除的标准设定为log10(CFU) < 0，其中CFU为细菌总数。每种给药方案的疗效以剂量组中达到成功标准的人群比例来衡量。这个效能指标，Pr{log10(CFU) < 0}，作为每个剂量组的时间函数进行跟踪。
在他们的模拟研究中，作者调查了两类患者的剂量方案的疗效:
中度感染(初始细菌计数= 1e4 CFU/ml)
严重感染(初始细菌计数= 1e7 CFU/ml)
在本例中，我们只复制严重感染病例的结果。请注意，您可以很容易地模拟另一种情况，即中度感染的患者，通过改变初始的细菌计数(即细菌计数)日益增长的种)，在模型中为1e4 CFU/ml。
% Preallocatecfu = nan(nTPoints,nPatients);log10CFU = cell(1,nDoseGrps);为i = 1:nDoseGrps disp([“模拟组”num2str(我),“……”])直接从每个剂量的现有剂量对象获取剂量表%给药方案。doseTable = getttable (doseRegimens(i));%模拟simdata = simfunc(phi{i}，[]，doseTable,tObs);每个病人生长和静息细菌计数的总和为j = 1:nPatients cfu(:，j) = sum(simdata(j).Data,2);结束存储每个剂量组的日志转换计数。log10CFU{i} = log10(cfu);结束%保存结果log10CFU_250bid = log10CFU{1};log10CFU_250tid = log10CFU{2};log10CFU_500bid = log10CFU{3};log10CFU_500tid = log10CFU{4};
模拟组1…模拟组2…模拟组3…模拟组4…
关闭并行池。
删除(gcp (“nocreate”));
细菌计数的时间过程剖面
我们绘制中位数(红色)和百分位数(阴影)的概况log10 (CFU)所有四种给药方案的水平。观察到，在所有四组中，中位时间病程剖面显示，细菌在治疗期(336小时)结束前完全根除。然而，从较高百分比的资料可以明显看出，治疗并不是对所有患者都成功。第95和90百分位曲线也表明，低剂量高频率(250 tid)比低频率高剂量(500 bid)更有效。
hax1(1) = subplot(2,2,1) plotCFUCount(tObs, log10CFU_250bid，“a.剂量250”) hax1(2) = subplot(2,2,2) plotCFUCount(tObs, log10CFU_250tid，'b.剂量250 tid') hax1(3) = subplot(2,2,3) plotCFUCount(tObs, log10CFU_500bid，“c.剂量500出价”) hax1(4) = subplot(2,2,4) plotCFUCount(tObs, log10CFU_500tid，d.剂量500 tid）链接子图轴linkaxes (hax1)
hax1 =带有属性的坐标轴:XLim: [0 1] YLim: [0 1] XScale: 'linear' YScale: 'linear' GridLineStyle: '-'位置:[0.1300 0.5838 0.3347 0.3412]单位:'normalized'使用GET显示所有属性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)，四种治疗策略的疗效差异显著。在这种情况下，500mg t.i.d.剂量比其他方案更有效。相比之下，涉及肾功能不全患者的模拟(CrCL< 30和30 <=CrCL< 50)时，治疗组间差异不大。这表明，对于肾功能不全的患者，较低强度或较低频率的给药策略几乎与较高频率或给药量的给药策略一样有效。
% PreallocateidCrCLGrp = false(nPatients, nDoseGrps);%线条样式Ls = {的双相障碍:，" b *:，“理查德·道金斯:”，的r *:};titleStr = {'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));%的阴谋plot(hax2(iCrCLGrp)， tObs, Pr, ls{iDoseGrp}，“MarkerSize”7)结束结束传奇(hax2 (4) {“250 b.i.d。”，“250 t.i.d。”，“500 b.i.d。”，“500 t.i.d。”})传说位置西北传说boxofflinkaxes (hax2)
f =图(1)与属性:数字:1名称:"颜色:[1 1 1]位置:[680 678 560 420]单位:'像素'使用GET显示所有属性
您也可以从以下列表中选择一个网站:
