聚合物前药载药性能的计算机模拟

罗雪莉， 郎美东

（上海市先进聚合物材料重点实验室，华东理工大学材料科学与工程学院，上海 200237）

聚合物前药的研究始于20 世纪50 年代中期，人们将某些药物与天然聚合物或合成聚合物以化学键的形式相连接形成聚合物前药，聚合物前药可以减少药物的肾排泄量，增强药物活性持续时间[1]。聚合物前药模型首先由Ringsdorf 提出并设计，该模型包括聚合物主链、药物、间隔基、靶向基团和增溶剂[2]。除了前药外，胶束也是常用的药物载体[3]，但是与前药以化学键连接不同，胶束是一种物理包封方法，药物通常进入胶束的疏水内核，因此其载药效率通常与药物和聚合物疏水段的相互作用力有关。

Flory-Huggins 相互作用参数(χ)是高分子物理中重要的参数之一，它反映了高分子与溶剂或者高分子与高分子混合时相互作用能的变化[4]。在药物载体领域，Flory-Huggins 相互作用参数可用于计算药物和载体之间的相容性，从而预测药物载体的载药性能。Flory-Huggins 相互作用参数的计算方法有实验法[5]、溶度参数法[6-7]、混合能法[8-9]、结合能法[10]等。其中溶度参数法只适用于没有氢键和极性基团的体系，混合能法通常通过计算机模拟方法得到，该方法适用于含有氢键和极性基团的体系，但是在计算中要求体系达到的平衡度高，计算时间长且计算过程复杂。

计算机模拟可以通过理论计算预测实验结果，已经被证实是一种强有力的科学研究手段[11]，被广泛应用于各领域的研究，在药物载体方面也发挥着重要的作用。

文献[12]运用耗散粒子动力学(DPD)模拟研究了pH 敏感的两亲性三嵌段共聚物聚(ε-己内酯)-b-聚(甲基丙烯酸二乙氨基乙酯)-b-聚(磺基甜菜碱甲基丙烯酸酯)/聚(乙二醇甲基丙烯酸酯)(PCL-b-PDEA-b-PSBMA/PEGMA)的微观自组装结构以及该聚合物加载和释放药物阿霉素的性质，发现随着聚合物摩尔浓度的不同，聚合物自组装形态呈现球形胶束、圆柱形结构、层状胶束等不同的状态，且胶束具有pH 响应性的药物释放性质。Gupt 等[13]基于分子动力学(MD)模拟，计算内聚能密度和溶度参数，预测了吲哚美辛在不同载体(聚环氧乙烷、葡萄糖、蔗糖)中的相容性，并用热分析的实验手段验证了模拟结果，表明了MD 模拟是一种确定物质溶度参数和预测药物与载体混溶性的有效技术。

二代卟啉类光敏剂5-氨基酮戊酸( 5-ALA)是内源性原卟啉IX(PpIX)的代谢前体物质，近年来在临床上被广泛应用于光动力治疗[14]。5-ALA 具有不良反应小、疗效确切等优点，但是其亲水性非常强，因此在局部应用时透皮效果不理想。对5-ALA 进行酯化可改善药物透皮性，获得性质稳定、靶组织选择性高、光动力效应强的衍生药物。5-氨基酮戊酸甲酯就是一种5-ALA 的酯化衍生药物，其商品名为Metvix，可用于光化性角化病和表面结节型基地细胞癌的治疗[15]。除了可通过酯化的方式改善5-ALA 的透皮效果外，胶束载药的方式也能改善其透皮性。

5-ALA 亲水性很强，通过普通载体进行物理包埋时包载效率很低，通过化学键接的方式时载药量有限，需要大量的活性基团。本文用5-ALA 修饰聚乙二醇单甲醚-b-聚(ε-己内酯-co-哌啶内酯)(mPEG-b-P(CL-co-PIL))得到聚合物前药载体聚乙二醇单甲醚-b-聚(ε-己内酯-co-哌啶内酯-g-5-氨基酮戊酸)(mPEGb-P(CL-co-APIL))，并以该前药载体包载5-ALA，该体系同时采用化学键接和物理包埋的方法协同载药，且通过药物接枝的方式进一步改善载体与药物之间的相容性，提高载体对药物的物理包埋效率。以上述体系为模型，采用MD 模拟的方法计算不同药物接枝率的聚合物前药载体中疏水链段与药物的Flory-Huggins 相互作用参数(χij)，预测载体的载药性能，同时采用DPD 模拟研究体系介观形态，通过粒子聚集程度来预测载药效率，从而为载体前期设计提供理论指导，减少实验探索的时间。本文提出的以聚合物前药作为载体包载药物的模型可提高载体对某些包载较困难的药物的包载效率，为药物载体的设计提供了一种新的思路。

