跨尺度方法研究聚合物載藥膠束的定量構(gòu)性關(guān)系

吳文勝，張燦陽(yáng)，李秀喜，徐驍，錢(qián)宇，章莉娟

（1 華南理工大學(xué)化學(xué)與化工學(xué)院，廣東 廣州 510641；2 肇慶學(xué)院化學(xué)化工學(xué)院，廣東 肇慶 526061）

引 言

“單元操作”和“三傳一反”被公認(rèn)為化學(xué)工程發(fā)展進(jìn)程中的第一、二階段的重要標(biāo)志。劉錚等[1-2]的文獻(xiàn)中曾提及普林斯頓大學(xué)的James Wei 教授認(rèn)為“產(chǎn)品工程”是化學(xué)工程學(xué)科發(fā)展第三階段的重要里程碑。以產(chǎn)品需求為導(dǎo)向，開(kāi)發(fā)滿(mǎn)足所需性能的化學(xué)產(chǎn)品已成為化學(xué)工程研究發(fā)展的一個(gè)重要趨勢(shì)。

化學(xué)產(chǎn)品大致可分為分子產(chǎn)品、配方產(chǎn)品和結(jié)構(gòu)化產(chǎn)品[3-5]，其設(shè)計(jì)與開(kāi)發(fā)的核心部分是結(jié)構(gòu)和性能的關(guān)系[6-9]。目前，對(duì)于復(fù)雜結(jié)構(gòu)化產(chǎn)品的研究和開(kāi)發(fā)，一般是通過(guò)實(shí)驗(yàn)合成許多結(jié)構(gòu)類(lèi)似的分子，然后進(jìn)行性能評(píng)價(jià)，直到最終獲得合乎性能要求的產(chǎn)品。這種經(jīng)驗(yàn)性的研究方法存在研發(fā)周期長(zhǎng)、風(fēng)險(xiǎn)大和費(fèi)用高的不足。而通過(guò)計(jì)算機(jī)模擬研究結(jié)構(gòu)-性能關(guān)系也還主要是定性層面的，比如分子模擬和介觀模擬方法定性描述聚合物或膠體的化學(xué)形態(tài)、微觀形貌、介觀相分離等[10-11]。如果能定量關(guān)聯(lián)結(jié)構(gòu)-性能關(guān)系，則可從結(jié)構(gòu)預(yù)測(cè)產(chǎn)品最終性能，也可通過(guò)優(yōu)化結(jié)構(gòu)得到性能最優(yōu)的產(chǎn)品，從而大大減少產(chǎn)品開(kāi)發(fā)周期和資源消耗。這種把微觀參數(shù)與宏觀性能直接相關(guān)聯(lián)的方法稱(chēng)為跨尺度方法[9,12-13]。

本研究以結(jié)構(gòu)化產(chǎn)品-聚合物載藥納米膠束體系為研究對(duì)象，基于宏觀到微觀的跨尺度方法，并結(jié)合實(shí)驗(yàn)數(shù)據(jù)，研究其定量構(gòu)性關(guān)系（QSPR）模型。聚合物納米膠束是一類(lèi)極具發(fā)展前景的新型藥物輸送系統(tǒng)[14]。但由于聚合物分子較大，可以形成不同的立體構(gòu)象，因而不能像有機(jī)小分子那樣預(yù)測(cè)其三維結(jié)構(gòu)，也難以精確地合成、表征和描述其結(jié)構(gòu)特征。因此，至今未見(jiàn)有關(guān)聚合物膠束QSPR 研究的文獻(xiàn)報(bào)道。為了更好地反映聚合物結(jié)構(gòu)特征對(duì)聚合物膠束性質(zhì)的影響，本研究中將采用自相關(guān)法來(lái)描述聚合物的結(jié)構(gòu)信息[15]，提出表征聚合物結(jié)構(gòu)的新描述符-嵌段單元自相關(guān)（block unit autocorrelation, BUA）描述符的計(jì)算方法。基于這些新描述符，聚合物可描述成數(shù)字的形式，然后通過(guò)遺傳函數(shù)算法-多元線性回歸方法（GFA-MLR）進(jìn)行特征選擇和建模。遺傳函數(shù)算法用于變量選擇，可以同時(shí)得到大量的模型群，而不只是產(chǎn)生一個(gè)單一的模型，從而防止丟失大量有用的結(jié)構(gòu)信息[16]。采用Friedman 的LOF 函數(shù)[17-18]作為評(píng)價(jià)函數(shù)，并通過(guò)內(nèi)部驗(yàn)證、外部驗(yàn)證等方法建立最佳QSPR 模型群。本研究將有助于深入認(rèn)識(shí)影響聚合物膠束載藥量的結(jié)構(gòu)因素和作用機(jī)制，為聚合物及其載藥膠束體系的設(shè)計(jì)和性能預(yù)測(cè)提供一定的理論依據(jù)，新的自相關(guān)描述符計(jì)算方法和跨尺度建模方法也有望指導(dǎo)復(fù)雜結(jié)構(gòu)化產(chǎn)品的設(shè)計(jì)和開(kāi)發(fā)。

1 材料與方法

1.1 數(shù)據(jù)集

選用15 種星型聚合物及其自組裝膠束為研究體系。15 種星型聚合物包括 6 種四雜臂(PCL)2(PDEA-b-PPEGMA)2， 4 種 六 雜 臂(PCL)3(PDEA-b-PPEGMA)3， 3 種 四 均 臂(PCL-b-PDEA-b-PPEGMA)4和 2 種 六 均 臂(PCL-b-PDEA-b-PPEGMA)6星型聚合物，其中PCL（聚已內(nèi)酯）、PDEA（聚N,N-二乙基胺基甲基丙烯酸乙酯）和PPEGMA（甲基丙烯酸聚乙二醇單甲醚酯）分別為疏水性嵌段、pH 響應(yīng)嵌段和親水性嵌段。聚合物結(jié)構(gòu)見(jiàn)圖1，圖中4ASP1-M/H 是4 種雜臂/均臂星型聚合物的縮寫(xiě)。

