HA-CS复合物界面作用与力学性能模拟

袁秋华，万 磊，石 鑫，林亚宁，徐安平，张自强，陈泽汇，张培新

1)深圳大学化学与环境工程学院，广东深圳 518060；2)中国电子科技集团公司第九研究所，四川绵阳621000

羟基磷灰石(hydroxyapatite, HA；Ca10(PO4)6OH2)是目前应用最广泛的骨修复材料，它是自然骨的主要成分，约占人体骨质量的60%[1]．其他的骨修复材料，比如不锈钢和钛合金等，由于其惰性不能很好地与生物体本身产生相容性，作为硬组织植入材料，它们与骨之间只是一种机械嵌连的骨整合，而非化学骨性结合，致使植入后与骨组织之间结合较差，常引起植入失效．羟基磷灰石具有良好的生物相容性和生物活性[2-3]，它在人体内可分解出Ca元素和P元素修复骨骼或牙齿．但羟基磷灰石机械性能不佳，抗弯强度和断裂韧性指标均低于人体致密骨，限制了它们单独在人体负重部位的使用．然而，羟基磷灰石可以与多种材料进行复合来改善它的力学性能[4-5]，因此，近年来羟基磷灰石及其复合物的研究受到了广泛关注．

壳聚糖(chitosan, CS)是生物医学领域应用十分广泛的一种材料，它具有优异的生物相容性、生物活性和生物可降解性[6]．同时它也是一种优秀的骨修复材料．LEE等[7]采用壳聚糖/磷酸三钙作为三维支架体外培养成骨细胞，扫描电镜检查显示骨细胞与材料吸附较好，细胞生长良好，证明了壳聚糖的复合物可作为骨细胞培养的支架材料．PAL等[8]以羟基磷灰石复合壳聚糖材料作为骨移植的生物材料修复兔胫骨缺损，实验证明，羟基磷灰石中加入壳聚糖可增加羟基磷灰石颗粒的赋形型．

分子动力学(molecular dynamics，MD)模拟为研究无机材料与高聚物的相互作用提供了一个有效的方法．MD是一种计算机辅助实验技术，它利用计算机以原子水平的分子模型来模拟分子的结构与行为，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质，是当前分子模拟的重要研究手段．ZHANG等[9]利用分子动力学模拟研究了在没有耦合剂作用的条件下，羟基磷灰石与其他高聚物组成的复合材料界面的相互作用，比较分析了聚合物在3个不同平面的约束能，对人工骨制造具有一定的指导作用．SANDOVAL等[10]运用分子动力学模拟研究CS分别与聚乙烯醇(polyvinyl alcohol,PVA)、聚甲基丙烯酸羟乙酯(hydroxyethyl methacrylate，HEMA)共混物中基团间的相互作用，以确定这两种共混物中PVA和HEMA的羟基在与CS的羟基相互作用中占主导地位，还与CS中的氨基作用占主导地位．

本文采用分子动力学计算方法，研究HA和CS复合物的界面作用及力学性质，以期发现一种优良的HA-CS复合结构．

1 实 验

1.1 羟基磷灰石-壳聚糖复合膜的制备

采用化学沉淀法合成羟基磷灰石粉体. 以Ca(NO3)2·4H2O和Na2HPO4为原料，按n(Ca)∶n(P)=1.67分别配成0.5 mol/L的溶液，用恒流滴定泵以2 mL/min滴定速度将Na2HPO4溶液滴加到Ca(NO3)2·4H2O溶液中，用1 mol/L NaOH溶液调整pH值在10～11，在40 ℃下磁力搅拌，搅拌速度为300 r/min，待反应完全，恒温搅拌2 h，静置沉化24 h，将沉淀物真空抽滤、洗涤，至pH值呈中性，得到的滤饼置于90 ℃烘箱中干燥12 h.烘干的样品经研磨后置于700 ℃马弗炉中煅烧3 h，得到羟基磷灰石粉体.合成HA的反应式为

10Ca(NO3)2·4H2O+6Na2HPO4+2NaOH=

