一种基于蒙特卡罗模拟的抗菌药物药效计算方法与流程

本发明属于体外药动学/药效学(Pharmacokinetic/pharmacodynamic, PK/PD)研究领域，涉及一种基于蒙特卡罗模拟（Monte Carlo Simulations, MCS）的抗菌药物药效计算方法。

背景技术：

蒙特卡罗模拟，或称计算机随机模拟方法，是一种利用随机数进行数值模拟的方法。蒙特卡罗的基本思想很早以前就被人们所发现和利用。17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。20世纪40年代中期随着科学技术的发展和电子计算机的发明，为了适应当时原子能事业的发展，有人提出了蒙特卡罗模拟,，这是一种以概率统计理论为指导的一类非常重要的数值计算方法。MCS是采用不同统计取样技术（如随机数字，伪随机数字等）来提供定量问题近似解决方案的随机模拟方法。蒙特卡罗模拟的基本思想是：当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征,并将其作为问题的解。简单来说，MCS是基于人工创造一个随机事件或“实验”（通常是通过计算机）的分析方法，当运行到指定的次数，可以获得任何特定目标的概率。

药动学（pharmacokinetics, PK）和药效学（pharmacodynamics, PD）是临床药理学的两个重要组成部分。PK是研究机体对药物的吸收、分布、代谢和排泄规律，主要反映药物血药浓度随时间变化的规律。PD是研究血药浓度和药物效应之间的关系。根据抗生素PK和PD相关性，可将其分为浓度依赖型和时间依赖型两类。当药物的药时曲线下面积（AUC）或药物的峰浓度（Cmax）对药物杀菌能力有显著影响时，该类抗生素称为浓度依赖型抗菌药物，如氨基糖苷类、氟喹诺酮类、酮内酯类、两性霉素B和达托霉素等。此类药物的杀菌作用与合理的AUC与M IC的比值（AUC/M IC）或Cmax与M IC的比值（Cmax/M IC）相关。因此，此类药物的给药目标是最大化地提高药物的AUC或者峰浓度来提高临床疗效，但剂量不能超过最低中毒剂量。第二类是指当药物浓度是MIC的2～4倍时，其杀菌速率与药物浓度超过MIC的时间（%T>M IC）密切相关，此类药物亦称为时间依赖型抗生素，如β-内酰胺类和林可霉素等。目前用于评价时间依赖型药物杀菌作用的PK/PD参数主要是%T>MIC，此类药物的给药目标是尽量延长药物浓度超过MIC的时间，而不出现药物不良反应。

过去，人们主要通过AUC的平均值以及MIC90来衡量药物在一个固定的给药剂量下的药效；这种评估方式忽略了药物在体内的个体差异以及细菌的MIC值的广泛分布使得简单的通过比较AUC/MIC的某个值的分析方法很难全面地评估给药后药物的在动物体内真实的药效。

技术实现要素：

本发明所要解决的技术问题是，提供一种同时能预测浓度依赖型和时间依赖型抗菌药物药效，并能给出给药方案的计算方法。

本发明所采用的技术方案是：一种基于蒙特卡罗模拟的抗菌药物药效计算方法，包括如下步骤：

1）确定抗菌药物属于浓度依赖型还是时间依赖型，选择最能预测该药物体内活性的PK/PD参数（AUC/MIC或%T>MIC）；以及确定体内产生杀菌活性且不会出现耐药性的PK/PD参数的目标值；

2）若药物属于浓度依赖型药物，PK/PD参数为AUC/MIC，AUC=Dose/CL，AUC/MIC=Dose /CL/MIC，其中，Dose是24h给药剂量（mg/kg b.w.），CL为体清除率（L/h/kg），MIC是最低抑菌浓度(μg/mL)；

3）若该药物属于时间依赖型药物，PK/PD参数为%T>MIC，对于静脉注射一室模型、一级吸收一室模型、静脉注射二室模型或一级吸收二室模型，先采用二分法求出浓度方程f(c,θ,t)-MIC=0的解t1和t2，再通过公式(t2-t1)/24求出%T>MIC的值；

