材料研究学报 2014 , 28 (2): 133-138 https://doi.org/10.11901/1005.3093.2013.325

羟基磷灰石/α-氰基丙烯酸正丁酯相互作用及力学性能的MD模拟^*

汪焰恩¹^,, 魏庆华¹^,², 杨明明¹^,², 魏生民¹^,²

1. 西北工业大学机电学院西安 710072
2. 西北工业大学现代设计与集成制造技术教育部重点实验室西安 710072

Molecular Dynamics Simulation of Mechanical Properties and Surface Interaction for HA/NBCA

WANG Yanen¹^,^**^,, WEI Qinghua¹^,², YANG Mingming¹^,², WEI Shengmin¹^,²

1. Mechatronics School, Northwestern Polytechnical University, Xi’an 710072
2. The Key Laboratory of Contemporary Design and Integrated Manufacturing Technology, Ministry of Education, Northwestern Polytechnical University, Xi’an 710072

中图分类号: TB321, O647

通讯作者: **To whom correspondence should be addressed, Tel: (029)88493232-329, E-mail:wangyanen@nwpu.edu.cn

收稿日期: 2013-05-25

修回日期: 2013-05-24

网络出版日期: --

基金资助: * 国家自然科学基金51175432、50905147, 教育部高校博士点专项科研基金20116102110046, 陕西省自然科学基金2011JQ7005, 西北工业大学博士论文创新基金201213, 西北工业大学研究生创业种子基金Z2013043 资助项目。

作者简介:

展开

摘要

基于羟基磷灰石陶瓷微球人工骨支架的制备工艺, 用分子动力学(MD)方法研究了羟基磷灰石(HA)的三个晶面(001)、(100)、(110)分别与生物胶黏剂α-氰基丙烯酸正丁酯(NBCA)相互作用后混合体系的结合能, 并计算分析了(110)晶面的力学性能和径向分布函数。结果表明, 三晶面所对应结合能关系为 $H A (110) > H A (100) > H A (001)$ ; HA(110)晶面与NBCA之间粘结强度高于HA(100)和HA(001)晶面; 对聚合度为40的NBCA聚合物与HA相互作用的MD计算表明, HA(110)/NBCA混合体系的力学性能比单组份HA体系有所下降, 但满足人工骨对力学性能的要求; 元素之间的径向分布函数揭示了混合体系HA/NBCA组分之间的相互作用机理, NBCA与HA(110)晶面存在较强的相互作用。其主要原因是, NBCA中的N原子和羰基中的O原子分别与HA中的H原子形成的氢键, 进而说明了HA对NBCA的强吸附作用。

关键词： 材料科学基础学科 ; 羟基磷灰石 ; α-氰基丙烯酸正丁酯 ; 分子动力学 ; 结合能 ; 力学性能 ; 径向分布函数

Abstract

Molecular dynamics simulation was applied to investigate the binding energy of NBCA on HA crystallographic planes (001), (100) and (110), and then the mechanical properties and radial distribution function of the HA(110)/NBCA mixed system were calculated and analyzed. The results show that HA (110) has the highest binding energy with NBCA because of its higher planar atom density than that of HA (001) and (100). The mechanical properties of HA(110)/NBCA mixed system is weaker than one-component HA system. However, the NBCA contains 40 monomers, its mechanical properties completely meet the requirements for the artificial bone scaffold. By calculating the radial distribution function of HA(110)/NBCA, the essence of its interface interaction were elucidated. There is a strong interaction between HA crystallographic plane (110) and NBCA, it mainly derives from the hydrogen bonds between O atoms which connect with C atoms of NBCA and H atoms in HA crystal, and a strong adsorption effect can be demonstrated between HA and NBCA.

Keywords： foundational discipline in materials science ; HA ; NBCA ; molecular dynamics ; binding energy ; mechanical properties ; radial distribution function

PDF (2756KB) 元数据多维度评价相关文章收藏文章

本文引用格式导出 EndNote Ris Bibtex