Ca10(PO4)6(OH)2+14NaNO3+

6HNO3+40H2O

(1)

采用溶液混合法制备HA-CS悬浮液.取1.0 g壳聚糖(脱乙酰度≥95%)溶于30 mL质量分数为2%的醋酸溶液中，50 ℃条件下磁力搅拌2 h.待CS完全溶解后，加入1.5 g经200目(约75 μm)筛子过筛后的HA粉末，超声分散0.5 h，恒温搅拌6 h，得到HA-CS悬浮液.

采用旋涂法制备HA-CS复合膜.将处理过的不锈钢板光滑面朝上吸附于SC-1B匀胶机旋转台上，用胶头滴管取配好的HA-CS悬浮液在钢板表面滴加1～2滴，以800和2 500 r/min的启动速度和恒定速度分别转动10和30 s.如此重复涂膜3～5次，然后取出钢板，分别在25和37 ℃条件下干燥12 h，最终得到成型的HA-CS复合膜．

1.2 扫描电子显微镜测试

采用日本日立S-3400N(Ⅱ)型扫描电子显微镜(scanning electron microscopy,SEM)观测复合膜表面形貌，结果如图1. 其中，图1(a)和(c)为25 ℃条件下干燥制得的HA-CS复合膜的不同放大倍率SEM照片；图1(b)和(d)是37 ℃条件下干燥制得HA-CS复合膜的SEM图. 图1(a)和(b)是放大500倍复合膜的表面形貌，可以看出CS与HA黏结较好， 结合致密， 复合膜表面无裂纹， 无机相均匀地分散在有机相中，且无明显界限.此外，图1(b)中复合膜表面较图1(a)更为平整、均匀，游离的HA颗粒相对较少.图1(c)和(d)是放大5 000倍的复合膜表面微结构图，可以清晰看到CS包覆HA颗粒形成的致密结构，HA颗粒团聚较少，并且复合膜表面形成许多微孔结构，有利于引导成骨细胞黏附和生长，这也是骨替代材料的一个必要特性.相比图1(c)，图1(d)中HA晶粒尺寸更小，在CS中分散更加均匀，并且复合膜表层的CS相对减少，更多融入内层，表明接近人体温下(37 ℃)的干燥条件更适宜理想的HA-CS复合膜成型.

图1 不同干燥条件下HA-CS复合膜的SEM图Fig.1 SEM images of HA-CS composite coatings under different drying conditions

2 模拟方法及模型

COMPASS(condensed-phase optimized molecular potential for atomistic simulation studies)力场具有较好的普适性，是第1个基于量子力学从头计算和经验数据参数化的力场，已经被证明能应用于研究有机和无机体系[11]．相关研究表明，COMPASS力场能够用于准确预测多种孤立分子和凝聚相的结构、物理和化学性质[12-13]．本研究通过原子模拟凝聚态优化分子势能力场(condensed-phase optimized molecular potentials for atomistic simulation studies, COMPASS)在Material Studio(MS)模拟平台(Accerly, San Diego, CA)上进行模拟．COMPASS力场中的相互作用如式(2)[14]．

Epot=E1+E2+…+E13

(2)

其中，

Epot为键合物的能量；E1为键伸缩能；E2为键弯曲能；E3为键扭转能；E4为键角面外弯曲能；E5～E11为相互耦合能；E12为库仑能；E13为范德华能相互作用；θ和θ0分别为键角和平衡键角；i为电荷受体，j为电荷给体；K、H、F、m和V为力场参数，b、h和φ分别为键长、弯曲角和扭转角[12];rij为相互作用电荷间的范德华半径．COMPASS力场已成功应用于各种材料研究，如LI等[14]对于Al2O3界面上的正构烷烃的吸附性能的研究；TUNG等[15]通过分子动力学模拟技术研究了溶剂种类对无规聚甲基丙烯酸甲酯链在溶液中的动力学的影响；PRATHAB等[16]研究的聚合物-金属氧化物的相互作用．同时，COMPASS已经成功地应用于纳米级系统的研究，如FERNANDO等[17]报道的氧化钛纳米粉体的相态变化，TOTH等[18]研究的聚合物纳米复合材料体系的相变．

