骨陷窝--骨细胞形状和方向对骨单元内液体流动行为的影响1)

于纬伦 武晓刚,2) 李朝鑫 孙玉琴 张美珍 陈维毅,3)

∗(太原理工大学生物医学工程学院，太原 030024)

†(太原理工大学体育学院，太原 030024)

引言

骨组织是一种非常智能的材料，它能够根据日常所受的生理载荷调整其骨量和微观结构.这种自适应的行为被认为是由骨细胞主导的.骨细胞是骨内所有细胞中对力学刺激最敏感的细胞[1]，且负责将力学刺激信号转化为成骨细胞和破骨细胞的调节信号，从而使骨适应环境变化[2].实验研究指出负荷引起的流体流动在骨细胞膜上产生剪切应力是骨重建过程中生化反应被激活的主要的刺激[3-5]，且与骨的力传导机制有关[6].因此，研究骨细胞在骨陷窝--骨小管系统内经历的流体流动行为的特征是非常重要的.

由于骨细胞嵌入矿化的细胞外基质，直接进行实验研究具有难度，而体外实验并不能真实呈现骨细胞在体内所受的力学环境[7].研究人员试图通过建立骨骼间质流体流动的数学模型来克服这个问题[8].渗透率是影响间质流体通过骨内孔隙流动的一个重要特性.在利用数值模拟方法来探索外力引起的骨内液体流动时，研究人员通过多孔弹性介质模型计算出骨内的液体流速、压力梯度、流体剪切力等参数，但这些模型都假定骨的渗透性是各向同性的(骨细胞为球型)[7-9]，不能准确反映骨内具体微观结构的(如骨细胞--骨小管系统的形状，方向和密度等)液体流动信息.

近年来，人们对骨细胞的空间特性(包括密度和形态)以及这些特性与疾病和衰老的潜在关系越来越感兴趣.体外实验研究表明，骨细胞几何形状影响其应变响应[10].Carter 等发现股骨的前侧，后侧，内侧和外侧的骨细胞的密度，形状和方向都有较大差别，其可能原因是与负载的局部变化有关[11].骨质疏松症已经成为世界范围内的一个严重问题，最近有研究表明患有骨质疏松的病人的骨细胞形状更不规则[12-13]，骨小管更弯曲[13]，而其骨细胞膜上感受的流体剪应力和流体速度也发生较大的变化[14].年龄也是影响骨细胞形状的一个重要因素，年龄增大会导致骨细胞面积减小且形状偏扁平[15-16].以前的研究一般都认为骨陷窝--骨细胞的方向与骨中胶原纤维的方向一致[9,17-19]，即渗透率主方向与坐标轴方向一致.事实上,二者方向会有不一致的情况[11,20-21]，这会使得骨内液体流动主方向(即渗透率主方向)与胶原纤维的方向成一定的夹角.骨陷窝--骨细胞--骨小管的这种形态和方向变化会导致骨组织内渗透率的各向异性行为.在进行有限元多孔弹性力学分析时，需要准确量化骨陷窝--小管孔隙的渗透性来捕捉骨的各向异性流体流动行为，才能更真实地预测孔隙流体压力和流速响应.

因此，为了更精确地阐明骨内液体流动的具体形式，本研究开发一种基于骨组织微观结构的多孔弹性力学有限元模型.首先根据骨陷窝的形状，方向和其他微观结构生理参数，计算得到骨的三维渗透率.然后，基于多孔弹性力学理论建立骨单元的有限元模型，计算在轴向载荷作用下骨单元内液体的压力和流速，以期为骨的力传导与骨的功能适应性机制提供一个更深入的理解.

1 基于骨单元微观结构计算其渗透率

1.1 骨单元渗透率的计算

图1(a)是整段骨组织，图1(b)是图1(a)中的骨单元.假设图1(b)中的骨陷窝--小管系统的排列规则且骨小管分布均匀，则骨单元可以看成是由一个个骨基质包裹着的骨陷窝--骨细胞--小管系统的正方体的周期单元(CPUC)组成(图1(c)).图1(d)是骨小管的微观结构，其中rc是骨小管半径，ro是骨细胞突触半径，a0是骨细胞突触周围纤维基质的半径.

图1 骨组织分层结构示意图Fig.1 The hierarchical structure for bone tissue

根据骨陷窝--骨细胞--骨小管的微观结构得出其渗透率klcp[22-23]