聚合物的設(shè)計(jì)、合成和膠束載藥實(shí)驗(yàn)均由本研究組完成。載藥膠束在兩種配比（阿霉素/聚合物=10 mg/40 mg 和20 mg/40 mg）、pH=7.4 和25℃下制備，得到30個(gè)載藥量（LC）數(shù)據(jù)，分別見(jiàn)表1和表2[19-20]。這里對(duì)不同配比下的載藥膠束分別建模，一是為了確保制備載藥膠束的條件一致，二是希望能比較不同配比濃度下影響膠束載藥量的結(jié)構(gòu)因素是否相同。依據(jù)訓(xùn)練集和測(cè)試集的數(shù)據(jù)點(diǎn)必須相互靠近，以及訓(xùn)練集化合物盡可能多樣性的準(zhǔn)則[21]，將配比10 mg/40 mg 和配比20 mg/40 mg 的LC 數(shù)據(jù)集分別分割成訓(xùn)練集（9 個(gè)數(shù)據(jù)）和測(cè)試集（6 個(gè)數(shù)據(jù)，不能少于5 個(gè)數(shù)據(jù)[22]）。

1.2 嵌段單元自相關(guān)（BUA）描述符計(jì)算方法

20 世紀(jì)80年代初，Moreau 等[15]首次將數(shù)學(xué)上的自相關(guān)函數(shù)引申提出了自相關(guān)拓?fù)渲笖?shù)（autocorrelation topological index）?；趦捎H性嵌段聚合物分子可看作是由嵌段單元排列組合成的鏈條，提出將自相關(guān)拓?fù)渲笖?shù)法應(yīng)用于描述兩親性嵌段聚合物分子的結(jié)構(gòu)信息，稱(chēng)為嵌段單元自相關(guān)（BUA）描述符，通過(guò)對(duì)分子中幾種嵌段單元的組成、電子結(jié)構(gòu)特征、物理化學(xué)特征等廣度性質(zhì)的特征量加權(quán)獲得?；谶@些描述符，兩親性嵌段聚合物分子可描述成數(shù)字的形式。通過(guò)自相關(guān)拓?fù)渲笖?shù)的分析能夠找出某一個(gè)位置上的嵌段單元的特征值是否和與其鄰近的嵌段單元的特征值相互獨(dú)立或者是否有關(guān)，如果存在相互關(guān)系，那么相關(guān)的BUA描述符可以包含空間自相關(guān)的信息，這種信息反映二級(jí)結(jié)構(gòu)并區(qū)分嵌段單元順序上的差異。將常數(shù)項(xiàng)、接近常數(shù)的項(xiàng)和具有高度相關(guān)的描述符預(yù)刪除，剩余描述符用于后續(xù)的變量選擇。

圖1 星型聚合物的分子結(jié)構(gòu)Fig.1 Molecular structures of star polymers

以四雜臂星型聚合物(PCL)2(PDEA-b-PPEGMA)2為例，介紹BUA 描述符的計(jì)算過(guò)程。首先，將3 個(gè)嵌段（PCL，PDEA 和PPEGMA）的單元在B3LYP/6-311+G(d,p)水平通過(guò)Gaussian09 程序包進(jìn)行幾何結(jié)構(gòu)全優(yōu)化[23]，然后，通過(guò)Hyperchem、Gamess 和Materials Studio 程序包計(jì)算3 個(gè)嵌段單元的組成、電子結(jié)構(gòu)和物理化學(xué)特征值后計(jì)算出嵌段的特征值，最后通過(guò)式（1）計(jì)算BUA 描述符。

表2 配比20 mg/40 mg 下各模型的ln(LC)實(shí)驗(yàn)值和預(yù)測(cè)值Table 2 Predicted and experimental values of ln(LC) for each model at ratio 20 mg/40 mg (drug/polymer)

式中，BUA(l)fk代表嵌段間隔為l、聚合物分子特征為fk時(shí)的BUA 描述符；fki和fkj分別為嵌段i和嵌段j 的特征k 的值；N 為求和計(jì)算中的單元個(gè)數(shù)。圖2 為聚合物(PCL)2(PDEA-b-PPEGMA)2的拓?fù)浣Y(jié)構(gòu)，其最長(zhǎng)鏈有4 個(gè)嵌段長(zhǎng)度，即最大的嵌段間隔為3，這在15 種聚合物分子中是最長(zhǎng)鏈上的最小的嵌段間隔數(shù)，因而，聚合物每個(gè)特征量的自相關(guān)描述符有3 種，分別為一級(jí)（lag=1）、二級(jí)（lag=2）和三級(jí)指數(shù)描述符（lag=3）。此外，還有針對(duì)每個(gè)嵌段單元特征進(jìn)行簡(jiǎn)單求和、平均所得的兩個(gè)描述符。每個(gè)聚合物有56 種特征，每個(gè)特征有5 個(gè)描述符，因此每個(gè)聚合物分子共有280 個(gè)描述符。

1.3 描述符選擇和建模方法

圖2 四雜臂星型聚合物(PCL)2(PDEA-b-PPEGMA)2 拓?fù)浣Y(jié)構(gòu)Fig.2 Topological structure of four-miktoarm star polymers (PCL)2(PDEA-b-PPEGMA)2