由于HA空间群为p21/b空间群，其晶胞参数为a=b=0.943 2 nm，c=0.688 1 nm[19]．在HA晶体中，OH-沿着c轴方向排成一列，每个OH-同3个Ca原子配位从而构成位于同一平面的三角形．该Ca原子三角形沿着c轴平移1/4轴长度后旋转60°即可与另一平面的Ca原子三角形重合，沿着c轴俯视HA晶体，可观察到Ca原子构成的六边形．基于COMPASS力场对研究中HA晶胞中所有原子力场参数进行赋值，参考文献[20]研究成果对HA晶胞中其余部分原子电荷进行赋值(正负电荷参数)如下：Ca2+为+2.0，P5+为+2.6，H+为+0.6，PO43-中O为-1.4，OH-中O为-1.6(单位晶胞含有10个Ca2+、6个 PO3-和2个OH-)．为验证模拟的准确性，将基于COMPASS力场的MD模拟得到的HA晶胞参数及磷原子的坐标值与文献报道理论值对比．MD模拟后的晶胞的晶格长度a=b=0.948 6 nm，c=0.677 8 nm；a与b之间夹角为 90°，a与c之间夹角为 120°，这与文献[19]报道的数值非常接近．MD模拟后HA晶胞中P原子的分子内坐标为(0.382,0.377,0.248)，与文献[19]报道的数据(0.399,0.369,0.250)也非常接近．该结果表明，利用COMPASS力场对HA的模拟是有效的．

在MD模拟中，沿着(001)、(100)和(110)三个平面切割，并将最小化周期边值的无键截断距离设置为0.95 nm，从而建立HA晶体表面模型[21]．最小化周期边值的过程中依次采用收敛标准为4.184×104kJ/(mol·nm)的最速下降法、收敛标准为4.184×102kJ/(mol·nm)共轭梯度法以及收敛标准为4.184 kJ/(mol·nm)的BFGS算法对模型进行优化．优化后所得HA表面尺寸参数分别为，HA(001)：3.2 nm×2.8 nm×1.1 nm；HA(100)为2.8 nm×2.0 nm×1.0 nm；HA(110)为3.2 nm×2.0 nm×1.1 nm．HA表面模拟结果如图2．

图3 优化后的CS的模型结构Fig.3 Optimized CS simulated structure

壳聚糖链使用Materials Studio 模拟平台中的Amorphous 模块建立．壳聚糖密度设为1.4 g/cm3，与壳聚糖的实际密度相符．壳聚糖链由80个单体组成，建好模型后，采用分子动力学对其进行优化以得到稳定的壳聚糖链式结构．在计算中采用等温等容系统，时间步长为0.1 fs，在300 K的温度下进行50 ps的分子动力学计算．图3展示了优化后的壳聚糖链模型．

在建立并优化好壳聚糖链和羟基磷灰石晶体基础上，建立羟基磷灰石-壳聚糖的复合结构．首先将壳聚糖链置于羟基磷灰石的表面上，然后在复合层的表面上添加10 nm厚的真空层，以保证壳聚糖链只与羟基磷灰石的一个面反应．建好复合模型后对其进行了分子动力学模拟，系综设为等温等容系统，时间步长为0.01 fs，总步数时间长为30 ps．图4和图5分别展示了HA-CS复合物进行MD模拟前后的状态．

图4 HA各晶面与CS在MD模拟前结构Fig.4 Interactions between HA surfaces and CS before MD simulation

图5 HA各晶面与CS在MD模拟后结构Fig.5 Interactions between HA surfaces and CS after MD simulation

3 结果与讨论

HA-CS复合物之间的附着力可通过计算HA晶体表面和CS分子之间的结合能得到．由式(3)计算得到的CS聚合物链和HA表面的结合能Ebinding， 如表1．

Ebinding=(EHA+ECS-EHA-CS)/N

(3)

其中，EHA为HA表面的能量；ECS为CS聚合链的能量；EHA-CS为HA-CS复合物的能量；N为每个聚合物链上的原子数．Ebinding越大，表示复合物之间的结合能越大[22]．研究发现，对于HA-CS体系，CS与HA(110)表面的结合能最大，与HA(001)表面的结合能最小．这是因为HA(110)面表面的原子密度最大，也导致其对CS的结合能最大．因此，HA(110)面是HA与CS复合的最优表面．

对固体材料而言，弹性性能是材料基础的固态性质，且与硬度、塑性和晶界能等参数密切相关，是了解材料“脾气”的重要参数．

表1 HA各晶面与壳聚糖复合后结合能的大小

应力分为正应力和切应力，根据合金受力情况的不同可以把弹性模量分为杨氏模量(E)、 剪切模量(G)和体积模量(B)． 所有模量的定义都采用比值的方式．把正应力(沿纵向拉应力和压应力的总称)与正应变的比值称为杨氏模量，所以，杨氏模量是纵向的弹性模量，表征了合金对压应力和拉应力的抵抗能力，其值越小说明材料越易拉伸．

切应力与切应变的比值定义为剪切模量． 它反映了材料对切应力抵抗能力的强弱．其值越小说明抵抗切应力的能力越弱，材料越易发生切应变．体积应力与体积应变的比值即为体积模量． 其中体积应力是指材料整体受到的压强，体积应变是指材料变化的体积与原体积的比值．体积模量较小，说明材料的体积应变越容易．

根据晶体弹性力学的理论，应力和应变满足矩阵关系．直观上看在弹性常数的矩阵中有36个弹性常数，由于晶体的结构和对称性密切相关，所以很多物理特性都与对称性有关，弹性常数也是这样，在弹性常数矩阵中Cij=Cji， 所以最多有21个独立的弹性常数，由于不同晶系的对称性不同，独立的弹性常数的个数也不相同．其独立的弹性常数的个数，随着晶系对称性的提高而减少．对六角晶系的HA而言，有5个独立的弹性常数C11、C12、C13、C33和C44． 计算得到的HA(001)、HA(100)和HA(110)面结构与HA-CS复合膜的弹性性能值如表2．式(4)为立方晶系弹性矩阵，其中，σ1～σ3和τ4～τ6为应力分量；ε1～ε3和γ4～γ6为应变分量．

表2 HA各晶面和HA-CS复合结构的弹性性能

(4)

Voigt-Reuss-Hill(V-R-H)是常用的计算材料弹性模量方法．对立方晶体来说，式(5)至(7)为V-R-H方法下的剪切模量表达式．

(5)

(6)

(7)

体积模量是一种十分重要的弹性模量，可以表征材料的多种性能．而且杨氏模量E、 泊松比ν与体积模量B及剪切模量G之间也存在一定的关系，见式(8)和式(9)．

(8)

(9)

HA各晶面及HA各晶面与CS结合后在三维方向上的杨氏模量如图6．

由于羟基磷灰石的机械性能不佳，抗弯强度和断裂韧性均较差，不能满足人体骨的强度需求．从羟基磷灰石的弹性性能分析，这是由于HA的剪切模量和杨氏模量相对较低．

从表1可见，3个面的杨氏模量均呈现较大的各向异性，不利于羟基磷灰石的应用．从表2中列出的羟基磷灰石不同表面和壳聚糖链复合模型的弹性性能分析，可以看出HA(001)、HA(100)和HA(110)3个面的弹性性能相似，HA(001)的剪切模量G较大，但是在x方向上杨氏模量偏小．当羟基磷灰石与壳聚糖复合后， HA(001)面出现了负的弹性常数(C13和C44)，说明该复合结构的机械结构不稳定．对杨氏模量来说，在x方向上的杨氏模量增大，但是z方向上E减小，整体的各向异性没有改善．对于HA(100)面，弹性性能出现了明显的下降，各项弹性常数减小50%以上，杨氏模量也大幅度下降，虽然各向异性较低，但是杨氏模量太小，机械性能极差．对于HA(110)面，发现复合结构的剪切模量由29.2 GPa升至42.4 GPa，杨氏模量由92.0、62.7和57.0 GPa分别升至95.6、102.2和97.5 GPa，3个方向上的杨氏模量差量最大值由46.7%到7.0%，因此可克服因材料各向异性而导致的材料缺陷．通过以上分析，可以看出CS能够对HA(110)面的力学性能起到明显的加强作用，即在应用HA-CS复合材料时，需要选择HA(110)面与CS进行复合，以期达到最好的效果．