其中，q是rc(0.23µm)和ro(0.1µm)之间的无量纲比(q=rc/ro)[7]，γ 是rc与纤维基质填充的单个骨小管的渗透性(kp)平方根的无量纲比γ=，其中kp=(∆/a0)2.377[17,22-23]，a0的半径是5 nm，∆是纤维与纤维的间距(7 nm)[23-24].

A1和B1由以下方程得出

其中，I0，K0，I1和K1是第1 类和第2 类的修正贝塞尔函数.

L表示两个骨陷窝之间的距离，也是CPUC 的边长，可由以下公式得出

式中，VL为单位体积，NLac为单位体积内骨陷窝的数量.NLac的范围在2.6 ∼9.0×104(NLac/mm3)[5,11,25-27]，本文选取单位体积的NLac值为3.7×104[28]，得到L的值为30µm.根据文献报告的平均测量值，每个骨陷窝周围的骨小管总数N=62[28].通过对几种不同生物的骨细胞的组织形态进行考察发现骨陷窝的形状类似于椭球形状[18].用椭球的标准方程来表示骨陷窝

其中，a,b,c分别为骨陷窝长半轴、中半轴和短半轴.则骨陷窝--小管的孔隙率ϕ 为

式中，Lc表示骨小管的平均长度.值得注意的是由于骨细胞与骨陷窝空间内的液体会进行液体交换，且与骨陷窝空间的渗透率相差不大[7,29]，所以我们把整个骨陷窝空间和骨细胞体都算作孔隙空间.ni(i=x,y,z) 是骨小管穿过CPUC 某个面的数目.如图2(a)所示，为了确定x，y和z方向的nx，ny和nz，我们假设骨小管基于骨陷窝椭圆体投影表面积三维分布，并可以用骨陷窝的投影面积比来测量不同方向骨小管的数量[17]

其 中，nx，ny和nz分别是平行于x，y和z轴穿过CUPC 各面的小管数；Sx，Sy和Sz分别是x，y和z方向上骨陷窝的投影表面积，N是一个骨陷窝周围的骨小管总数.如图2(b)所示，当骨陷窝绕x轴偏转θ 角时，骨陷窝在，和轴上的投影面积和骨小管的三维分布同样可由式(6)求出.

图2 (a)骨陷窝沿着x 轴，y 轴和z 轴的投影面积(S x，S y 和S z).(b)骨陷窝绕x 轴旋转θ 角Fig.2 (a)Projected surface areas of the lacuna-osteocyte in the x,y,and z directions(S x,S y and S z)were used to determine the distribution of canaliculi;(b)The lacuna-osteocyte rotates around x-axis by θ

1.2 骨陷窝形状

为了更清晰描述骨陷窝形状，首先规定了3 个特征值EV1，EV2 和EV3(EV1 是长半轴的平方，EV2 是中半轴的平方和EV3 是短半轴的平方).然后，从颗粒形状的研究中修改了3 种指标，平衡度(degree of equancy)(EV3:EV1)，伸长度(degree of elongation)(1–EV2:EV1)和扁平度(flatness)(1–EV3:EV2)[11,30]，以定义形状的差异程度.如表1 所示，取几组不同形状的骨陷窝来观察骨陷窝形状的差异对其内部液体流动的影响.图3 展示了这3 种指标对骨陷窝形状影响的重要性.

图3 骨陷窝形状的代表性示例.骨陷窝在不同平面的投影形状(a) x-y 面，(b)y-z 面，(c)x-z 面Fig.3 Representative cases of lacunae shapes.The project of lacunae shapes are shown schematically in(a)the x-y plane,(b)the y-z plane,and(c)the x-z plane

当3 个特征值相等(EV1=EV2=EV3)时，骨陷窝的形状即为球形(reference).当前两个特征值相等，第3 个特征值较小(EV1=EV2>EV3)时，骨陷窝的形状呈扁平形(Case 1 和Case 2)；第1 个特征值较大，后两个特征值相等(EV1>EV2=EV3)时，骨陷窝的形状呈伸长形(Case 3 和Case 4)；3 个特征值都不相同(EV1>EV2>EV3)，形状偏向取决于伸长度和扁平度的大小(Case 5 和Case 6).

表1 骨陷窝代表性示例的几何参数和指标Fig.1 Geometrical and degree of representative cases

根据表1 中骨陷窝形状参数计算得到Case1-Case 6 的渗透率分别为K1，K2，K3，K4，K5和K6

参考模型(reference)模型的渗透率是各向同性的其值计算得：1.05×10−20m2.

1.3 骨陷窝方向