采用遺傳函數(shù)算法（GFA）選擇描述符[24-25]，同時(shí)建立全局的多元線性回歸（MLR）模型，全部 計(jì)算在Materials Studio 程序包中完成。在各種描述符選擇方法中，GFA 因其具有自動(dòng)選擇優(yōu)化描述符參數(shù)，能同時(shí)建立多個(gè)模型，可以避免局部?jī)?yōu)化而達(dá)到全局優(yōu)化等優(yōu)點(diǎn)，目前已成為了一種應(yīng)用最為廣泛、高效、可靠的描述符選擇方法[24]。MLR 具有形式簡(jiǎn)單、表達(dá)清晰、易于解釋等特點(diǎn)，并且能夠在描述符和實(shí)驗(yàn)值之間建立明確的方程，是QSPR 研究中最為常用的建模方法之一[26]。

GFA-MLR 基本過(guò)程：首先，通過(guò)GFA 按照所需選擇的描述符個(gè)數(shù)產(chǎn)生一系列的初始種群，包含多個(gè)初始隨機(jī)解，即描述符的初始解；然后分別建立初始解的MLR 模型，并根據(jù)不同的評(píng)價(jià)函數(shù)對(duì)模型進(jìn)行打分，根據(jù)打分情況進(jìn)行選擇、交叉和變異等操作產(chǎn)生子代種群，得分高的解將有更多機(jī)會(huì)進(jìn)入子代種群，得分低的則被淘汰。GFA 將最好的幾個(gè)個(gè)體保存下來(lái)作為精華種群。經(jīng)過(guò)交叉和變異操作后，逐一比較新種群的個(gè)體和精華種群中的個(gè)體，如果新種群中存在更好的個(gè)體，就將其復(fù)制到精華種群中。這個(gè)過(guò)程循環(huán)往復(fù)直到優(yōu)化收斂[27]。

1.4 模型的性能參數(shù)與驗(yàn)證方法

建模時(shí)遵循經(jīng)濟(jì)合作與發(fā)展組織（OECD）的QSAR 建模五法則[28]：① 確定的終點(diǎn)；②明確的算法；③定義應(yīng)用域；④ 適當(dāng)驗(yàn)證模型擬合度、穩(wěn)健性和預(yù)測(cè)能力；⑤ 如果可能，進(jìn)行機(jī)理解釋。

模型的內(nèi)部預(yù)測(cè)能力通過(guò)留一法交叉驗(yàn)證（leave-one-out cross-validation）和均方根誤差（root mean squared error, RMSE）來(lái)評(píng)估。留一法交叉驗(yàn)證技術(shù)是在每一次循環(huán)時(shí)，隨意從訓(xùn)練集中刪去一個(gè)樣本，然后用剩余的樣本建立模型，然后用所建模型去預(yù)測(cè)刪除的樣本的性質(zhì)。重復(fù)這個(gè)過(guò)程直到所有的樣本被刪除一次，從而計(jì)算得到留一法交叉驗(yàn)證相關(guān)系數(shù)（）[29]。根據(jù)式（2）計(jì)算RMSE，用來(lái)衡量實(shí)驗(yàn)值與預(yù)測(cè)值之間的偏差。

式中，yi表示第i 個(gè)樣本的實(shí)驗(yàn)值，i表示模型的預(yù)測(cè)值；n 表示訓(xùn)練集中的樣本數(shù)。

為了避免偶然相關(guān)，進(jìn)一步采用Y-隨機(jī)性驗(yàn)證（Y-randomization）技術(shù)進(jìn)行模型的穩(wěn)健性評(píng)估[30]。該方法是將響應(yīng)矢量Y 隨機(jī)打亂順序，而保持描述符的順序不變，然后建立新的模型，重復(fù)500 次分別得到擬合相關(guān)系數(shù)平方（R2）和的平均值和。R2是訓(xùn)練集樣本實(shí)驗(yàn)值和預(yù)測(cè)值之間的相關(guān)系數(shù)平方，可作為評(píng)估模型擬合能力的指標(biāo)。如果（）遠(yuǎn)小于原來(lái)模型的性能參數(shù)R2（），則認(rèn)為原樣本數(shù)據(jù)中存在真正的QSPR 關(guān)系，所建模型穩(wěn)健性較好，原模型不存在“偶然相關(guān)”現(xiàn)象。

對(duì)所建模型進(jìn)行內(nèi)部驗(yàn)證后還必須進(jìn)行外部驗(yàn)證，外部驗(yàn)證能為模型預(yù)測(cè)能力提供一個(gè)更加嚴(yán)厲的評(píng)估。使用不參與建模的外部測(cè)試集（表1 和表2）進(jìn)行外部驗(yàn)證。使用驗(yàn)證參數(shù)Q2ext評(píng)估模型的外部預(yù)測(cè)能力，計(jì)算式如式（3）所示

1.5 定義應(yīng)用域

使用 leverage 方法定義模型的應(yīng)用域[31]。Leverage 值（h）定義如下

式中，xi表示待預(yù)測(cè)樣本i 的特征矩陣（即描述符矩陣），X 表示訓(xùn)練集樣本的特征矩陣，n 表示待預(yù)測(cè)樣本數(shù)。標(biāo)準(zhǔn)leverage 值定義為為模型中描述符的個(gè)數(shù)加l，m 為訓(xùn)練集樣本數(shù)），當(dāng)樣本的h 值高于h*時(shí)，則認(rèn)為該樣本在模型的應(yīng)用域之外（即X 例外點(diǎn)），那么所建模型對(duì)該樣本的預(yù)測(cè)不一定可靠，是模型外推的結(jié)果。當(dāng)樣本的標(biāo)準(zhǔn)殘差（σ）大于±3σ時(shí)，則認(rèn)為該樣本是Y 例外點(diǎn)。

