中药小分子模拟扩散系数中模拟时间和并行模拟问题_吴平

84
15 卷
辽
宁
中
医
药
大
学
学
报
1
计算方法
模拟体系分为 POPC 膜和镶嵌于其中的碳酸 钙 （CaCO3） 分子两个部分。磷脂膜式由上下两层 各 119 个 POPC 组成， 每个 POPC 分子具有相同的初 构 始结构。再将结构中插入 6 个碳酸钙 （CaCO3） 成的非晶态固体体系， 分子动力学模拟， MSD 分析 300K， NPT 都在 GROMACS[5]。主要参数设置如下 ： 系综， 压力设定为 105 Pa， 采用周期性边界条件， 时 间步长 0.001 ps， 每 10 fs 作结构取样。对范德华力 和库仑力均采用组截断方法 （group-based） ， 截断 半径为 8.5 Å。 本文连续作多次分子动力学模拟， 每次 300 ps。 后一次模拟均是采用上一次模拟的最终结构为起始 结构。由于数据点过于密集， 为了显示效果， 以下各 图 （图 2~4） 是在原始数据的基础上每 500 个数据留 下一个数据作图所得。 体 系 出 现 正 常 的 Fick 扩 散 的 时 候 ： 满 足 log （MSD） =αlog（t） ， （α=0.9~1.1） ， 从进入该时刻后， 记录模拟得到的轨迹 （即获得 MSD 随时间变化的数 据） ， 然后根据公式 （1） 即可求得所需的扩散系数。 2 1 d N D lim ri t r0 （1） t 6N dt i 1 其中距离均方 MSD
85
辽
宁
中
医
药
大
学
学
报
15 卷
图5
曲线 B-S 对应扩散系数比较示意图
给体系中各个原子一定的初始速度。因此初速度向 量对模拟的结果有着直接并且严重的影响。因而对 相同模型在相同的模拟参数下的 “并行模拟” ， 会有 不同的结果。 为 此， 我们在拷贝最终结构时设法保留速度 数据。具体做法是 ： 将模型按常规方法构建好， 进 行短时间的模拟 （原则上只要能产生轨迹文件就可 以） ， 然后将模拟停止， 保存好轨迹文件。在接下来 的并行模拟中， 导入刚才保存好的轨迹文件， 截取最 后一个 Frame（原则上任何时刻的一个都行） ， 即获 取带有速度的结构。 我们任意取一个轨迹文件 （如模拟编号 H） ， 截 取最后的一个 Frame， 然后作编号 I（1） 的模拟。同 样， 又将此轨迹文件再次截取最后的一个 Frame， 作 编号 I（2） 的模拟。结果见图 8。曲线 I（1） 和 I（2） 完全重合。我们又任意取一个轨迹文件 （如模拟编 号 K） ， 作同上述同样的处理， 即作保留速度数据的 并行模拟， 称为 L（1） 和 L（2） 。结果见图 9。曲线 L（1） 和 L（2） 完全重合。 将图 7 和图 8、 9 比较， 可见前者的 “并行模拟” 结果不同， 其原因是没有保留速度数据、 只有坐标数 据造成的。所以在接下来的任何并行模拟 （例如考 察温度、 压强等对该体系扩散系数的影响） 中一定 要注意数据的完整。
Simulation Time and Parallel Method in Calculate Diffusion Coefficient during TCM Simulation WU Qiong1， ZHU Ruixin2， KANG Tingguo1 WU Ping1， （1.College of Pharmacy， Liaoning University of Traditional Chinese Medicine， Dalian 116600， Liaoning， China ； 2.School of Life Sciences and Technology， Tongji University， Shanghai 200092， China） In this study， we employed molecular dynamics simulation to research the diffusion coefficient Abstract ： we of CaCO3-POPC system.CaCO3 is the main active ingredients of Abalone shell in TCM.During research， found the diffusion coefficient of simulation system fluctuated within a narrow range over a period of time. Besides， irregular simulation-curve is always here.Our article also gives a reasonable parallel method for simulation.This study lays the foundation for the diffusion simulation of TCM in biological systems. TCM small molecular ； molecular dynamics simulation ； diffusion ； simulation time ； Key words ： parrarel method 近年来随着中医药现代化的进程不断推进， 关 于中药小分子活性机理的研究受到广泛关注。科技 发展日新月异， 计算机性能逐步提高， 使得运用分子 动力学方法模拟中药小分子的疏运性质成为可能， 其中扩散系数就是一种重要的运输性质。 小分子在无定型聚合物膜中的扩散是一类重 要的扩散现象， 它在膜分离、 食品包装、 敏感器的应 用等方面具有广泛应用。实验技术不能在分子水平 上给出详细的描述， 在宏观水平上获得的信息不能 够解释受化学基团的类型、 化学链和基团的空间关 系， 以及它们的热力学波动等因素影响的气体运动 它 机理， 而分子动力学模拟 [1-2] 恰能弥补这一缺陷， 既能避免昂贵的实验， 又能从分子水平设计改良聚 合物膜的渗透和选择性等问题。到目前为止， 运用 的分子动力学模拟并且根据 Einstein-Smoluchowski 理论计算扩散系数一直受到模拟时间的影响， 并且 相同条件下的模拟结果也不同， 也就是无法对相同 体系进行并行的研究。 中药石决明为鲍科动物杂色鲍、 皱纹盘鲍、 耳 鲍、 羊鲍等的贝壳。本品为常用中药， 《名医别录》 列 为上品。有平肝清热、 明目去翳的功效。石决明的 药用价值在古代就被人们所熟悉， 在 《雷公炮炙论》 《本草经集注》 《日华子本草》 《本草纲目》 等古代医 药书籍中均有记载。石决明的有效成分为碳酸钙。 在 本文以中药石决明有效成分碳酸钙 （CaCO3） 中的扩散为模拟体系， 对目 细胞膜 POPC[3-4]（图 1） 前模拟中存在模拟时间和模拟结果重复 （或并行） 问题进行了研究。POPC 全名是棕榈酰油酰磷脂酰 胆碱， 是组成细胞膜的磷脂双分子层的一种。另外 本文利用本研究组提出的变分原理， 对模拟得到的 MSD 曲线进行了处理， 得到了精确可靠的扩散系 数。
在实际模拟的过程中发现， 即使对于早已进入 符合爱因斯坦方程的正常扩散过程， 扩散系数在不 同的模拟时间中依然存在震荡 （图 2果的单独 比较示意图 （从曲线 L 到 S）
图 2 以时间顺序每 300 ps 对上述体系进行 模拟的结果示意图 （时间顺序为从 B 到 S， 共 16 次模拟， 每次 300 ps， 总计 4.800 ns）
σ2 ▲ 0.15 0.16 σ2 ◆ 0.03 0.07 D ▲ 2.12 2.54 D ◆ 0.45 0.84
从各曲线求得的扩散系数 D（cm2/s×10-4）
B C D E F G H I K L N O P Q R S D▲ D◆ 2.40 1.87 0.73 0.60 0.52 0.45 0.28 0.44 0.68 0.51 0.52 0.49 0.33 0.42 0.31 0.40 0.68 0.48 2.73 1.43 1.03 0.57 0.39 0.50 0.19 0.21 0.91 0.43 0.34 0.49 0.24 0.31 0.24 0.60 0.66 0.46 注： a： 差分－区限变分法 ； b： 常规处理方法 ； ▲： 数据 B-S ； ◆： 数据 D-S。
图1
POPC 单体分子结构图
收稿日期 ： 2013-01-26 基金项目 ： 国家自然科学基金项目 （30973611） ； 中国高等教育博士研究基金计划项目 （20100072120050） ； 上海中药现代化基金项目 （09dZ1972800） 作者简介 ： 吴平 （1987-） ， 男， 江苏泰州人， 硕士研究生， 研究方向 ： 计算机辅助药物设计与天然产物结构修饰。 通讯作者 ： 康廷国 （1955-） ， 男， 内蒙古赤峰人， 教授， 博士研究生导师， 研究方向 ： 中药鉴定与品质评价。
1
从图 2 可以得出如下结论 ： 前 2 次 300 ps 模 拟 （图中曲线 B 和 C） 结果比其余的曲线偏离相对较大， 这是由于非晶态固体相对于晶体的结构过于复杂， 系统处于一种亚稳态。因此认为该 2 次的模拟结果 尚未进入稳态。 从图 3 和 4 可以得出如下结论 ： 曲线的斜率并 没有按照模拟时间的先后呈现单调变化的趋势， 即
表1 编号 Da Db
n 1 i 1
Di
n
D
2
， 其中 D

1
n
Di ； △ D Dmaxi
1
n
Dmin。 图 5、 图 6 分别是表 1 数据所作的图， 图5和6说 明进入稳态后本例的扩散系数涨落还是不小的。本 文的研究再次证明 ： 如图 2 所示， MSD 曲线尾部异 常现象的存在是由于计算 MSD 的算法造成的， 而不 是模拟时间不够的问题 [12]。
对 于 扩 散 系 数 的 获 得， 文献中大多采用对 MSD 曲 线 直 接 求 斜 率 的 方 法 [6-8] 或 者 人 为 舍 弃 曲 但 是 由 于 存 在 MSD 计 算 误 线 不 规 则 的 部 分 [9-11]， [12] 为 此， 本 研 究 组 曾 经 根 据 Einstein差的原因 ， Smoluchowski 方程的物理意义， 提出了利用 “差分－ 区限变分法” 从各 MSD vs.t 曲线求算扩散系数 [13-14]。 将从曲线 B 到曲线 S 分别采用 “差分－区限变分法” 和常规方法处理， 得到如下的扩散系数。见表 1。 σ2 为每次 300 ps 模拟求得的扩散系数 Di 的涨 落， 即 2