1 模拟方法

1.1 分子模型及计算体系的构建

采用Materials studio 6.1(MS，Accelrys)构建聚合物两种前药载体聚(ε-己内酯-co-哌啶内酯-g-5-氨基酮戊酸)P(CL-co-APIL)和聚(ε-己内酯-co-聚羧基化己内酯-g-5-氨基酮戊酸甲酯)(P(CL-co-mACCL))。聚羧基化己内酯(PCCL)、聚己内酯(PCL)、聚哌啶内酯(PPIL)、聚哌啶内酯-g-5-氨基酮戊酸（PAPIL）、聚羧基化己内酯-g-5-氨基酮戊酸甲酯（PmACCL)分子结构见图1。

图 1 聚合物的分子结构Fig. 1 Molecular structure of polymers

为了研究不同聚合物前药与药物分子之间的相容性，采用MD 模拟方法计算聚合物与药物之间的χij。利用Materials Studio 中的Discovery 模块对建立好的分子模型进行能量最小化，设定最小化的最长时间为20 000 fs，再利用Forcite 模块中的Geometry optimi-zation 对分子模型进行结构最优化。利用Amorphous cell 模块建立聚合物前药、聚合物亲水段mPEG、药物分子5-ALA、溶剂水的纯体系无定形模型以及聚合物/药物、聚合物/H2O、聚合物/聚合物、药物/水等混合体系的无定形模型，纯体系的无定形模型格子大小保持在3 nm×3 nm×3 nm 左右，混合体系的无定形模型格子大小则控制在3 nm×3 nm×3 nm～4 nm×4 nm×4 nm之间，纯体系分子的物质的量和混合无定形体系分子的物质的量之比见表1，部分无定形模型示意见图2。以上所有步骤均在COMPASS 力场[16]下进行。

1.2 MD 模拟方法

在COMPASS 力场下对无定形结构进行能量最小化，上限时长为20 ps；然后，进行正则系综（NVT）动力学模拟，模拟时长为150 ps；最后进行等温等压（NPT）动力学模拟，模拟时长为1.5 ns。依据密度和势能是否收敛适当延长计算时间，每5 ps 取一个轨迹，取最后的30 个格子轨迹进行计算，所有的模拟过程中步长为1 fs，计算得到每个格子的密度、内聚能密度、范德华力内聚能密度和静电力内聚能密度，并进行χij的计算。

χij及混合能的计算如式（1）～式（2）所示[8]：

式中：V 为格子体积；Vr为参考体积；ΔEmix为混合能；Ecoh为内聚能；Ecoh/V 为内聚能密度(以ρcoh表示)；φi、φj为各成分的体积分数，满足φi+φj=1。

表 1 纯无定形体系中分子的物质的量和混合无定形体系中各分子的物质的量之比Table 1 Molecular molar numbers of the compounds in pure amorphous system and molar ratio of the different compounds in mixed amorphous systems

图 2 部分无定形模型的示意图Fig. 2 Schematic diagrams of partial amorphous models

Vr一般选择小分子或聚合物结构单元的摩尔体积。为了计算更加精准，选择混合体系的平均摩尔体积作为Vr值，其计算式如式(3)所示：

式中：ni、nj为混合体系中各成分的物质的量；Mi、Mj为混合体系中各成分的摩尔质量；ρij为混合体系的密度。

溶度参数是表征分子间相互作用强度特征的重要参数，可预测聚合物之间或聚合物与其他小分子之间的相容性，其概念是Hildebrand 等[17]在20 世纪中期提出，定义为物质内聚能密度的平方根，Hansen[18]在Hildebrand 溶度参数的理论基础上将液体的内聚能视为色散力、极性力和氢键这3 种分子间作用力的贡献之和，建立了一个三维溶解度参数体系：