以leverage 值為橫坐標(biāo)，標(biāo)準(zhǔn)殘差為縱坐標(biāo)構(gòu)成Williams 圖。圖中，兩條水平線表示3 倍標(biāo)準(zhǔn)殘差，垂直線表示h*值，只有落在y 軸與這3 條線圍成的區(qū)域內(nèi)的樣本才認(rèn)為符合模型的應(yīng)用域，其預(yù)測(cè)值才是可靠的。

2 結(jié)果與討論

2.1 數(shù)據(jù)集分析

通常，對(duì)于任何QSPR 模型來(lái)說(shuō)，數(shù)據(jù)集的質(zhì)量對(duì)所建模型的性能有很大的影響。為了建立一個(gè)通用的QSPR 模型，對(duì)數(shù)據(jù)集進(jìn)行預(yù)處理是必要的。從兩個(gè)方面對(duì)數(shù)據(jù)集進(jìn)行預(yù)處理：數(shù)據(jù)集（響應(yīng)值）的正態(tài)分布情況分析；數(shù)據(jù)集的離群值（outlier）分析。一般來(lái)說(shuō)，要求建立QSPR 模型的響應(yīng)值需符合正態(tài)分布規(guī)律。將訓(xùn)練集的LC 值轉(zhuǎn)化成倒數(shù)、平方根、平方、立方和自然底對(duì)數(shù)（ln）等多種形式后進(jìn)行Kolmogorov-Smirnov 檢驗(yàn)，發(fā)現(xiàn)LC 值轉(zhuǎn)化成ln(LC)時(shí)更符合正態(tài)分布，因此，建模以ln(LC)為響應(yīng)值。

表3 模型的性能參數(shù)Table 3 Performance parameters of models

對(duì)于整個(gè)數(shù)據(jù)集，應(yīng)用GFA-MLR 方法建立了大量的模型，通過(guò)分析模型的應(yīng)用域，沒(méi)有發(fā)現(xiàn)超出3 倍標(biāo)準(zhǔn)偏差的疑似響應(yīng)值離群（Y 離群）的樣本，因此不存在需要剔除的樣本。

2.2 QSPR 模型群的建立與性能分析

基于回歸分析的QSPR 是一種統(tǒng)計(jì)技術(shù)，通過(guò)增加模型中的變量可以提高模型的統(tǒng)計(jì)擬合，但也可能導(dǎo)致虛假、不穩(wěn)定或無(wú)效的模型。因此，QSPR模型的變量應(yīng)該有嚴(yán)格的要求，響應(yīng)值數(shù)與描述符數(shù)之比應(yīng)該越高越好，至少要達(dá)到5:1[33]。本研究中訓(xùn)練集是9 個(gè)樣本，因此下面研究的模型都是兩變量的模型。

分別對(duì)兩種配比下的ln(LC)訓(xùn)練集進(jìn)行建模。參數(shù)設(shè)置：種群大小為1000，最大迭代數(shù)為5000，精華種群大小為100，初始模型允許的最大變量數(shù)固定為2，變異概率為0.1，評(píng)價(jià)函數(shù)選擇Friedman 的LOF 函數(shù)。對(duì)于每種配比下產(chǎn)生的100 個(gè)模型分別計(jì)算出其各種性能參數(shù)，根據(jù)R2＞0.8，Q2ext＞0.6的條件，得到每種配比下由4 個(gè)最佳模型組成的模型群，共8 個(gè)模型（詳細(xì)性能驗(yàn)證參數(shù)和模型方程分別見(jiàn)表3 和表4）。

由表3 可知，8 個(gè)模型的回歸與擬合（R2＞0.854）和內(nèi)部驗(yàn)證（Qloo-cv＞0.651）均很理想。且RMSE＜0.089，說(shuō)明8 個(gè)模型有很好的內(nèi)部擬合度和預(yù)測(cè)能力。另外，R2Y-rand和Q2Y-rand都低于0.1，排除偶然相關(guān)的可能性，說(shuō)明模型具有很好的穩(wěn)健性能[30,34-35]。

由表3 中的外部驗(yàn)證參數(shù)Q2ext和表1、表2 中的模型預(yù)測(cè)結(jié)果可知，在測(cè)試集中的樣本都能很好地被預(yù)測(cè)（Q2ext＞0.629），表明8 個(gè)模型對(duì)測(cè)試集樣本有很好的預(yù)測(cè)能力[30]。另外，圖3 顯示，8 個(gè)模型的預(yù)測(cè)數(shù)據(jù)和實(shí)驗(yàn)數(shù)據(jù)之間也有較好的擬合曲線（這里僅顯示了兩個(gè)模型的圖，其他模型的擬合曲線類(lèi)似），模型的預(yù)測(cè)值ln(LC)均分布在實(shí)驗(yàn)值附近。

表4 不同配比下的最佳模型群Table 4 Best model groups at two different ratios

表5 列出了模型中各描述符和響應(yīng)值ln(LC)之間的相關(guān)系數(shù)矩陣。模型中變量間的相關(guān)系數(shù)都很小，說(shuō)明模型中描述符間不存在共線性和冗余信息。8 個(gè)模型中都有一個(gè)相同的描述符X1，另外分別為X2、X3、X4和X5，說(shuō)明不同配比下影響膠束載藥量的結(jié)構(gòu)因素是相同的。X2、X3、X4和X5之間的相關(guān)系數(shù)都很高，說(shuō)明它們之間具有相似的描述符空間，從圖3 中模型擬合曲線很相似也可以證明。X1出現(xiàn)在每個(gè)模型中，與ln(LC)間的相關(guān)系數(shù)分別為0.7763（配比10 mg/40 mg）和0.8364（配比20 mg/40 mg），比X2、X3、X4和X5前的系數(shù)大很多，說(shuō)明X1對(duì)ln(LC)的影響非常大。

