PK/PD建模与仿真指导抗生素的给药策略
这个例子演示了如何对抗菌剂的药代动力学/药效学(PK/PD)模型进行蒙特卡洛模拟。这个例子还展示了如何使用SimBiology®SimFunction对象并行执行参数扫描。
这个例子需要Statistics和Machine Learning Toolbox™。如果您有并行计算工具箱™软件，则可以提高性能。
背景
Katsube等人[1]使用PK/PD建模和模拟方法来确定碳青霉烯类抗生素多力培南的最有效剂量方案。他们的研究目标是:
建立多力培南对铜绿假单胞菌的PK/PD模型
使用蒙特卡洛模拟来比较四种常见抗生素剂量方案的效果，并确定最有效的剂量策略
观察肾功能对抗菌治疗效果的影响
在本例中，我们将在SimBiology®中实现Katsube等人[1]开发的抗菌PK/PD模型，并复制他们工作中描述的蒙特卡洛模拟结果。
参考文献
Katsube、Yano、T. Wajima、Y. Yamano和M. Takano。确定多力培南有效剂量方案的药代动力学/药效学建模和模拟。制药科学杂志(2010) 99(5), 2483 - 91。
PK / PD模型
Katsube等人假设了一个由中央室线性消除的双室输注模型来描述多立培南的药代动力学。对于细菌生长模型，他们假设总细菌群由药物敏感的生长细胞和药物不敏感的静息细胞组成。通过简单的Emax模型，将药物的抗菌作用纳入对细菌的杀灭率:
![$$ kill Rate = {Kmax*[Drug]*[Growing]}{KC50 + [Drug]} $$](http://www.tatmou.com/uk/help/examples/simbio/win64/antibioticdemo_eq06411219986050949670.png)
在哪里(药物)药物在中央腔室的浓度(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毫克，一天两次(每日两次)
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);%单位:毫升/分钟性别比例%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 = logrnd (k12_mu, k12_sigma, nPatients, nDoseGrps);%单位:1 /小时k21 = logrnd (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 = 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等人假设了这个值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];%使用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”};
选择反应。
可见= {“细菌生长模型。日益增长的，．..“细菌生长模型。休息的};
设定一个模板剂量。
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)];结束
集群计算
这个示例使用预先配置到本地机器的本地集群配置文件。您还可以搜索其他在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)
在本例中，我们只复制严重感染病例的结果。注意，你可以很容易地模拟另一种情况，中度感染的患者，通过改变初始的细菌计数(the日益增长的种)，在模型中为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)“a. 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毫克每分钟的剂量比其他方案更有效。与此相反，模拟涉及肾功能不全患者(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));%的阴谋图(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 = Figure (1) with properties: Number: 1 Name: " Color: [1 1 1] Position: [680 678 560 420] Units: 'pixels'使用GET显示所有属性
你也可以从以下列表中选择一个网站:
