进行PK/PD建模和仿真，指导抗生素给药策略
这个例子展示了如何对一种抗菌剂的药代动力学/药效动力学(PK/PD)模型进行蒙特卡罗模拟。本例改编自Katsube等。本例还展示了如何使用SimBiology®SimFunction对象并行执行参数扫描。
本示例需要统计和机器学习工具箱™。如果您有并行计算工具箱™软件，则可以提高性能。
背景
Katsube等人[1]使用PK/PD建模和模拟方法确定了碳青霉烯类抗生素多oripenem的最有效剂量方案。他们的研究目标是:
建立PK/PD模型来描述多oripenem对几种铜绿假单胞菌的抗菌作用
使用蒙特卡罗模拟来比较四种常用抗生素给药方案的疗效，并确定最有效的给药策略
观察肾功能对药物抗菌效果的影响
在本例中，我们将实现Katsube等人在SimBiology®中开发的抗菌PK/PD模型，并复制他们工作中描述的蒙特卡罗模拟结果。
参考文献
[1]胜部T.，矢野Y.，和岛T.，山野Y.，高野M.。确定多oripenem有效给药方案的药代动力学/药效学建模和模拟。药理学杂志(2010) 99(5), 2483 - 91。
PK / PD模型
Katsube等人假设了一个从中央室室线性消除的双室输注模型来描述多oripenem的药代动力学。对于细菌生长模型，他们假设细菌总数由药物敏感的生长细胞和药物不敏感的静止细胞组成。通过一个简单的Emax型模型，将药物的抗菌作用计入细菌的杀死率:
![$$杀伤率= \frac{Kmax*[药物]*[增长]}{KC50 +[药物]}$$](http://www.ru-cchi.com/help/examples/simbio/win64/antibioticdemo_eq16197705998259887470.png)
在哪里(药物)药物的浓度(ug/ml)是否在中央隔室(增长)为生长菌群数，单位为CFU/ml (CFU =菌落形成单位)。Kmax最大杀伤速率是常数(1/小时)和吗KC50为Michaelis-Menten速率常数(ug/ml)。
该模型的SimBiology实现的图形视图如下所示。
%负荷模型sbioloadproject (“AntibacterialPKPD.sbproj”，“m1”);
剂量方案
Katsube等人使用四种常见的抗生素剂量策略模拟了该模型。
每天两次250mg(每日一次)
每日三次，每次250毫克。
每日两次500毫克(每日一次)
每日三次，每次500毫克。
四种给药方案均采用输液给药，输液时间设置为30分钟。在SimBiology中，这些剂量方案作为剂量对象实现。
%选择模型中的剂量对象doseNames = {“250毫克报价”，“250毫克tid”，“500毫克报价”，“500毫克tid”};为iDoseGrp = 1:length(doseNames) doseRegimens(iDoseGrp) = sbioselect(m1，“名字”, doseNames {iDoseGrp});结束
虚拟人口的描述
根据人口统计变量的分布和PK/PD参数生成个体的虚拟种群。分布类型和分布参数的值基于日本早期多oripenem临床试验的数据。
注:[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²/√(v+m²));= @(m,v)√(log(v/m²+1));m = @(typicalValue);v = @(typevalue,CV) typevalue ^2*CV^2;病人的人口百分比rng (“默认”）;Wt = normrnd(51.6, 11.8, nPatients, nDoseGrps);%单位:千克年龄= 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);%单位:毫升/分钟性别比例%id = 1:nPatients*nDoseGrps;idFemale = randsample(id, round(0.26*nDoseGrps*nPatients));% 26%为女性
肌酐清除(使用Cockcroft-Gault方程)
CrCL =(140 -年龄).*Wt./(Scr*72);%单位:毫升/分钟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 - σ = σ (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);%单位:毫升/分钟
药效学(PD)参数分布:
从增长到静止的变换速率常数，k1而且k2样本为对数正态分布，典型值分别为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 = 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．选择生长和静止的细菌数量，日益增长的而且休息，作为响应，即您希望在改变这些输入参数时观察的模拟结果。
选择要更改的参数。
params = {“中央”，“k12的”，k21里面的，“氯”，“k1”，“k2”，“Kmax”，“KC50”};
选择反应。
可见= {“细菌生长模型。日益增长的，.．.“细菌生长模型。休息的};
设置模板剂量。
tempdose = sbiodose (“剂量”）;tempdose。目标=的中央。药物的；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)];结束
集群计算
本示例使用预先配置到本地计算机的本地集群配置文件。您也可以搜索其他的MATLAB®并行服务器™集群运行在Amazon EC2®。在首页选项卡中环境部分中,选择平行>发现集群．要访问这些集群，您必须提供MathWorks®帐户登录信息。有关详细信息,请参见发现集群并使用集群概要文件(并行计算工具箱)．
如果不存在并行池，则创建并行池。
如果parpool isempty (gcp);结束
使用“本地”配置文件启动并行池(parpool)…连接到并行池(工作人员数量:6)。
对于所有给药方案，模型一直模拟到t =第一次给药后的2周。在整个给药过程中，每24小时(每天一次)采样一次细菌总数CFU。
tObs = 0:24:336;%小时nTPoints =长度(tObs);%采样点的个数
重症感染患者的蒙特卡罗模拟
用不同的PK/PD指标可以测定药物的抗菌效果。Katsube等人将细菌消除的标准设定在log10 (CFU) < 0，其中CFU为细菌总数。每种剂量方案的疗效以剂量组中达到成功标准的人群比例来衡量。这个功效指标,公关{log10 (CFU) < 0}，作为每个剂量组的时间函数进行跟踪。
在他们的模拟研究中，作者调查了两类患者的剂量方案的疗效:
中度感染(初始细菌计数= 1e4 CFU/ml)
严重感染(初始细菌计数= 1e7 CFU/ml)
在本例中，我们只复制严重感染病例的结果。注意，您可以很容易地模拟另一种情况，即中度感染的患者，通过改变初始的细菌计数日益增长的种)，在模型中为1e4 CFU/ml。
% Preallocatecfu =南(nTPoints nPatients);log10CFU = cell(1,nDoseGrps);为i = 1:nDoseGrps disp([“模拟组”num2str(我),“……”])%直接从每个剂量的现有剂量对象获取剂量表%给药方案。doseTable =可以获得的(doseRegimens (i));%模拟simdata = simfunc(φ{我},[],doseTable则);%每个病人的生长和静止细菌计数之和为j = 1:nPatients cfu(:，j) = sum(simdata(j).Data,2);结束存储每个剂量组的对数转换计数。log10CFU{我}= 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 bid”) hax1(4) = subplot(2,2,4) plotCFUCount(tObs, log10CFU_500tid，'d.剂量500 tid'）链接子图轴linkaxes (hax1)
XLim: [0 1] YLim: [0 1] XScale: 'linear' YScale: 'linear' GridLineStyle: '-' Position: [0.1300 0.5838 0.3347 0.3412] unit: 'normalized'使用GET显示所有属性hax1 = 1×2 Axes array: Axes Axes hax1 = 1×3 Axes array: Axes Axes hax1 = 1×4 Axes array: Axes Axes Axes
肾功能对抗菌活性的影响
最后，作者比较了不同剂量方案对肾功能的影响。他们根据肌酐清除率将患者分为四个肾功能组(CrCL)：
肌酐清除组1:
CrCL< 30肌酐清除组2:30 <=
CrCL< 50肌酐清除率3组:50 <=
CrCL< 70肌酐清除组4:
CrCL> = 70
下图为肾功能(肌酐清除率)对四种给药方案抗菌效果的影响。观察肾功能正常组(CrCL>= 70)， 4种治疗策略的疗效差异显著。在这种情况下，500毫克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)属性:Number: 1名称:"颜色:[1,1]位置:[440 378 560 420]单位:'pixels'使用GET显示所有属性
您也可以从以下列表中选择网站: