PK / PD建模和仿真指导抗生素给药策略
这个例子演示了如何对抗菌剂的药代动力学/药效学(PK/PD)模型进行蒙特卡洛模拟。这个例子还展示了如何使用SimBiology®SimFunction对象并行执行参数扫描。
此示例需要统计和机器学习工具箱™。如果您有并行计算工具箱™软件，则可以提高性能。
背景
Katsube等人[1]使用PK/PD建模和模拟方法来确定碳青霉烯类抗生素多力培南的最有效剂量方案。他们的研究目标是:
建立多力培南对铜绿假单胞菌的PK/PD模型
使用Monte Carlo模拟来比较四种常见的抗生素剂量制度的疗效，并确定最有效的给药策略
观察肾功能对抗菌治疗效果的影响
在本例中，我们将在SimBiology®中实现Katsube等人[1]开发的抗菌PK/PD模型，并复制他们工作中描述的蒙特卡洛模拟结果。
参考文献
Katsube、Yano、T. Wajima、Y. Yamano和M. Takano。确定多力培南有效剂量方案的药代动力学/药效学建模和模拟。制药科学杂志(2010) 99(5), 2483 - 91。
PK / PD模型
Katube等人。假设具有从中央隔间的线性消除的双隔室输液模型，以描述多翅膜的药代动力学。对于细菌生长模型，他们认为总细菌群由药物易感生长细胞和药物不敏感静息细胞组成。通过简单的Emax型模型将药物的抗菌作用包含在细菌的杀伤率中：
![$$ kill Rate = {Kmax*[Drug]*[Growing]}{KC50 + [Drug]} $$](http://www.tatmou.com/jp/help/examples/simbio/win64/antibioticdemo_eq06411219986050949670.png)
在哪里[药品]药物在中央腔室的浓度(ug/ml)和(增长)为CFU/ml中不断增长的细菌数量(CFU =菌落形成单位)。Kmax是最大杀伤率常数（1 /小时）和KC50.为Michaelis-Menten速率常数(ug/ml)。
模型的SimBiology实现的图形化视图如下所示。
％加载模型sbioloadproject（“AntibacterialPKPD.sbproj”，'m1'） ;
剂量方案
Katsube等人使用四种常见的抗生素剂量策略模拟了该模型。
每日两次250毫克(每日一次)
每天250毫克三次（T.I.D.）
500毫克，一天两次(每日两次)
500毫克，一天3次(t.i.d)
四种给药方案均采用输注给药，输注时间设置为30分钟。在SimBiology中，这些剂量方案已经作为剂量对象实施。
%在模型中选择剂量对象doseNames = {“250毫克报价”，“250毫克tid”，“500毫克报价”，“500毫克tid”};为iDoseGrp = 1:length(doseNames)“名字”, 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²/√(v+m²))= @(m,v)√(log(v/m²+1))m = @(typicalValue) typicalValue;v = @(typicalValue,CV) typicalValue^2*CV^2;％患者人口统计学rng (“默认”）;Wt = normmrnd (51.6, 11.8, nPatients, 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);％单位：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/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 = logrnd (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 /分钟
药效学（PD）参数的分布：
生长-静止转换速率常数，K1和K2，用典型值为5.59E-5和0.02971 /小时的典型值取样，每次具有20％的CV。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 = logrnd (k1_mu, k1_sigma, nPatients, nDoseGrps);％单位：1 /小时k2 = lognrnd（k2_mu，k2_sigma，npatients，ndosegrps）;％单位：1 /小时Kmax = logrnd (Kmax_mu, Kmax_sigma, nPatients, nDoseGrps);％单位：1 /小时
Katsube等人假设了这个值K1，K2和Kmax与被处理的菌株无关。的价值bet，净生长速率常数，固定在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];%使用randsample从离散分布中取样MIC值MIC = 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'};
选择反应。
可见= {'[细菌生长模型] .gring'，．..“细菌生长模型。休息的};
设定一个模板剂量。
tempdose = sbiodose (“剂量”）;tempdose。目标='central.drug'；tempdose。AmountUnits =毫克的；tempdose.timeUnits =.“小时”；tempdose。DurationParameterName =“TDose”；
创建一个SimFunction对象。集UseParallel为the启用并行计算。
simfunc = createImfunction（m1，params，观察到，tempdose，“UseParallel”,真正的);
创建一个输入矩阵披每一剂量组。
φ=细胞(1、nDoseGrps);为i = 1: nDoseGrpsφ{我}=[中央(:,i)、k12(:,我),k21(:,我)．..CL(:,我),k1 (:, i), k2(:,我)．..kmax（：，i），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是细菌总数。每个剂量方案的有效性是通过在该剂量组中达到成功标准的人口比例来衡量的。这个功效指标,公关{log10 (CFU) < 0}，以时间为函数跟踪每个剂量组。
在他们的模拟研究中，作者研究了剂量方案对两类患者的疗效：
中度感染(初始细菌计数= 1e4 CFU/ml)
严重感染（初始细菌计数= 1E7 CFU / mL）
在此示例中，我们将仅对严重的感染案件复制结果。注意，通过改变细菌数量的初始量（该方法，您可以轻松模拟其他情况，中度感染的患者（日益增长的种)，在模型中为1e4 CFU/ml。
％preallocatecfu =南(nTPoints nPatients);log10CFU = cell(1,nDoseGrps);为i = 1:nDoseGrps disp([“模拟组”，num2str（i），'...']）%直接从现有剂量对象获取剂量表％给药方案。doseTable =可以获得的(doseRegimens (i));%模拟simdata = simfunc（phi {i}，[]，装饰，tobs）;%每名患者生长和静止细菌计数之和为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）=子图（2,2,1）plotcfucount（tobs，log10cfu_250bid，'一种。剂量250出价') hax1(2) = subplot(2,2,2)'b.剂量250 tid') hax1(3) = subplot(2,2,3)“c. 500剂量出价”) hax1(4) = subplot(2,2,4)d. 500剂量tid）链接子图轴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 mg t.i.d.剂量比其他方案更有效。相比之下，涉及肾功能障碍患者的模拟（CrCL< 30和30 <=CrCL< 50)，我们没有看到治疗组之间的太大差异。这表明，对于肾功能不全的患者，低强度或低频率的给药策略与高频率或高剂量的给药策略几乎一样有效。
％preallocateidCrCLGrp = false(nPatients, nDoseGrps);％线条风格ls = {'BD：'，" 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));%的阴谋图(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。”})传说位置西北传奇boxofflinkaxes (hax2)
f =图（1）具有属性：Number：1名称：''颜色：[1 1 1]位置：[680 678 560 420]单位：'像素'使用可以显示所有属性
您还可以从以下列表中选择一个网站：