表5 響應(yīng)值與模型群中自變量之間的相關(guān)系數(shù)矩陣Table 5 Correlation coefficient matrix between response value and variable in model groups

2.3 模型的應(yīng)用域

所建模型的應(yīng)用域能夠提供關(guān)于預(yù)測(cè)可靠性的信息，只有對(duì)樣本的預(yù)測(cè)值落在應(yīng)用域里才被認(rèn)為是可靠的。用Williams 圖分析模型的應(yīng)用域，見(jiàn)圖4（僅顯示兩個(gè)模型的Williams 圖，其他模型的圖類(lèi)似）。Williams 圖顯示，8 個(gè)優(yōu)化模型對(duì)外部測(cè)試集均具有較好的擬合能力和預(yù)測(cè)結(jié)果，其樣本都落在其應(yīng)用域內(nèi)，表明預(yù)測(cè)是準(zhǔn)確可靠的。

2.4 描述符分析和機(jī)理解釋

通過(guò)對(duì)優(yōu)化模型的內(nèi)部和外部性能評(píng)估，確定了8 個(gè)模型對(duì)訓(xùn)練集和測(cè)試集都有很好的擬合度、穩(wěn)定性、魯棒性和較強(qiáng)的預(yù)測(cè)能力，為模型的進(jìn)一步分析和解釋提供了可靠的保證。下面將深入討論和分析模型中特征的定義、描述符的含義，以及模型的構(gòu)效關(guān)系和機(jī)理闡釋。

2.4.1 模型中特征的定義及描述符分析 模型中出現(xiàn)了5 種描述符，分別為X1（FI(-)-M-O(mean)）、X2（FI(+)-M-O(lag=3)）、X3（HBD(lag=3)）、X4（Adsorp-E-DOX(lag=3)）和X5（HE(lag=3)），分別對(duì)應(yīng)5 種不同的特征量：羰基氧親電福井指數(shù)（FI(-)-M-O）、羰基氧親核福井指數(shù)（FI(+)-M-O）、氫鍵贈(zèng)體數(shù)（ HBD ）、 阿霉素吸附能（Adsorp-E-DOX）和水合能（HE）?！癿ean”是指各嵌段特征的平均值；“l(fā)ag=3”是指嵌段間隔為3。

圖3 模型的擬合曲線Fig.3 Fitting curves of models

圖4 模型的Williams圖Fig.4 Williams plots of models

FI(-)表示當(dāng)分子中某原子得到一個(gè)電子時(shí)，該原子上電荷的變化，定義為FI(-)=q-(r) -q(r)。q-(r)、q(r)分別表示帶一個(gè)負(fù)電荷分子、中性分子（不帶電荷）中的r 原子上Mulliken 自然電荷。FI(-)絕對(duì)值越大，說(shuō)明r 原子在作用過(guò)程中越易得到電子，具有強(qiáng)的親電作用。對(duì)于聚合物嵌段單元來(lái)說(shuō)，如果嵌段單元上某原子的FI(-)絕對(duì)值越大，說(shuō)明該原子具有越強(qiáng)的親電作用。對(duì)應(yīng)描述符X1是指聚合物中各嵌段羰基氧的FI(-)平均值，反映的是聚合物的整體親電能力。FI(+)表示當(dāng)分子中某原子失去一個(gè)電子時(shí)，該原子上電荷的變化，定義為FI(+)= q(r) -q+(r)。q+(r)表示帶一個(gè)正電荷的分子中r 原子的Mulliken自然電荷。FI(+)絕對(duì)值越大，說(shuō)明r 原子越容易失去電子，具有強(qiáng)的親核作用[36]。對(duì)于聚合物嵌段單元來(lái)說(shuō)，如果嵌段單元上某原子的FI(+)絕對(duì)值越大，說(shuō)明該原子親核作用越強(qiáng)。對(duì)應(yīng)描述符X2是指聚合物中各嵌段羰基氧的FI(+)的三級(jí)BUA 描述符，說(shuō)明聚合物各嵌段親核能力在三級(jí)自相關(guān)拓?fù)浣Y(jié)構(gòu)上是有相互關(guān)系的，包含了空間自相關(guān)信息，反映其二級(jí)結(jié)構(gòu)上嵌段順序上的差異。

HBD 是表示嵌段單元中氫鍵的贈(zèng)體數(shù)量。氫鍵贈(zèng)體數(shù)越多，結(jié)合受體數(shù)形成氫鍵的數(shù)量就越多，對(duì)應(yīng)描述符X3是指各嵌段氫鍵供體數(shù)的三級(jí)BUA描述，說(shuō)明聚合物各嵌段的氫鍵的贈(zèng)體數(shù)量在三級(jí)自相關(guān)拓?fù)浣Y(jié)構(gòu)上是有相互關(guān)系的。Adsorp-E-DOX為嵌段單元吸附阿霉素的能量大小，其絕對(duì)值越大，則該分子對(duì)阿霉素的吸附作用力越大，對(duì)增加阿霉素的載藥量越有利，對(duì)應(yīng)描述符X4是指各嵌段阿霉素吸附能量的三級(jí)BUA 描述，說(shuō)明聚合物各嵌段的阿霉素吸附能量在三級(jí)自相關(guān)拓?fù)浣Y(jié)構(gòu)上是有相互關(guān)系的。HE 是物質(zhì)溶于大量水中成為無(wú)限稀釋溶液時(shí)所釋放的能量，通常該值越大，其水溶性越強(qiáng)，對(duì)應(yīng)描述符X5是指各嵌段水合能的三級(jí)BUA描述，說(shuō)明聚合物各嵌段的水合能在三級(jí)自相關(guān)拓?fù)浣Y(jié)構(gòu)上是有相互關(guān)系的。X3、X4和X5都包含了空間自相關(guān)信息，反映其二級(jí)結(jié)構(gòu)上嵌段順序上的差異。