式中：δd是色散力对溶度参数的贡献；δp是极性力的贡献；δh是氢键的贡献。

借鉴Hansen 的工作，本文将χij分解为范德华力（vdw）分量相互作用参数(χij-vdw)和静电力（e）分量相互作用参数(χij-e)，其计算式如式(5)～式(9)所示：

1.3 DPD 模拟方法

DPD 模拟是一种介观尺度的模拟方法，该方法基于粗粒化模型，而MD 模拟基于原子量级，计算过程需要关注原子细节，过程更加复杂，耗时更长，难以达到平衡。对于高分子等分子量较大的体系，粗粒化模型有利于减少计算时间，被大量用于高分子自组装等复杂相领域[19]。DPD 的粗粒化模型将几个原子划分为一个珠子，同一个分子中的珠子与珠子之间通过弹簧相连，珠子间的作用力分别用保守力、耗散力与随机力表示。保守力（）的计算如式(10)所示[20]：

式中：aij为珠子i 与珠子j 之间的排斥参数，rij为两个珠子i、j 之间的距离，为其向量。当体系密度等于3 时，aij满足式(11)：

mPEG-b-P(CL-co-APIL)、5-ALA 和H2O 的粗粒化模型如图3 所示。粗粒化的原则是保证每个粒子的质量接近，因此本文粗粒化过程如下：mPEG 的3 个结构单元作为一个珠子，用E 表示(橙色)，聚己内酯的一个结构单元为一个珠子，命名为C(绿色)；聚哌啶内酯接枝5-ALA 的一个结构单元分成两个珠子，命名为P(蓝色)；一个5-ALA 为一个珠子，命名为A(红色)；8 个水分子为一个珠子，命名为W(棕色)。

计算并设置好排斥参数后，对不同拓扑结构的聚合物前药进行DPD 模拟的参数设置。为了得到一个完整的胶束形态，设置模拟格子大小为20rc×20rc×20rc以及30rc×30rc×30rc，其中rc为模拟体系的截断半径，聚合物、药物与水的物质的量之比为8∶2∶90，采用周期性边界条件，体系密度选为3，设置弹簧常数为4.0，耗散力参数设为4.5，模拟步长为20 000 步。

2 结果与讨论

2.1 聚合物与药物及水分子之间的相互作用参数

PAPIL20分子动力学模拟过程平衡图见图4。从图4(a)可以看出，在16 000 fs 之后体系趋于平衡，势能不再变化，说明在本文所用的计算方法及格子体积中，选择20 ps 的时长进行能量最小化足以使系统达到平衡。PAPIL 纯体系在NPT 模拟过程中的密度变化见图4(b)，从图中可以看出，在模拟时长达到1 000 ps 以后，密度变化基本上收敛，说明本文所用的计算方法及格子体积中，选择NPT 模拟时长为1.5 ns 并根据密度收敛情况适当延长模拟时间的方法是合理的。

图 3 mPEG-b-P(CL-co-APIL)(a)，5-ALA(b)和H2O(c)的粗粒化模型Fig. 3 Coarse-grained models of mPEG-b-P(CL-co-APIL)(a), 5-ALA(b) and H2O(c)

图 4 PAPIL20 分子动力学模拟过程平衡图Fig. 4 Molecular dynamics simulation balance diagrams of PAPIL20

纯无定形体系的密度和内聚能密度的实验结果见表2。对比聚合物PCL20、PPIL20、PAPIL20、PCCL20、PmACCL20可以看出，PCL20的静电力内聚能密度最小，因为它的碳链上没有活性基团；PAPIL20的静电力内聚能密度最大，原因是在5-ALA 接枝到哌啶内酯上以后，分子间存在大量的极性基团，且有裸露的氨基，容易形成分子间氢键。PCCL20的范德华力内聚能密度最大，PmACCL20的最小，原因可能是后者分子量最大，分子体积最大。对比5-ALA 和m-ALA 可以看出，两者的范德华力内聚能密度相差不大，而静电力内聚能密度相差较大，导致总的内聚能密度相差较大。从两者的结构考虑很容易理解这个结果，因为m-ALA 上只有一个裸露的氨基，其羧基被羧基甲酯取代，因此5-ALA 相比m-ALA 有更多的极性基团且更容易产生分子间氢键[21]。