汪焰恩, 魏庆华, 杨明明, 魏生民. 羟基磷灰石/α-氰基丙烯酸正丁酯相互作用及力学性能的MD模拟^*[J]. , 2014, 28(2): 133-138 https://doi.org/10.11901/1005.3093.2013.325

WANG Yanen, WEI Qinghua, YANG Mingming, WEI Shengmin. Molecular Dynamics Simulation of Mechanical Properties and Surface Interaction for HA/NBCA[J]. 材料研究学报, 2014, 28(2): 133-138 https://doi.org/10.11901/1005.3093.2013.325

随着生物增材技术的应用, 人们用快速原形技术将生物陶瓷粉体、颗粒或者浆体材料直接制造人工骨支架产品, 以期减少支架的污染。羟基磷灰石(hydroxyapatite, HA)是人体骨组织的主要无机成分, 也是人体内骨和牙齿的主要无机成分, 是最适合做人体组织支架的生物活性材料^[1-3]。羟基磷灰石骨支架植入人体后有承力作用, 因此其力学性能是重要的评价指标。羟基磷灰石骨支架能与病患正常骨的弹性模量较好匹配, 可避免高弹性模量骨支架植入体内引起的应力遮挡效应^[4-5]。若弹性模量太小, 则承力作用相对弱化, 易导致二次骨折或骨裂。在制造骨支架过程中, 生物胶黏剂用量和成分对骨支架的力学性能有明显影响。因此有必要探索其作用机理, 以便探寻最优生产工艺参数。

制造人工骨支架所采用的基体材料羟基磷灰石和生物胶黏剂α-氰基丙烯酸正丁酯(α-n-butyl cyanoacrylate, NBCA)表面相互作用, 特别是在微观分子层面的机理, 有待深入研究。从微观角度研究HA与生物胶黏剂相互作用机理, 有助于深入了解人工骨支架微细结构的力学性能及其骨材料相互作用机制, 对宏观制造工艺参数有重要影响。

分子动力学(Molecular dynamics, MD)计算为材料微观性能研究提供了可行方法, 能揭示物质在分子层面相互作用机理。MD是一种基于原子力场的计算方法, 模拟不同物质之间的力学、电学等性质, 计算空间距离从几个埃到几个纳米, 以解释物质之间相互作用机理。Hong ping Zhang et al^[6]用MD模拟研究了在没有耦合剂作用条件下人工骨组成材料羟基磷灰石以及复合材料界面的相互作用, 比较分析了聚合物在三个不同平面的约束能, 对人工骨制造有一定的指导作用。潘海华等^[7]研究了人工骨材料羟基磷灰石(HA)晶面上的水分子行为, 发现HA晶面附近的结构化水层可阻止溶液离子自由进入晶面, 对HA颗粒在水溶液中的动力学稳定性有重要的影响。李敏等^[8]利用MD技术研究了HA及复合材料的制备与性能, 特别是人工骨的制备工艺及弹性模量。本文用MD方法探究人工骨支架制造过程中HA晶体表面与生物胶黏剂(α-氰基丙烯酸正丁酯)的界面相互作用, 计算它们之间的结合能以及力学特性。本文基于羟基磷灰石陶瓷微球人工骨支架的制备工艺, 用分子动力学(MD)方法研究了羟基磷灰石(HA)的三个晶面(001)、(100)、(110)分别与生物胶黏剂α-氰基丙烯酸正丁酯(NBCA)相互作用后混合体系的结合能, 并计算分析(110)晶面的力学性能和径向分布函数。

1 模拟方法与模型

1.1 力场

以Materials Studio(MS)为模拟平台, 采用原子模拟研究的凝聚态优化分子势场COMPASS(Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies, COMPASS)。它是第一个参数化可用于凝聚态物质研究的计算力场, 适用于有机材料与无机材料之间作用性能研究^[9]。COMPASS力场函数可表达为^[^10]:

$E_{p o t} = \sum^{} [K_{2} (b - b_{0})^{2} + K_{3} (b - b_{0})^{3} + (b - b_{0})^{4}]$ (a)

$+ \sum^{} H_{2} (θ - θ_{0})^{2} + H_{3} (θ - θ_{0})^{3} + H_{4} (θ - θ_{0})^{4}$ (b)

$+ \sum^{} υ_{1} [1 - c o s (Φ - Φ_{1}^{0})] + υ_{2} [1 - c o s (2 Φ - Φ_{2}^{0})]$

$+ υ_{3} [1 - c o s (3 Φ - Φ_{3}^{0})]$ (c)

$+ \sum^{} K_{χ} χ^{2}$ (d)

$+ \sum^{} \sum^{} F_{b b'} (b - b_{0}) (b' - b'_{0})$ (e)

$+ \sum^{} \sum^{} F_{θ θ'} (θ - θ_{0}) (θ' - θ'_{0})$ (f)

$+ \sum^{} \sum^{} F_{b θ} (b - b_{0}) (θ - θ_{0})$ (g)

$+ \sum^{} \sum^{} (b - b_{0}) [V_{1} c o s Φ + V_{2} c o s 2 Φ + V_{3} c o s 3 Φ]$ (h)

$+ \sum^{} \sum^{} (b' - b'_{0}) [V_{1} c o s Φ + V_{2} c o s 2 Φ + V_{3} c o s 3 Φ]$ (i)

$+ \sum^{} \sum^{} (θ - θ_{0}) [V_{1} c o s Φ + V_{2} c o s 2 Φ + V_{3} c o s 3 Φ]$ (j)

$+ \sum^{} \sum^{} \sum^{} K_{Φ θ θ'} c o s Φ (θ - θ_{0}) (θ' - θ'_{0})$ (k)

$+ \sum^{} \frac{q_{i} q_{j}}{ε γ_{i j}}$ (l)

$+ \sum^{} [\frac{A_{i j}}{γ_{_{i j}}^{9}} - \frac{B_{i j}}{γ_{_{i j}}^{6}}]$ (m)

方程(1)包含键伸缩项(a)、键角弯曲项(b)、双面角扭曲项(c)、离平面振动项(d)、交叉作用项(e-k)、静电作用项(l)和范德华力作用项(m)。 $K$ 、 $H$ 、 $F$ 、 $υ$ 、 $V$ 为力场参数; $b$ 、 $θ$ 、 $Φ$ 分别表示键长、弯曲角、扭转角。COMPASS力场能够准确地给出凝聚态体系的结构和各种性质^{[11, 12]}, 有效改善了以往力场难以清晰描述力场势能的缺陷^[13]。

1.2 HA及其表面模型

HA(化学式 $C a_{10} (P O_{4})_{6} {(O H)}_{2}$ )^[14]晶体结构的空间群为P63/m, 属六方晶系, 单位晶胞含有10个 $C a^{2 +}$ , 6个 $P O_{4}^{3 -}$ 和2个 $O H^{-}$ 。10个 $C a^{2 +}$ 占据2种位置, 4个 $C a^{2 +}$ 占 $Z = 0$ 和 $Z = 1 / 2$ (或0.5)的位置, 6个 $C a^{2 +}$ 中3个位于 $Z = 1 / 4$ (或0.25)、另3个位于 $Z = 3 / 4$ (或0.75) 的位置。 $O H^{-}$ 位于晶胞的4个角上, 6个 $P O_{4}^{3 -}$ 四配位体分别位于 $Z = 1 / 4$ (或0.25)和 $Z = 3 / 4$ (或0.75)的平面上, 这些 $P O_{4}^{3 -}$ 四面体的网络使得HA的结构具有较好的稳定性。晶格参数为: a=b=0.9432 nm, c=0.6881 nm, $α = β = 90^{0}$ , $γ = 120^{0}$ ^[15], 根据这些参数可构建HA的晶胞结构。

沿(0 0 1), (1 0 0), (1 1 0)三个剪切平面, 分别建立HA晶体表面模型, 将其非键截断距设置为1.55 nm, 用于与生物胶黏剂相互作用界面模拟的HA超级晶胞大小均为 $3 \times 2 \times 1$ 的规格, 如图1所示。

显示原图| 下载原图ZIP| 生成PPT

图1 HA表面的球棍模型, (a) HA(0 0 1), (b) HA(1 0 0), (c) HA(1 1 0)本论文中, 所有图内氧原子, 氢原子, 氮原子, 碳原子, 钙原子, 磷原子分别用红色, 灰白色, 深蓝色, 灰黑色, 绿色和粉红色表示

Fig.1 Ball and stick models of the HA surfaces, the dimension unit is angstrom, (a) HA(0 0 1), (b) HA(1 0 0), (c) HA(1 1 0) Color codes: hydrogen atom, white; calcium atom, green; oxygen atom, red; and phosphorus atom, purple

1.3 α-氰基丙烯酸正丁酯分子模型

α-氰基丙烯酸正丁酯(NBCA)是最常用瞬干胶的主要成分和有效成分, 其化学式为CH₂=C(CN)-COO-C₄H₉, 是最早用化学方法合成的医用胶, 也是临床应用最多的医用胶之一, 其主要成分是长链酯单体。在弱盐基性物质如水(湿气)、醇等环境下, 迅速发生阴离子聚合, 瞬间发挥其强粘结作用, 故本文用其作为制造人工骨支架的生物胶黏剂。

由于最终与HA发生界面相互作用的是NBCA高聚物, 本文使用MS软件平台对其进行建模, 构建聚合度为40的NBCA聚合物。在正则系综(NVT), 温度为298 K的环境下进行30飞秒的MD运算, 得到优化后的高分子立体三维化学结构, 如图2所示。

显示原图| 下载原图ZIP| 生成PPT

图2 聚合度为40的α-氰基丙烯酸正丁酯NBCA

Fig.2 Three-dimensional ball-and-stick models and the chemical structures of the NBCA chain for forty monomers

1.4 HA/NBCA相互作用

在图1所示的HA模型表面上分别放置聚合物NBCA链, 以构造出HA/NBCA模型, 同时在高聚物表面上方设置3 nm的真空层, 以便高分子聚合物仅作用于HA的一个侧表面。HA/NBCA相互作用的界面模型在正则系综(NVT)下进行分子动力学模拟计算, 其详细模拟计算参数设置如表1所示。

表1 HA/NBCA 分子动力学模拟的参数设置

Table 1 MD Simulation Parameters for HA/NBCA

Forcefield	Non- bond	Summation method	Ensemble	Temperature	Cut-off distance	Thermostat	Barostat	Equilibrium	Production	Time Step	Frame output
COMPASS	Vdw,Coulomb	Atom based, Ewald	NVT	310 K	1.55 nm	Andersen steps	Parrinello steps	300,000	100,000	0.1fs	Every 300 steps

新窗口打开

应用分子动力学模拟计算HA/NBCA的前、后力学模型分别表示为图3中的(a)、(b)、(c)。通过分子动力学模拟, 分析HA与NBCA的结合能、原子的态密度、杨氏模量以及径向分布等, 从而研究HA/NBCA的相互作用, 预测羟基磷灰石人工骨支架的力学性能, 以探究增材制造HA人工骨支架的工艺参数。

显示原图| 下载原图ZIP| 生成PPT

图3 HA(100)/(001)/(110)与NBCA混合体系分子动力学计算前/后模型

Fig.3 NBCA on the HA/(100)/(001)/ (110) surface before and after the MD simulation

2 结果及讨论

2.1 结合能

计算HA晶体表面和高聚物NBCA分子链之间的结合能, 以得到HA晶体表面和粘结剂NBCA之间的黏结力, 其计算公式为

(2) $E_{i n t e r} = E_{H A + N B C A} - (E_{\begin{matrix} H A & s u r f a c e \end{matrix}} + E_{N B C A})$

(3) $E_{b i n d i n g} = \frac{- E_{i n t e r}}{M} = \frac{E_{\begin{matrix} H A & s u r f a c e \end{matrix}} + E_{N B C A} - E_{H A + N B C A}}{M}$

其中 $E_{\begin{matrix} H A & s u r f a c e \end{matrix}}$ 表示HA的表面能量, $E_{N B C A}$ 为高聚物NBCA链的能量, $E_{H A + N B C A}$ 为HA/BNCA相互作用后混合体系的能量, M为高聚物NBCA所含原子个数, 该模拟计算体系为882个。通过MD模拟计算, 得到的相互作用能及对应的结合能如表2和图4所示(相互作用能见表2, 结合能见图4)。

表2 结合能能的计算 (compass)

Table 2 Binding energy calculation (compass)

	$E_{\begin{matrix} H A & s u r f a c e \end{matrix}}$ ( $K J ∙ m o l^{- 1}$ )	$E_{N B C A}$ ( $K J ∙ m o l^{- 1}$ )	$E_{H A + N B C A}$ ( $K J ∙ m o l^{- 1}$ )	- $E_{i n t e r}$ ( $K J ∙ m o l^{- 1}$ )
(0 0 1)	29485532.92	80290.54	29241730.82	324092.64
(1 0 0)	29780245.51	94178.91	29506537.85	367886.57
(1 1 0)	56216546.17	99364.56	55910123.82	405786.91

新窗口打开

显示原图| 下载原图ZIP| 生成PPT

图4 HA与NBCA的结合能柱状图

Fig.4 The binding energies between the NBCA and HA surfaces

从表2可以清楚地看出HA与NBCA在三个平面的表面能、体系能和NBCA链能值, 从而得到HA与NBCA在三个剪切面的结合能值。图4表示HA与NBCA在三个剪切面所对应结合能关系, $H A (110) > H A (100) > H A (001)$ , 结合能的大小反应了HA与NBCA之间粘结强度的大小, 结合能越大其粘结强度就越大^[16-17], 即HA(110)表面与NBCA的相互作用后的牢固性更强。高聚物NBCA与HA晶体表面结合能不同, 其原因与HA表面的原子密度D相关。原子密度的计算公式为

(4) $D = \frac{N}{A}$

其中N和A分别表示HA表面原子数和HA平面的表面积。对HA三个平面的Ca、 $P O_{4}$ 和 $O H$ 进行平面密度计算, 结果如表3所示。HA(110)的原子密度最大, 而其所对应的结合能也最大, 故原子密度的大小可以反应结合能的大小。即HA平面原子密度越大, 结合能也越大, 这也预示着HA(110)的界面作用最大。所以对HA/NBCA混合体系机械性能的研究, 只在HA(110)界面进行MD计算分析。

表3 HA表面的原子密度 (原子(×0.1 nm)^-2)

Table 3 The atom density of the HA surfaces (atom(×0.1 nm)^-2)

原子团	HA(001)	HA(110)	HA(100)
$C a$	0.02	0.23	0.21
$P O_{4}$	0.11	0.12	0.11
$O H$	0.06	0.072	0.072

新窗口打开

2.2 力学性能

为了探究HA/NBCA混合体系的加工工艺参数和机械性能, 分析经过分子动力学运算达到平衡后的HA/NBCA无定形结构, 研究其混合体系的力学性能。在Material Studio平台上用Forcite 模块中的Mechanical Properties功能对其进行MD计算分析。通过计算即可获得HA(110)/NBCA相互作用体系的杨氏模量、泊松比等力学参数(其计算结果列于表4)。

表4 杨氏模量E与泊松比λ

Table 4 Young's modulus and Poisson's ratio

Axis	Young Modulus/GPa		Poission Ratios
	HA	HA/NBCA	Poission Ratios
X	54.3	50.5	Exy	0.28	Exz	0.29
Y	54.3	51.4	Eyx	0.34	Eyz	0.31
Z	62.5	38.2	Ezx	0.21	Ezy	0.26

新窗口打开

表4列出了HA(110)/NBCA混合模型的力学参数。可以看出, 该混合模型在三个方向的弹性模量并不一致。这与羟基磷灰石的各向异性相符, 即使添加胶黏剂NBCA并不改变材料的各向异性。相互作用后材料的杨氏模量E比单一HA的杨氏模量E在Z方向上有所减小, 但完全能够满足人工骨对材料的力学性能3~20 GPa的最低要求, 该混合体系的泊松比也与自然骨的泊松比较为接近, λ约为0.3左右^[18]。

2.3 径向分布函数——界面相互作用的本质

径向分布函数 g(r)又称对相关函数(pair correlation function), 是反映材料微观结构的特征物理量。它表示在一个分子周围半径为r距离的范围处, 出现另一个分子的概率密度相对于随机分布概率密度的比值。对整个平衡轨迹进行径向分布函数g(r)分析, 由径向分布函数揭示出 HA(110)/NBCA界面相互作用的方式和本质。

NBCA分子中的原子类型有C原子、H原子、N原子、羰基中的O原子和羰基相连的O原子。其中主链C原子周围有H原子等环绕, 难以与HA分子界面产生相互作用, 故只考虑其余3类原子, 依次表示为N、Oa₂和Oa₁。类似地, HA中含有Ca原子、P原子、磷酸根中与P相连的O原子、H原子和羟基中的O原子, 依次表示为Ca、P、O_b1、H_b、O_b2。本文计算了NBCA中3类原子与HA中5类原子之间的g(r)图, 图5(a～f)给出有重要相互作用的g(r)～r关系。

显示原图| 下载原图ZIP| 生成PPT

图5 HA(110)/NBCA的径向分布函数

Fig.5 The radial distribution function of HA(110)/NBCA

分子间的径向分布函数g(r), 能揭示非键原子之间的相互作用方式和本质。分子间相互作用分为氢键(r＝0.15～0.31 nm)和vdw作用(r＝0.31～0.5 nm)^[19], 大于0.5 nm的vdw作用很微弱。图5(a, b): $N - H_{b}$ 、 $O_{a 2} - H_{b}$ 径向分布函数均在r=0.2 nm附近取得最高峰值, 比O、H 原子的共价半径之和0.11 nm略大, 表明NBCA中的N原子和羰基中的O原子分别与HA中的H原子形成了氢键; 而图5(c-f): $N - C a$ 、 $N - O_{b 1}$ 、 $O_{a 2} - C a$ 、 $O_{a 2} - O_{b 1}$ 的径向分布函数在r=0.2~0.25 nm处取得第一峰值, 对应于以分子状态吸附在HA表面的第一水合层。这验证了HA的表面对NBCA具有强大的吸附性, 其中 $N - C a$ 、 $O_{a 2} - C a$ 的最高峰值分别在r=0.215 nm、r=0.21 nm处取得, 比他们所对应的离子键略短, 表明HA中的Ca原子与高聚物中的N原子、 $O_{a 2}$ 原子形成了离子键; 因此, 由对g(r)的分析表明, HA(110)/NBCA界面相互作用主要源于NBCA中N原子和羰基中的O原子分别与HA中的H原子之间较强的氢键以及 $N - C a$ 、 $O_{a 2} - C a$ 之间的离子键, 表明HA对NBCA具有强大的吸附性。

3 结论

1. 高聚物NBCA与HA晶体表面结合能的大小关系 $H A (110) > H A (100) > H A (001)$ , 与表面原子密度相关。表面原子密度越大, 其结合能也越大, 即原子密度的大小可以反应结合能的大小。

2. 聚合度为40的NBCA聚合物与HA相互作用后其混合体系的力学性能较单组份HA的力学性能略有降低, 但完全满足人工骨对全力学性能的要求。

3. NBCA在HA(110)晶面存在较强的相互作用, 主要源自于NBCA中的N原子和羰基中的O原子分别与HA中的H原子形成的氢键。HA对NBCA具有很强的吸附性, 与HA材料具有较强吸附性相符。

〈

〉

羟基磷灰石/α-氰基丙烯酸正丁酯相互作用及力学性能的MD模拟*

Molecular Dynamics Simulation of Mechanical Properties and Surface Interaction for HA/NBCA

1 模拟方法与模型

1.1 力场

1.2 HA及其表面模型

1.3 α-氰基丙烯酸正丁酯分子模型

1.4 HA/NBCA相互作用

2 结果及讨论

2.1 结合能

2.2 力学性能

2.3 径向分布函数——界面相互作用的本质

3 结论

羟基磷灰石/α-氰基丙烯酸正丁酯相互作用及力学性能的MD模拟^*