结 语

本研究采用分子动力学模拟的方法对HA-CS复合膜体系的性质进行研究．分析发现，羟基磷灰石不同表面与壳聚糖的结合能大小不同，其大小关系为HA(110)>HA(100)>HA(001)，结合能大小与表面原子密度相关：表面原子密度越大，其结合能也越大，即原子密度的大小可以反应结合能的大小．

HA(110)面与壳聚糖复合后力学性能最佳，复合结构的剪切模量由29.2 GPa升高到42.4 GPa，杨氏模量由92.0、62.7和57 GPa分别升高到95.6、102.2和97.5 GPa，并且3个方向上的杨氏模量差量最大值由46.7%到7.0%，各向异性较低．因此，壳聚糖组分能够对羟基磷灰石的力学性能起到明显的加强作用，且 HA-CS混合体系在各个方向的力学性能较单组分HA更为接近，从而可以克服因材料各向异性而导致的材料缺陷．综上可见，本研究可对羟基磷灰石和壳聚糖的复合结构的设计起到借鉴和理论指导作用．

基金项目：深圳市科技基础研究计划资助项目(JCYJ20150525092941007)；国家自然科学基金资助项目(21471102)

作者简介：袁秋华(1967—)，男，深圳大学副教授、博士．研究方向：无机仿生材料．E-mail：yuanqiuh@szu.edu.cn

引文：袁秋华，万 磊，石 鑫，等．HA-CS复合物界面作用与力学性能模拟[J]. 深圳大学学报理工版，2018，35(3)：299-306.

参考文献/References：

[1] HENCH L L. Bioceramics[J]. Journal of the American Ceramic Society, 1998,81(7):1705-1728.

[2] 巫剑波，秦草坪，袁秋华，等．湿化学法合成羟基磷灰石晶体及其表征[J]．深圳大学学报理工版，2015, 32(5)：465-472．

WU Jianbo, QIN Caoping, YUAN Qiuhua,et al.Synthesis and characterization of hydroxyapatite crystals using wet chemistry methods[J].Journal of Shenzhen University Science and Engineering,2015,32(5):465-472.(in Chinese)

[3] 袁秋华，巫剑波，秦草坪，等．羟基磷灰石-聚乳酸复合膜的制备与表征[J]．深圳大学学报理工版，2016, 33(1)：10-17．

YUAN Qiuhua, WU Jianbo, QIN Caoping,et al.Preparation and characterization of hydroxyapatite-polylactic acid HA-PLA composite film[J]. Journal of Shenzhen University Science and Engineering,2016,33(1):10-17.(in Chinese)

[4] MANO J F, SOUSA R A, BOESEL L F, et al. Bloinert, biodegradable and injectable polymeric matrix composites for hard tissue replacement: state of the art and recent developments[J]. Composites Science and Technology, 2004,64(6):789-817.

[5] DARDER M, ARANDA P, RUIZ-HITZKY E. Bionanocomposites: a new concept of ecological, bioinspired, and functional hybrid materials[J]. Advanced Materials, 2007, 19(10):1309-1319.

[6] 袁秋华，张自强，徐安平，等．羟基磷灰石-壳聚糖复合膜旋涂法制备及表征[J]．深圳大学学报理工版，2017, 34(5)：451-456．

YUAN Qiuhua, ZHANG Ziqiang, XU Anping,et al.Preparation and characterization of hydroxyapatite-chitosan composite coatings using a spin-coating technique[J].Journal of Shenzhen University Science and Engineering,2017,34(5):451-456.(in Chinese)