不同分子之间的Flory-Huggins 相互作用参数见表3，参数越小代表两种分子间的相容性越好。对比并分析表3 中5-ALA/mPEG、5-ALA/H2O 和m-ALA/mPEG、m-ALA/H2O 的数据，发现相比5-ALA 分别与mPEG 和H2O 的相容性，m-ALA 与两者的相容性都变差，且变差的主要原因是静电力相互作用参数变大，该结果同样与极性基团的减少有关。

对比PCL20、PPIL20、PAPIL20、PCCL20、PmACCL205 种聚合物与5-ALA 的相互作用参数发现：PPIL20与PAPIL20相比（接枝5-ALA 前后），聚合物与5-ALA 的相互作用参数从−0.84 下降为−2.65，其中范德华力相互作用参数从−0.70 下降为−2.59，静电力相互作用参数从−0.14 略有上升至−0.06，说明其相容性的改变主要依赖于范德华力的变化，而范德华力的变化主要是聚合物侧链与药物分子都有裸露的氨基，相似程度高，因此静电力相互作用的变化并不大。PCCL20与PmACCL20相比（接枝m-ALA 前后），接枝药物分子后，聚合物与5-ALA 的相互作用参数从−0.54 上升至1.35，其中范德华力相互作用参数从−0.49上升至0.14，静电力相互作用参数从−0.05 升至1.21，说明在PCCL20上接枝m-ALA 并不能提高5-ALA的包载效率，虽然m-ALA 接枝后与5-ALA 的结构也相似，但是接枝后原本裸露的羧基转变为酰胺键，使得聚合物与药物形成氢键的能力大大减弱，且在氢键较少的情况下，接枝后聚合物的体积变大，同样不利于药物与聚合物之间的相互作用。相比PCL20与5-ALA 之间的相互作用参数，PAPIL20与5-ALA 的相互作用参数更小，而PmACCL20与5-ALA 的相互作用参数变化较小。

表 2 纯无定形体系的密度和内聚能密度Table 2 Density and cohesive energy density of pure amorphous system

表 3 分子之间的Flory-Huggins 相互作用参数Table 3 Flory-Huggins interactional parameters of molecules

对比PCL20、PPIL20、PAPIL20、PCCL20、PmACCL205 种聚合物与m-ALA 的相互作用参数变化发现：在PPIL20上接枝5-ALA 也能提高两者的相容性，相容性变化的主要原因是范德华力的变化，而静电力变化相对较小；在PCCL20上接枝m-ALA 后，聚合物与m-ALA 之间的相容性几乎不变，这是范德华力相互作用参数与静电力相互作用参数的变化相互抵消的结果。与PCL20相比，PPIL20与PCCL20两种前药聚合物接枝药物后，和m-ALA 的相容性均变好。

对比PCL20、PPIL20、PAPIL20、PCCL20、PmACCL205 种聚合物与H2O 的相互作用参数发现：与H2O 的相容性最好的聚合物是PCCL20(χij=1.44)，参考mPEG与H2O 的相互作用参数(0.36)和PCL20与H2O 的相互作用参数(8.85)，PCCL20是一种水溶性较好的聚合物，因此虽然它与m-ALA 和5-ALA 的相互作用参数都比较好，但是不适用于载药胶束的疏水端。在PPIL20上接枝5-ALA 和在PCCL20上接枝m-ALA 后，聚合物的亲水性都变差，主要是由于静电力相互作用参数变化的原因。

综上分析，认为在PCCL20上接枝m-ALA 不利于载体对m-ALA 和5-ALA 的包载，在PPIL20上接枝5-ALA 能有效提高载体对5-ALA 的包载能力，因此选取PAPIL20对聚合物前药进行分析。为了保证聚合物疏水端的疏水性，采用PAPIL20与PCL20的无规共聚物作为疏水端。

2.2 药物接枝率不同的聚合物前药载体与药物之间的相互作用参数

以P(CL-co-APIL)为主要对象，研究药物接枝率对前药载体载药效率的影响。前药载体设计为：不同物质的量之比的PCL 与PPIL 的无规共聚物，PPIL 上的氨基均与5-ALA 中的羧基相连形成带有药物接枝的重复单元PAPIL 以及药物接枝率不同的聚合物前药载体，药物接枝率为非0 和非100%时采用无规的方式建立模型，部分模型见图5。