由表4 可知，模型群的形式一樣，只是模型中描述符前的系數(shù)不同，說(shuō)明不同配比下聚合物結(jié)構(gòu)對(duì)載藥量的影響因素沒(méi)有改變。另外，每個(gè)模型都由與載藥量正相關(guān)的描述符X1和另外一個(gè)與載藥量負(fù)相關(guān)的描述符（分別為X2、X3、X4或X5）兩個(gè)自變量組成。X1前的系數(shù)比其他描述符前的系數(shù)大，說(shuō)明分子中羰基氧的親電性強(qiáng)弱對(duì)阿霉素的載 藥量影響最大。由于在模型中X1與載藥量是正相關(guān)，則X1增大，載藥量增加。從阿霉素分子的靜電勢(shì)圖（圖5）分析可知，氧和氮原子周?chē)患舜罅康呢?fù)電荷，而X1反映的是聚合物分子的親電性強(qiáng)弱，因此，X1越大，聚合物分子與阿霉素分子的靜電作用越強(qiáng)，越有利于聚合物膠束載藥量增加。

圖5 阿霉素分子的靜電勢(shì)圖Fig.5 Electrostatic potential diagram of DOX

模型中X2、X3、X4和X5描述符前的系數(shù)為負(fù)，表明其與阿霉素載藥量負(fù)相關(guān)，即增大描述符，載藥量減小。隨著X2增大，即羰基氧的親核性增強(qiáng)，不利于靜電作用吸附負(fù)電比較密集的阿霉素分子；阿霉素分子自身有很多氫鍵贈(zèng)體，它需要更多的氫鍵受體，因此，隨著聚合物氫鍵贈(zèng)體數(shù)X3增加，不利于與阿霉素的氫鍵形成；聚合物吸附阿霉素的能量X4減小，不利于藥物的吸附；隨著X5增大，即水化能增大，說(shuō)明聚合物水溶性增強(qiáng)，不利于疏水性阿霉素的包載。另外，這些特征量對(duì)載藥量的影響與聚合物的拓?fù)浣Y(jié)構(gòu)密切相關(guān)，將在后面進(jìn)行分析。

2.4.2 模型的構(gòu)性關(guān)系分析與機(jī)理闡釋 為了考察嵌段單元對(duì)載藥量的影響，表6 列出了聚合物嵌段單元的5 個(gè)特征值?！癈L”嵌段單元的“FI(-)-M-O”值為0.3871，遠(yuǎn)大于“DEAEMA”和“PEGMA”嵌段單元單體的0.0242 和0.0187，因此，“CL”嵌段單元結(jié)構(gòu)對(duì)X1影響最大，所以增加“CL”嵌段單元的數(shù)量有利于阿霉素載藥量的增加?！癋I(+)-M-O”、“HBD”和“Adsorp-E-Dox”3 個(gè)特征值比較接近，故3 種嵌段單元單體對(duì)X2、X3和X4的影響差別較小。對(duì)X5影響較大的是“DEAEMA”和“PEGMA”，兩者的“HE”特征值分別為8.01和6.91，遠(yuǎn)大于“CL”的1.49，這是由于“DEAEMA”中有氨基（氨基有利于增加其水合能），以及“PEGMA”較強(qiáng)親水性有關(guān)，因此增加這兩種嵌段單元數(shù)原則上是不利于載藥量增加。

表6 嵌段單元的特征值Table 6 Feature values of block units

聚合物拓?fù)浣Y(jié)構(gòu)對(duì)載藥量的影響將根據(jù)表7 中聚合物分子嵌段的聚合度、總聚合度和膠束載藥量數(shù)據(jù)進(jìn)行探討。對(duì)于四雜臂聚合物1-6 和六雜臂聚合物10-13，“PCL”嵌段越長(zhǎng)，“PPEGMA”嵌段越短，阿霉素的包載量越高，這主要是由于“PCL”嵌段的“FI(-)-M-O”值較大，為疏水性，有利于疏水性阿霉素的包載，而“PPEGMA”嵌段的“HE”值較大，為親水性，因而不利于載藥量增加。例如，6 個(gè)四雜臂聚合物中聚合物4 的“PCL”嵌段最長(zhǎng)（76 個(gè)），“PPEGMA”嵌段最短（18 個(gè)），其膠束載藥量最大（15.7%）。“DEAEMA”嵌段單元的增加也有利于載藥量的增加，例如，聚合物12 和13的“PCL”嵌段的聚合度相同，盡管聚合物12 的“PPEGMA”嵌段比聚合物13 的短，但由于聚合物13 的“DEAEMA”嵌段比12 的長(zhǎng)很多（前者63，后者36），結(jié)果是聚合物13 膠束的載藥量比聚合物12 的稍高。

對(duì)于四均臂聚合物7-9 和六均臂聚合物14-15，同樣是“PCL”嵌段越長(zhǎng)，“PPEGMA”嵌段越短，越有利于疏水性阿霉素的包載，但是“DEAEMA”嵌段單元的增加不利于載藥量的增加。例如，聚合物14 和15 的“PCL”嵌段和“PPEGMA”嵌段長(zhǎng)度相同，但“DEAEMA”嵌段長(zhǎng)度不同，分別為54 和90 個(gè)單元，而載藥量分別為9.7%和8.6%。主要可能是在雜臂聚合物中“DEAEMA”嵌段沒(méi)有與“PCL”嵌段相連成臂，而在均臂聚合物，“DEAEMA”嵌段是與“PCL”嵌段相連成臂，這種拓?fù)浣Y(jié)構(gòu)減弱了“DEAEMA”嵌段的疏水作用。由此可知，聚合物中不同嵌段在拓?fù)浣Y(jié)構(gòu)上位置的改變也會(huì)影響膠束的載藥量。