4）对参数初始化，设置蒙特卡罗模拟运行次数（通常5 000次或10 000次）；

5）对动学参数和细菌MIC的分布进行假设，通常假设药动学参数服从对数正态分布(Log Gaussian distribution)，MIC根据统计学的结果服从自定义分布(custom distribution)；

6）运行模拟，对药效指标AUC/MIC或者%T>MIC是否达标进行判断，并记录结果；

7）重复模拟，记录每次循环指标的判断结果，直到循环终止；

8) 根据每次的判断结果，计算达标率TAR，并输出其它参数，以及对结果作图表。

所述步骤3）中，静脉注射一室模型的浓度方程为：

；

一级吸收一室模型的浓度方程为：

；

静脉注射二室模型的浓度方程为：

；

一级吸收二室模型的浓度方程为：

。

本发明所提供的基于蒙特卡罗模拟的浓度依赖型抗菌药物药效计算方法，其特征在于，包括如下步骤：

1）根据药物A的杀菌性能，确定用于蒙特卡罗分析的AUC/MIC目标值T_A，根据A药物的药动学特征，确定其体清除率CL值为，对临床分离细菌的最小抑菌浓度MIC的统计分析，获得细菌MIC的概率分布；

2）确定初始参数，确定给药剂量Dose，定义循环次数i为N（如10000）次，初始值为1，循环增加1次，i=i+1，k为AUC/MIC>= T_A的次数，初始值为0，若第i次循环时AUC/MIC[i]>=T，k=k+1，j为AUC/MIC< T_A的次数，初始值为0，若第i次循环时AUC/MIC[i] < T_A，j=j+1，根据细菌临床分离株对药物A的MIC的分布概率，定义在[0,1]均匀分布的取值区间对应的细菌MIC值；

3）在Matlab软件中对CL进行第一次模拟，得到一个CL[1]的值，根据公式AUC=Dose/CL，得到一个AUC[1]的值；

4）在[0,1]区间中按均匀分布规则，随机抽取一个值Random[1]，判断Random[1]所在区间，赋予一个MIC[1]的值；

5）根据公式AUC/MIC[1] = AUC[1]/MIC[1]，求出AUC/MIC[1]的值，若AUC/MIC[1] >=T_A，k=k+1，反之，j=j+1；

6）对步骤3）~5）进行重复模拟，记录每一个参数，如CL[i]、AUC[i] 、 MIC[i], 、AUC/MIC[i]的值，直到i=N次，循环终止；

7）根据公式TAR=k/i，得到TAR的结果，并输出其它参数，以及对CL、AUC、MIC、AUC/MIC的结果作图表。

本发明所提供的基于蒙特卡罗模拟的时间依赖型抗菌药物药效计算方法，其特征在于，包括如下步骤：

1）根据药物B的杀菌性能，确定用于蒙特卡罗分析的%T>MIC的目标值为T_B；对于B药物的药动学模型，采用一级吸收二室模型，其公式为，A、B、α、β为混杂参数，k_a为吸收速率常数，并确定各参数的值，对临床分离细菌的最小抑菌浓度MIC的统计分析，获得细菌MIC的概率分布；

2）确定初始参数，确定给药剂量Dose，定义循环次数i为N（如10000）次，初始值为1，循环增加1次，i=i+1，k为[(T>MIC)/24]>= T_B的次数，初始值为0，若第i次循环时[(T>MIC)/24]>= T_B，k=k+1，j为[(T>MIC)/24]< T_B的次数，初始值为0，若第i次循环时[(T>MIC)/24]< T_B，j=j+1，根据细菌临床分离株对药物A的MIC的分布概率，定义在[0,1]均匀分布的取值区间对应的细菌MIC值；

3）在Matlab软件中对A、B、α、β、ka进行第一次模拟，分别得到一个A[1]、B[1]、α[1]、β[1]、ka[1]的值，根据公式，可得

；