聚合物前药载体的药物接枝率（Y）定义为Y=y/(x+y)×100%，其中x 为PCL 的重复单元数，y 为PAPIL 的重复单元数。药物m-ALA 接枝率不同时P（CL-co-APIL）体系的密度和内聚能密度见图6，可以看出，随着药物m-ALA 接枝率的增大，聚合物的密度总体呈现逐渐上升的趋势。内聚能密度静电力分量逐渐上升，一方面是由于接入药物分子引入更多极性基团，另一方面是因为随着药物接枝率的上升，密度增加引起内聚能密度增加。内聚能密度的范德华力分量的变化较小，说明此时聚合物的内聚能变小，与接枝后聚合物主链之间的接触能力变小有关。但总的内聚能密度随着药物接枝率的上升而上升，说明其决定因素主要是静电力。

不同接枝率的聚合物前药与药物5-ALA之间的相互作用参数见图7。观察图7 可以看出无论是χij-vdw还是χij-e均随着5-ALA 接枝率的升高呈现下降的趋势。χij-vdw下降的主要原因可能是接枝药物分子后，聚合物侧链与药物的结构相似，使得其相容性变好。接枝药物后极性基团的数量大大增加，因此，χij-e随着药物接枝率的上升而下降，并且下降趋势几乎是线性的。两者共同作用使得χij也随着药物接枝率的上升而下降。

2.3 DPD 模拟结果

图 5 分子结构式转化为分子模型的示意图Fig. 5 Schematic diagrams of the transformation of molecular structure into molecular models

图 6 药物m-ALA 接枝率不同时P(CL-co-APIL)的密度(a)及内聚能密度(b)Fig. 6 Density (a) and cohesion energy density (b) of P(CL-co-APIL) with different grafting ratios of m-ALA

为了更加直观地观察聚合物前药胶束的形态和载药状态，采用DPD 模型进行介观模拟，排斥参数的计算与前文的Flory-Huggins 相互作用参数有关，珠子之间的Flory-Huggins 相互作用参数及排斥参数分别见表4、表5。

不同模拟步数下聚合物的DPD 模拟的自组装过程分别如图8～10 所示，模拟过程图显示不同接枝率的聚合物前药均能形成球形载药胶束，胶束形态差别并不大，说明在该体系范围内，亲水性的增大不影响胶束的形成。由0～20 000 步的变化过程可以看出，胶束首先是形成一些较小的聚集体，聚集体相互接触融合形成一个大的胶束，15 000 步以后，体系基本平衡，胶束形态稳定。对比不同接枝率的聚合物前药载体可以看出，在药物接枝率为0 的mPEG-b-PCL20/5-ALA 组(图8)，平衡后药物的分布仍然比较均匀，胶束外有大量的药物分子(红色)，平衡前后药物分布状态差别不大。药物接枝率为40%的mPEGb-P(CL12-co-APIL8)/5-ALA 组(图9)，在0 步时药物分子的分布比较分散，与接枝率为0 时相似；平衡后药物分子在胶束外的分布明显减少；当药物接枝率为100%时(图10)，平衡后药物分子在胶束外几乎没有分布，且在5 000 步时，药物分子几乎已经全部包载进胶束。以上结果说明药物接枝率的提高能有效改善载体的载药效率。

5-ALA 接枝率不同时聚合物前药载体DPD 模拟自组装形成的胶束形态及剖面形态见图11。由总体形态图和截面图可以看出，在5-ALA 接枝率为0 时，形成了球形载药胶束，5-ALA 在胶束外占较大体积，在胶束内部分布较少，且主要分布在胶束亲水段mPEG 上，这是因为mPEG 与5-ALA 有较好的相容性(χij=0.302 1)；在5-ALA 接枝率为40%时，依然形成球形胶束，在胶束外5-ALA 尽管有分布，但相对接枝率为0 时分布仍较少，在胶束内核5-ALA 则分布均匀。当5-ALA 接枝率达到100%时，形成的胶束与其他组差别不大，但5-ALA 在胶束外几乎没有分布，同样均匀地分布在胶束内核，且相比其他组分布更加集中。因此随着5-ALA 接枝率的上升，胶束的载药效率升高。