表7 聚合物分子中的嵌段單元數(shù)和總數(shù)及膠束載藥量Table 7 Number of block units, total and LC in polymers

3 結(jié) 論

以聚合物載藥膠束體系為案例，采用新的表征聚合物結(jié)構(gòu)特征的嵌段單元自相關(guān)（BUA）描述符的計(jì)算方法，建立了結(jié)構(gòu)-性能關(guān)系最佳模型群。模型經(jīng)內(nèi)外檢驗(yàn)參數(shù)評(píng)估、應(yīng)用域定義，確認(rèn)了模型的可靠性，確保了BUA 描述符對(duì)載藥機(jī)制的合理解釋。聚合物中羰基氧的親電作用對(duì)阿霉素載藥量的影響最大，親電作用越強(qiáng)，阿霉素的載藥量越大，反映在嵌段單元上是“CL”越多，阿霉素載藥量越大；對(duì)于星型聚合物，增加“PCL”嵌段長(zhǎng)度或減小“PPEGMA”嵌段長(zhǎng)度，都有利于聚合物膠束阿霉素載藥量的增加；增加星型雜臂聚合物的“PDEAEMA”嵌段長(zhǎng)度有利于載藥量增加，而對(duì)于均臂聚合物則不利于載藥量增加。通過(guò)對(duì)模型中BUA 描述符的分析和機(jī)理闡釋?zhuān)M(jìn)一步闡釋了聚合物嵌段單元結(jié)構(gòu)和拓?fù)浣Y(jié)構(gòu)對(duì)阿霉素載藥量的影響。研究顯示，BUA 描述符能很好地在模型中反映聚合物的嵌段結(jié)構(gòu)和拓?fù)浣Y(jié)構(gòu)對(duì)載藥量的影響，因此在設(shè)計(jì)目標(biāo)聚合物時(shí)，可以根據(jù)需求調(diào)控分子結(jié)構(gòu)，也可以通過(guò)模型預(yù)測(cè)載藥量。該研究為定量設(shè)計(jì)和開(kāi)發(fā)合適載藥量的聚合物及其膠束提供一種新的思路和理論方法，有望指導(dǎo)其他復(fù)雜結(jié)構(gòu)化產(chǎn)品的設(shè)計(jì)和開(kāi)發(fā)。

[1]Liu Zheng(劉錚), Jin Yong(金涌), Wei Fei(魏飛), Li Yourun(李有潤(rùn)), Luo Guangsheng(駱廣生), Yuan Naiju(袁乃駒).Recent development of chemical engineering science [J].Chemical Industry and Engineering Progress(化工進(jìn)展), 2002, 21(2): 87-91.

[2]Li Jinghai (李靜海).Perspectives on chemical engineering in the 21st century [J].Journal of Chemeical Industry and Engineering (China) (化工學(xué)報(bào)), 2008, 59(8): 1879-1883.

[3]Qian Yu(錢(qián)宇), Pan Jizheng(潘吉錚), Jiang Yanbin(江燕斌), ZhangLijuan(章莉娟), Ji Hongbing(紀(jì)紅兵).Chemical product engineering theories and technologies [J].Chemical Industry and Engineering Progress(化工進(jìn)展), 2003, 22(3): 217-223.

[4]Li Jinghai(李靜海), Hu Ying(胡英), Yuan Quan(袁權(quán)), He Mingyuan(何鳴元).Perspectives on the 21st Chemical Engineering (展望21世紀(jì)的化學(xué)工程)[M].Beijing: Chemical Industry Press, 2004.

[5]Cussler E L, Moggridge G D.Chemical Product Design[M].Cambridge: Cambridge University Press, 2011.

[6]Sun Hongwei(孫宏偉).Focus of chemical engineering-understanding space-time multi-scale structure [J].Chemical Industry and Engineering Progress(化工進(jìn)展), 2003, 22(3): 224-227.

[7]Yuan Weikang(袁渭康).Miniaturization of target scale: refreshing chemical engineering philosophy [J].Chemical Industry and Engineering Progress(化工進(jìn)展), 2004, 23(1): 9-11.

[8]Guo X D, Zhang L J, Qian Y.Systematic multiscale method for studying the structure-performance relationship of drug-delivery systems [J].Ind.Eng.Chem.Res., 2012, 51: 4719-4730.

[9]Yang Ning(楊寧), Li Jinghai(李靜海).Mesoscience in chemical engineering and virtual process engineering: analysis and perspective [J].CIESC Journal(化工學(xué)報(bào)), 2014,65(7): 2403-2409.

[10]Zheng L S, Yang Y Q, Guo X D, Sun Y, Qian Y, Zhang L J.Mesoscopic simulations on the aggregation behavior of pH-responsive polymeric micelles for drug delivery [J].J.Colloid Interf.Sci., 2011, 363(1): 114-121.

[11]Chen H, Ruckenstein E.Formation and degradation of multicomponent multicore micelles: insights from dissipative particle dynamics simulations [J].Langmuir, 2013, 29(18): 5428-5434.

[12]Li Bogeng(李伯耿), Luo Yingwu(羅英武).Production engineering : new developing space for chemical reaction engineering [J].Chemical Industry and Engineering Progress(化工進(jìn)展), 2005, 24(4): 337-340.