4）在[0,1]区间中按均匀分布规则，随机抽取一个值Random[1]，判断Random[1]所在区间，赋予一个MIC[1]的值；

5）采用二分法求非线性方程C[1]-MIC=0，即

的两个解t1和t2；

6）利用公式%T>MIC=(t2-t1)/24，求出%T>MIC [1]的值，若%T>MIC [1]>=0.35，k=k+1，反之，j=j+1；

7）对步骤3）~6）进行重复模拟，记录每一个参数，如A[i]、B[i]、α[i]、β[i]、k_a[i]、MIC[i]、%T>MIC[i]的值，直到i=N次，循环终止；

8）根据公式TAR=k/i，得到TAR的结果并输出其它参数，以及对%T>MIC的统计学参数及概率统计等相关结果作图表。

本发明的有益效果为：本发明所提供的抗菌药物药效计算方法，是以蒙特卡罗模拟方法为基础，蒙特卡罗模拟( Monte Carlo simulation, MCS) 为抗菌药物经验治疗和个体化治疗提供了发现高效的抗菌治疗方案的方法，通常表达为对细菌群体MIC的目标获得概率（target attainment rate, TAR），能达到最高TAR或能达到90%TAR的给药方案可能是提高抗菌药物经验治疗的合理选择。本发明能同时预测浓度依赖型和时间依赖型抗菌药物药效，以便能提供合理的给药方案。

附图说明

图1为实施例1蒙特卡罗模拟得到的CL分布图；

图2为实施例1蒙特卡罗分析得到的AUC分布图；

图3为实施例1 MIC的统计值（左列）和蒙特卡罗模拟值（右列）；

图4为实施例1蒙特卡罗模拟得到的AUC/MIC分布图；

图5为实施例1剂量为10mg/kg b.w. （body weight，体重）的药物对不同MIC值的细菌的达标率（TAR）。

具体实施方式

以下结合实施例对本发明作进一步的说明，但并不局限于此，本领域技术人员在不脱离本发明的基本技术思想前提下，所作出的修改、替换、变更均落入本发明的保护范围。

实施例1

对浓度依赖型药物A，具体步骤为：

（1）根据药物A的杀菌性能，确定用于蒙特卡罗分析的AUC/MIC目标值120 h；根据A药物的药动学特征，确定其体清除率CL值为0.421±0.058 L/h/kg（）；对临床分离细菌的最小抑菌浓度MIC的统计分析，获得细菌MIC的概率分布（见表1）。

细菌MIC的概率分布如表3-1所示。

表1 细菌临床分离株对药物A的MIC的分布概率

（2）确定初始参数

给药剂量Dose=10mg/kg。定义循环次数i为10000次，初始值为1，循环增加1次，i=i+1。k为AUC/MIC>=120h的次数，初始值为0，若第i次循环时AUC/MIC[i]>=120h，k=k+1，j为AUC/MIC<120h的次数，初始值为0，若第i次循环时AUC/MIC[i] <120h，j=j+1。根据细菌临床分离株对药物A的MIC的分布概率，定义在[0,1]均匀分布的取值区间对应的细菌MIC值（表2）。

表2 均匀分布的取值区间对应的MIC值

（3）在Matlab软件中对CL进行第一次模拟，得到一个CL[1]的值，根据公式AUC=Dose/CL，得到一个AUC[1]的值

（4）在[0,1]区间中按均匀分布规则，随机抽取一个值Random[1]，判断Random[1]所在区间，赋予一个MIC[1]的值 。

（5）根据公式AUC/MIC[1] = AUC[1]/MIC[1]，求出AUC/MIC[1]的值。若AUC/MIC[1] >=120h，k=k+1，反之，j=j+1。

（6）对步骤3~5进行重复模拟，记录每一个参数，如CL[i], AUC[i], MIC[i], AUC/MIC[i]的值，直到i=10000次，循环终止。