在以前的研究中通常假设骨组织孔隙流体的渗流速度方向与坐标轴方向相同.事实上,由于骨陷窝方向的不同会带动渗透率方向的偏转，如果还按原来的方法，则计算结果会产生较大误差.骨陷窝的轴长一般在几微米到十几微米之间[2,7,11,14,31]，当如图2(b)所示的骨陷窝(a=15µm，b=9µm，c=3.9µm)发生θ 角的偏转时，骨组织渗透率可以用全方向的渗透率张量表示

当θ 角分别为0◦，30◦，45◦，60◦和90◦时,对应于o-xyz坐标下渗透率分别为

2 骨单元控制方程及有限元模型的建立

2.1 控制方程

本文将骨单元看成由CUPC 单元组成的固--液耦合多孔弹性材料，其本构方程可由以下方程来描述[33-35]

其中，σ 是整个多孔介质对应的应力张量，ε 为应变张量.M是多孔介质脱水的4 阶弹性矩阵，α 和M分别是Biot 有效应力系数张量和模量.p为多孔介质孔隙内液体压力，ξ 为单位体积内液体含量的改变量.tr()表示矩阵的迹算子.

在不考虑体积力的情况下的平衡方程[8]

液体质量守恒方程

达西定律可写为[36]

式中，V是流速矢量，k和µ分别是渗透张量和液体的黏度.将式(8)代入式(10)，式(9)和式(12)代入式(11)得到骨单元的控制方程

2.2 有限元模型的建立

在这项工作中，使用COMSOL Multiphysics 软件探讨轴向压缩载荷作用下骨单元内的流体与固体相互作用的多孔弹性行为.如图4(a)，数值模拟过程中骨单元模型被定义为一个由CPUC 组成的横观各向同性的中空的柱体(忽略了哈弗氏管)，其部分材料和几何参数设置如表2 所示.

表2 几何和材料常数Table 2 Geometrical and material constants used in osteon model

考虑到载荷频率较低，哈弗氏管像蓄水池一样起到储蓄的作用，能够维持流体正常的流进流出，其内的压力可设置为0，用作参考压力.哈弗氏管的孔隙尺寸远大于骨陷窝--小管孔隙，骨单元壁的孔隙压力可以由此释放.哈弗氏管表面的应力可以忽略不计，这也导致在哈弗氏管处应力为0 的边界条件[37].骨单元外壁设置刚性约束(位移为0)，表示骨单元受周围骨间质约束不能移动[33].如图4(a)和方程(14)，轴向对称的循环载荷代表纵向压缩的生理运动，其谐波幅值(w)和频率(f)分别为0.5µm 和1 Hz.在一个周期内的最大轴向应变(ε)发生在0.5 s 时，其幅度为0.001[33].

图4 骨单元模型的建立(a)和网格划分(b)Fig.4 Establishment of osteon model(a)and mesh generation(b)

如图4(b)所示，借助扫掠的网格剖分方法来生成精确有效的有限元网格，网格包含20 880 个单元.最小单元质量为0.601 5；平均单元质量为0.847 4.

3 结果

3.1 骨陷窝形状的影响

根据本文的加载形式，当t=0.25 s 时，应变率达到最大值0.001 s−1，骨单元内孔隙液体压力和流速的变化规律有最大响应[12,33].如图5 绘制了t=0.25 s时考虑不同形状骨陷窝CUPC 的骨单元孔隙压力及流速的分布图.在轴向的对称载荷下，参考模型呈现的是对称分布与之前的研究结果吻合[12,33,36].由于内部骨陷窝形状的不同，Case 1 ∼Case 6 的液体的压力和流速(靠近哈弗氏管部位较为明显)的大小与分布规律均与参考模型有较大差别，且最大的孔隙压力值均大于参考模型.这说明在骨陷窝体积相同的情况下，形状的差异程度会影响压力和流速的大小和分布规律.Case 1，Case 2 和Case 5 的骨陷窝形状都是扁平形的，Case 3,Case 4 和Case 6 的骨陷窝形状都是伸长形的.Case 1，Case 2 和Case 5 的压力和流速在分布上与参考模型有些类似，沿着y−z面和x−z面的分布差异不大.Case 3，Case 4 和Case 6 的模型的压力和流速在y−z面和x−z面的分布差异较大.