图 7 药物接枝率不同时P(CL-co-APIL)聚合物前药与5-ALA 的Flory-Huggins 相互作用参数Fig. 7 Flory-Huddins interaction parameters of P(CL-co-APIL) and 5-ALA calculated by MD with different drug grafting ratios

表 4 各珠子之间的Flory-Huggins 相互作用参数Table 4 Flory-Huggins interaction parameters between beads

表 5 各珠子之间的排斥参数Table 5 Repulsive parameters between beads

图 8 不同模拟步数下mPEG-b-PCL20/5-ALA 的DPD 模拟自组装过程(红色为5-ALA，绿色为PCL，橙色为mPEG；Y（5-ALA）=0）Fig. 8 DPD simulation self-assembly of mPEG-b-PCL20/5-ALA at different simulation steps(Red: 5-ALA, Green: PCL, Orange: mPEG;Y(5-ALA)=0）

图 9 不同模拟步数下mPEG-b-P(CL12-co-APIL8)/5-ALA 的DPD 模拟自组装过程(红色为5-ALA，绿色为PCL，蓝色为PAPIL，橙色为mPEG； Y(5-ALA）=40%）Fig. 9 DPD simulation self-assembly of mPEG-b-P(CL12-co-APIL8)/5-ALA at different simulation steps(Red: 5-ALA, Green: PCL, Blue:PAPIL, Orange: mPEG; Y(5-ALA)=40%）

5-ALA 的密度分布(图12)和归一化密度分布曲线(图13)更能具体体现接枝率不同时5-ALA 的集中分布情况。图13 中的x 表示胶束轴向的坐标。从图12 中可以看出，在5-ALA 的接枝率为0 和20%时，5-ALA 在整个格子中的分布相对较均匀；随着5-ALA 接枝率的增加，5-ALA 在胶束内的密度越来越大，而在胶束外的密度变得越来越小。这个结果说明了5-ALA 接枝率的增加可以增大聚合物前药载体的载药效率。

图 10 不同模拟步数下mPEG-b-PAPIL20/5-ALA 的DPD 模拟自组装过程(红色为5-ALA，蓝色为PAPIL，橙色为mPEG；Y （5-ALA）=100%）Fig. 10 DPD simulation self-assembly of mPEG-b-PAPIL20/5-ALA at different simulation steps(Red: 5-ALA, Green: PCL, Orange:mPEG;Y(5-ALA）=100%）

图 11 5-ALA 接枝率不同时聚合物前药的DPD 模拟的形态图及剖面图(红色为5-ALA，绿色为PCL，橙色为mPEG)Fig. 11 DPD simulation morphology and section view of drug-loaded micelle with different 5-ALA grafting ratios (Red: 5-ALA, Green: PCL,Orange: mPEG)

图 12 5-ALA 接枝率不同时聚合物前药中5-ALA 的密度分布Fig. 12 Density distribution of 5-ALA in polymeric prodrug with different 5-ALA grafting ratios

图 13 5-ALA 接枝率不同时的归一化密度分布曲线Fig. 13 Normalized density distribution of 5-ALA with different grafting ratios

3 结 论

采用分子动力学模拟计算了不同聚合物前药与药物分子5-ALA 之间的Flory-Huggins 相互作用参数，将Flory-Huggins 相互作用参数分解为范德华力分量以及静电力分量，以此具体分析接枝药物对聚合物前药的影响。选用mPEG-b-P(CL-co-APIL)作为聚合物前药研究药物接枝率对药物载体载药效率的影响，并采用DPD 模拟研究聚合物前药载体的自组装过程。

（1）分子动力学模拟结果显示：在5-ALA 与P(CL-co-APIL) 的体系中，随着药物接枝率的增加，载体与药物之间的相容性变好。药物接枝一方面引入极性基团，对载体与药物之间的静电力产生影响；另一方面，药物接枝使得药物与聚合物前药之间的结构相似程度变高，对载体与药物之间的范德华力产生影响。

（2）DPD 模拟结果说明了药物接枝率不同的聚合物前药都能形成载药胶束，且随着药物接枝率的增加，药物分子的分布越来越集中于胶束的内核；亲水性的增加并不影响胶束的形成，因此药物接枝率的增加总体有利于载药量的提高。