[7] LEE Y M, PARK Y J, LEE S J,et al.Tissue engineered bone formation using chitosan/tricalcium phosphate sponges[J].Journal of Periodontology,2000,71(3):410-417.

[8] PAL A K, PAL T K, MUKHERJEE K,et al.Animal experimentation with tooth derived Calcium hydroxyapatite based composites as bone-graft substitute biomaterials[J].Biomedical Sciences Instrumentation,1997,33(1):561-566.

[9] ZHANG Hongping, LU Xiong, LENG Yang,et al.Molecular dynamics simulations on the interaction between polymers and hydroxyapatite with and without coupling agents[J].Acta Biomaterialia,2009,5(4):1169-1181.

[10] SANDOVAL C, CASTRO C, GARGALLO L,et al.Specific interactions in blends containing Chitosan and functionalized polymers. Molecular dynamics simulations[J]. Polymer,2005,46(23):10437-10442.

[11] 汪焰恩，李鹏林，杨明明，等．快速成形中PAN基碳纤维/HA力学性能模拟研究[J]．机械工程学报，2013, 49(23)：6-11．

WANG Yanen, LI Penglin, YANG Mingming,et al.Molecular dynamics simulation to mechanical properties of Pan/HA in rapid prototyping[J].Journal of Mechanical Engineering,2013,49(23):6-11.(in Chinese)

[12] SUN H. Compass: an abintio force-field optimized for condensed phase application-overview with details on alkane and benzene compounds[J]. Journal of Physical Chemistry B, 1998, 102(38):7338-7364.

[13] 张红平．几种常规生物医用材料表面界面的计算机模拟[D]．成都：西南交通大学，2011.

ZHANG Hongping.Computer simulation of surface interface of several conventional biological medical materials[D].Chengdu:Southwest Jiaotong University,2011.

[14] LI Chunli, CHOI P. Molecular dynamics study of the adsorption behaviorof normal alkanes on a relaxed α-Al2O3(0001) surface[J]. Journal of Physical Chemistry C, 2007, 111(4):1747-1753.

[15] TUNG K L, LU K T, RUAAN R C,et al.Molecular dynamics study of the effect of solvent types on the dynamic properties of polymer chains in solution[J].Desalination,2006,192(1/2/3):380-390.

[16] PRATHAB B, SUBRAMANIAN V, AMINABHAVI T M. Molecular dynamicssimulations to investigate polymer-polymer and polymer-metal oxide interactions[J]. Polymer 2007, 48(1):409-416.

[17] FERNANDO A R, YOSADARA R M. Ab initio molecular dynamics calculations of the phase transformation mechanism for the formation of TiO2titanate-type nanosheets from anatase[J]. Chemistry of Materials, 2007, 19(12):2947-2959.

[18] TOTH R, FERRONE M, MIERTUS S, et al. Structure and energetics of biocompatible polymer nanocomposite systems: a molecular dynamics study[J].Biomacromolecules,2006,7(6):1714-1719.

[19] KAY M L, YOUNG R A. Crystal structure of hydroxylapatite[J]. Nature, 1964, 204(496):1050-1052.

[20] HAUPTMANN S, DUFNER H, BRICKMANN J, et al. Potential energy function for apatites[J]. Physical Chemistry Chemical Physics, 2003, 5:635-639.

[21] 魏庆华，汪焰恩，杨明明，等．HA/PLA复合材料界面相互作用及其力学性能的MD模拟[J]．功能材料，2013, 44(21)：3089-3093．

WEI Qinghua, WANG Yanen, YANG Mingming, et al.Molecula dynamics simulation of mechanical properties and surface interaction for HA/PLA[J]. Journal of Functional Maerials, 2013, 44(21):3089-3093.(in Chinese)

[22] BHOWMIK R, KATTI K S, KATTI D.Molecular dynamics simulation of hydroxyapatite-polyacrylic acid interfaces[J].Polymer,2007,48(2):664-674.