由图6 可以看出，在t=0.25 s 时，Case 1--参考模型的压力和流速分别为3.21 × 105Pa，4.29 ×105Pa，3.9×105Pa，3.43×105Pa，3.57×105Pa，3.84× 105Pa，2.51 × 105Pa 和6.34 × 10−8m/s 和6.57 ×10−8m/s,6.5 × 10−8m/s，6.81 × 10−8m/s，6.74 ×10−8m/s，6.4 × 10−8m/s，6.72 × 10−8m/s 和5.85 ×10−8m/s.在不同Case 模型中，Case 2 和Case 4 模型沿y轴的压力差距最大，增大了86%，Case 5 和参考模型沿y轴流速差距最大，增大了18%.在单个Case 模型上，Case 1，Case 2，Case 5 和参考模型沿x轴和y轴的最大压力和流速差别不大，而Case 3，Case 4 和Case 6 沿x轴和y轴的压力有明显区别，分别增大了62%，47%和58%，而流速差别不大.

图5 骨陷窝形状影响下的骨单元孔隙压力(p)及流速(v)的分布图Fig.5 The pore pressure(p)and flow velocity(v)distribution of osteon due to the influence of lacunae shape

图6 骨陷窝形状对骨单元孔隙压力和流速沿着不同方向的影响Fig.6 Effect of the lacunae shape on the pore pressure and flow velocity along the different directions of osteon

3.2 骨陷窝方向的影响

图7 绘制了t=0.25 s 时考虑具有不同方向骨陷窝CUPC 的骨单元孔隙压力及流速的分布图.在我们所研究参数范围内，0◦∼90◦模型的液体的压力和流速的大小与分布规律均有较大差别.随着骨陷窝旋转角度的增加其沿着y−z面和x−z面的压力分布差距越大.30◦，45◦，60◦模型的流速在哈弗氏管处出现流速局部增大的情况.

图7 骨陷窝方向影响下的骨单元孔隙压力(p)及流速(v)的分布图Fig.7 The pore pressure(p)and flow velocity(v)distribution of osteon due to the influence of lacunae direction

由图8 可以看出，骨陷窝在0◦∼90◦的旋转过程中导致骨单元的模型的压力值逐渐减小.在y轴方向上30◦，45◦和60◦的模型流速出现明显的不均匀分布.渗透率的主值方向与坐标轴方向出现偏转，即液体流动主方向(,和方向)与骨胶原纤维方向(x,y和z) 成θ 角时会导致骨组织内压力和流速大小和分布规律的变化.根据所选择的骨陷窝(a=15µm，b=9µm，c=3.9µm)及其偏转的角度(绕x轴偏转0◦∼90◦)，不同模型之间的最大的压力和流速有着较大差异，其中0◦沿y轴的压力比90◦增大了125%；30◦沿y轴的流速比90◦增大了56%.同一模型中由于偏转角度的不同，沿不同坐标轴方向最大压力和流速也存在这较大差异，其中60◦模型x轴上的最大压力比y轴上增大了58%;30◦模型沿y轴上的流速比x轴上增大了50%.可见，骨陷窝角度的不同，可明显使得局部的压力和流速发生改变.

图8 骨陷窝方向对骨单元孔隙压力和流速沿着不同方向的影响Fig.8 Effects of the lacunae direction on the pore pressure and flow velocity along the different directions of osteon

图8 骨陷窝方向对骨单元孔隙压力和流速沿着不同方向的影响(续)Fig.8 Effects of the lacunae direction on the pore pressure and flow velocity along the different directions of osteon(continued)

4 讨论

骨组织的渗透率是骨陷窝--小管流体流动空间的渗透率，不能脱离流体流动通道独立存在，也不可简单矢量合成.本研究首先通过估计不同形状和方向的骨陷窝中骨小管的数量，其次利用这些骨小管的三维分布以及骨组织其他微结构参数(骨陷窝密度，骨小管半径和骨细胞突触的半径等)来计算骨组织的渗透率和孔隙率等参数，最后根据计算所得的参数建立骨单元的多孔弹性有限元模型.结果表明，由于骨陷窝形状和方向的不同对骨单元内液体压力和流速的大小和分布规律有着明显的影响.

孔隙流体压力是骨陷窝--小管系统中的一种重要负荷诱导现象，它在很大程度上干预骨细胞生长，分化和化学物质运输[7,38].骨陷窝的形状和方向发生变化时，孔隙压力行为有着明显的改变.具体来说，对于形状改变的模型，Case 3，Case 4 和Case 6 在x和y方向上的压力差距最大，而Case1 Case2，Case5 和参考模型差距不大.这种情况和骨陷窝形状的改变引起的渗透率的各向异性有关，可以发现，Case 3，Case 4 和Case 6 模型在x和y方向的渗透率差了一个数量级，而Case1 Case2，Case5 和参考模型在x和y方向的渗透率在一个数量级且差距不大.而在z方向上所有的模型的压力都没有明显的变化，这是因为负载诱导产生流体压力的机制使得液体通过骨陷窝--小管流入哈弗氏管，使孔隙压力得到释放，即骨单元内液体的主要流动是在骨单元外壁和哈弗氏管之间，而在z方向上几乎没有流动.