[13]Li Jinghai(李靜海), Hu Ying(胡英), Yuan Quan(袁權(quán)).Mesoscience: exploring old problems from a new angle [J].Scientia Sinica: Chimica(中國(guó)科學(xué):化學(xué)), 2014, 44(3): 277-281.

[14]Zhang C Y, Yang Y Q, Huang T X, Zhao B, Guo X D, Wang J F, Zhang L J.Self-assembled pH-responsive MPEG-b-(PLA-co-PAE) block copolymer micelles for anticancer drug delivery [J].Biomaterials, 2012, 33(26): 6273-6283.

[15]Moreau G, Broto P.The auto-correlation of a topological-structure—a new molecular descriptor [J].Nouveau Journal De Chimie-New Journal of Chemistry, 1980, 4(6): 359-360.

[16]Shahlaei M.Descriptor selection methods in quantitative structure- activity relationship studies: a review study [J].Chem.Rev., 2013, 113(10): 8093-8103.

[17]Rogers D, Hopfinger A J.Application of genetic function approximation to quantitative structure-activity relationships and quantitative structure-property relationships [J].J.Chem.Inf.Comp.Sci., 1994, 34(4): 854-866.

[18]Friedman J H.Multivariate adaptive regression splines (with discussions) [J].Ann.Stat., 1991, 19: 1-67.

[19]Yang Y Q, Zhao B, Li Z D, Lin W J, Zhang C Y, Guo X D, Wang J F Zhang L J.pH-Sensitive micelles self-assembled from multi-arm star triblock co-polymers poly(ε-caprolactone)-b-poly(2-(diethylamino) ethyl methacrylate)-b-poly(poly (ethylene glycol) methyl ether methacrylate) for controlled anticancer drug delivery [J].Acta Biomater, 2013, 9(8): 7679-7690.

[20]Lin W J, Nie S Y, Zhong Q, Yang Y Q, Cai C Z, Wang J F, Zhang L J.Amphiphilic miktoarm star copolymer (PCL) 3-(PDEAEMA-b- PPEGMA) 3 as pH-sensitive micelles in the delivery of anticancer drug [J].J.Mater.Chem.B, 2014, 2(25): 4008-4020.

[21]Gramatica P, Giani E, Papa E.Statistical external validation and consensus modeling: a QSPR case study for KOCprediction [J].J.Mol.Graph.Model., 2007, 25(6): 755-766.

[22]Golbraikh A, Tropsha A.Beware of q2! [J].Journal of Molecular Graphics and Modelling, 2002, 20(4): 269-276.

[23]Carbonniere P, Thicoipe S, Pouchan C.Theoretical strategy to build structural models of microhydrated inorganic systems for the knowledge of their vibrational properties: the case of the hydrated nitrate aerosols [J].J.Phys.Chem.A, 2013, 117(18): 3826-3834.

[24]Xi L L, Sun H J, Li J Z, Liu H X, Yao X J, Gramatica P.Prediction of infinite-dilution activity coefficients of organic solutes in ionic liquids using temperature-dependent quantitative structure-property relationship method [J].Chem.Eng.J., 2010, 163(3): 195-201.

[25]Cho S J, Hermsmeier M A.Genetic algorithm guided selection: variable selection and subset selection [J].J.Chem.Inf.Comp.Sci., 2002, 42(4): 927-936.

[26]Lei B L, Ma Y M, Li J Z, Liu H X, Yao X J, Gramatica P.Prediction of the adsorption capability onto activated carbon of a large data set of chemicals by local lazy regression method [J].Atmos.Environ., 2010, 44(25): 2954-296.

[27]Choi J H, Ahn H, Han I.Utility-based double auction mechanism using genetic algorithms [J].Expert Syst.Appl., 2008, 34(1): 150-158.

[28]Chakravarti S K, Saiakhov R D, Klopman G.Optimizing predictive performance of CASE Ultra expert system models using the applicability domains of individual toxicity alerts [J].J.Chem.Inf.Model., 2012, 52(10): 2609-2618.

[29]Schüürmann G, Ebert R U, Chen J, Wang B, Kühne R.External validation and prediction employing the predictive squared correlation coefficient test set activity mean vs training set activity mean [J].J.Chem.Inf.Model., 2008, 48(11): 2140-2145.

[30]Roy K, Kabir H.QSPR with extended topochemical atom (ETA) indices: modeling of critical micelle concentration of non-ionic surfactants [J].Chem.Eng.Sci., 2012, 73: 86-98.

[31]Tropsha A, Gramatica P, Gombar V K.The importance of being earnest: validation is the absolute essential for successful application and interpretation of QSPR models [J].QSAR Comb.Sci., 2003, 22(1): 69-77.

[32]Qin Litang(覃禮堂), Liu Shushen(劉樹(shù)深), Xiao Qianfen(肖乾芬), Wu Qingsheng(吳慶生).Internal and external validtions of QSAR model: review [J].Environmental Chemistry(環(huán)境化學(xué)), 2013, 7(32): 1205-1211.

[33]Yan D, Jiang X, Yu G, Zhao Z, Bian Y, Wang F.Quantitative structure-toxicity relationships of organophosphorous pesticides to fish (Cyprinuscarpio) [J].Chemosphere, 2006, 63(5): 744-750.

[34]Zhu H S, Shen Z M, Tang Q L, Ji W C, Jia L J.Degradation mechanism study of organic pollutants in ozonation process by QSAR analysis [J].Chem.Eng.J., 2014, 255: 431-436.

[35]Roy P P, Roy K.On some aspects of variable selection for partial least squares regression models [J].QSAR Comb.Sci., 2008, 27(3): 302-313.

[36]Parr R G, Yang W.Density-functional Theory of Atoms and Molecules[M].Oxford: Oxford University Press, 1989.