（7）根据公式TAR=k/i，得到TAR的结果为77.3%，即10 mg/kg b.w. （body weight，体重）的给药剂量在临床的治疗中有77.3%的生物体的AUC/MIC值>=120h，达到较好的治疗效果。运行模拟后，得到的AUC/MIC的统计学参数如表3所示。 对MIC值不同的细菌，AUC/MIC大于30、40、50、60、70、80、90、100、110、120的概率统计如表4所示。

表3 描述AUC/MIC分布的统计参数

表4 细菌在各个MIC范围内达到特定的AUC/MIC比值的概率结果

（8）经过10000次从各自的概率密度函数中随机抽样后，得到的CL、AUC、MIC、AUC/MIC的分布图分别如图1~4所示。图5为不同MIC值所对应的AUC/MIC的达标率分布图（用箱形图表示），AUC/MIC>=120h为达标。从图中可以看到，当MIC为0.125μg/mL或以下时，药物对细菌的杀菌效果非常好，且0.008~0.0625的达标率都是100%，随着细菌MIC的增大，达标率也相应的变小。

实施例2

对于时间依赖型药物，具体步骤为：

（1）根据药物B的杀菌性能，确定用于蒙特卡罗分析的%T>MIC的目标值为35%；对于B药物的药动学模型，本实施例采用较为复杂的一级吸收二室模型，其公式为

，A、B、α、β为混杂参数，k_a为吸收速率常数。其中各参数的分别为：A=10.40±2.64，B=7.39±2.22，α=1.98±0.81，β=0.39±0.03，k_a=5.78±0.96。对临床分离细菌的最小抑菌浓度MIC的统计分析，获得细菌MIC的概率分布（见表5）。

表5 细菌临床分离株对药物B的MIC的分布概率

（2）确定初始参数

给药剂量Dose=5mg/kg。定义循环次数i为10000次，初始值为1，循环增加1次，i=i+1。k为[(T>MIC)/24]>=0.35的次数，初始值为0，若第i次循环时[(T>MIC)/24]>=0.35，k=k+1，j为[(T>MIC)/24]<0.35的次数，初始值为0，若第i次循环时[(T>MIC)/24]<0.35，j=j+1。根据细菌临床分离株对药物A的MIC的分布概率，定义在[0,1]均匀分布的取值区间对应的细菌MIC值（表6）。

表6 均匀分布的取值区间对应的MIC值

（3）在Matlab软件中对A、B、α、β、k_a进行第一次模拟，分别得到一个A[1]、B[1]、α[1]、β[1]、k_a[1]的值，根据公式，可得

。

（4）在[0,1]区间中按均匀分布规则，随机抽取一个值Random[1]，判断Random[1]所在区间，赋予一个MIC[1]的值 。

（5）采用二分法求非线性方程C[1]-MIC=0，即

的两个解t1和t2。具体方法包括：取初始区间[a, b]，当f(a)*f(b)<0时，取中点x= (a+b)/2代替a或者b，当区间长度小于1e-6时，停止计算，得到t1和t2两个近似解。

（6）利用公式%T>MIC=(t2-t1)/24，求出%T>MIC [1]的值，若%T>MIC [1]>=0.35，k=k+1，反之，j=j+1。

（7）对步骤3~6进行重复模拟，记录每一个参数，如A[i]、B[i]、α[i]、β[i]、k_a[i]、MIC[i]、%T>MIC[i]的值，直到i=10000次，循环终止。

（8）根据公式TAR=k/i，得到TAR的结果为84.1%，即5 mg/kg b.w.（body weight，体重）的给药剂量在临床的治疗中有84.1%的生物体的%T>MIC值>=0.35，达到较好的治疗效果。运行模拟后，得到的%T>MIC的统计学参数如表7所示。 对MIC值不同的细菌，%T>MIC大于0.3、0.35、0.4、0.5、0.6、0.7、0.8、0.9的概率统计如表8所示。

表7 描述%T>MIC分布的统计参数

表8 细菌在各个MIC范围内达到特定的AUC/MIC比值的概率结果

上述实施例和说明书中描述的只是说明本发明的原理和最佳实施例，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。