对于角度不同的模型，设定的是骨陷窝绕x轴旋转θ 角.实际上，两个坐标系中x轴与x轴是重合的，因此骨单元在x方向上的渗透率是相同的.如图8(a)和图8(c)可以看到，沿x轴上的压力变化远小于沿y轴的变化.在y方向的渗透率变化大约在一个数量级左右.如图8(c)，在y轴方向上0◦∼90◦的模型的渗透率是逐渐增大的，其压力是逐渐减小的.在30◦，45◦和60◦时，骨中胶原纤维的方向与骨陷窝的方向不在一个方向，会导致y-z方向的渗流行为，引起流速局部增大的情况.总体来看，压力大小的变化与渗透率有很大相关性，随着渗透率的减小，压力有增大的趋势.

虽然骨陷窝方向的改变，使得30◦，45◦和60◦模型局部的流速发生较大的变化，但就整体来看骨陷窝形状和方向的改变对流速的影响都不超过一个数量级.这种现象可以解释为，渗透性由1.0×10−20m2降低到1.0×10−21m2的过程中抑制了流体流动，但孔隙压力和压力梯度相应的增加维持了流体速度的大小，从而维持了骨陷窝--小管内正常的对流运输，因此骨内达西流速变化不大.

在典型的人类日常生理活动的载荷(1.0×103µε)和频率(1 Hz)下[39]，分别计算出了各向异性骨陷窝--小管孔隙的最大孔隙压力和流速响应，压力范围从2.6×105Pa 到4.8×105Pa，流速范围从4.7×10−8m/s 到9.8×10−8m/s.压力的大小足以驱动流体跨皮质流动(大于6.0×103Pa)[7].流速大小也足以直接刺激骨细胞(大于2.0×10−8m/s)[36].每一个骨单元(哈弗氏系统)都可以充当其自身的排水系统，并通过哈弗氏管释放因骨组织变形产生的流体压力.尽管单个骨单元在整体骨组织中受到骨外膜到骨内膜产生的跨皮质压力梯度的影响，但一个哈弗氏系统产生的局部的压力梯度完全可以控制本身骨陷窝--骨小管流体的流动.因此，本文未考虑相邻骨单元之间的耦合程度，选用不透水的粘合线计算流体压力和速度大小是符合生理值要求的.

我们研究的一个局限性是骨小管被理想化为直管.实际上，骨细胞突触是通过弯曲的小管从骨陷窝延伸到CUPC 表面[40].这项研究没有考虑到骨小管弯曲的影响，因为精确计算弯曲的骨小管非常困难.这个限制可以在模型的下一步改进中解决.另一个局限性是本文把骨单元看成是由一个个相同的CUPC 组成，但是实际上每个CUPC 内的骨陷窝形状和方向会有区别，这将引起局部的液体流动的变化.虽然理论上模型需要确定每个CUPC 内骨陷窝的形状和方向，但据实验观察了解到，在骨组织的一定区域内骨陷窝的形状是类似的[11]，且骨中胶原纤维的方向与骨陷窝的方向有一定的相关性[18-19].这样的区域大到包含一个或者几个骨单元大小[11]，因此本研究的结果可用于建立骨单元的各向异性多孔弹性有限元模型.因此，只要确定所分析的特定区域内骨陷窝形状和方向，就可以把本研究中计算的各向异性渗透率应用到有限元模型中，以求解载荷诱导的流体压力和速度场.

5 结论

本研究提出了一种根据骨陷窝形状和方向估计骨单元内骨陷窝--骨小管孔隙各向异性渗透性的方法，并通过确定的骨陷窝--小管渗透率来描述骨单元内液体各向异性的流动.结果表明，骨陷窝的形状和方向以及骨小管的三维分布情况是决定骨陷窝--小管孔隙内液体流动的各向异性程度的重要参数.形状不同的模型最大压力和流速的差异分别可达86%和18%；方向不同的模型最大压力和流速的差异分别可达125%和56%.本文的研究对进一步理解骨力传导机制和一些像骨质疏松之类的骨质疾病具有重要参考价值.