研究人体颈椎及相关疾病生物力学机制,创建一个解剖精细、非线性的正常人全颈椎(C0~T1)三维有限元模型。模型基于1例女性健康志愿者颈椎CT数据建立,采用MIMICS13.1、Hypermesh11.0、Abaqus 6.12-1等有限元软件依次创建、预处理、运算和分析。在生理静态载荷下分别模拟颈椎活动(前屈、后伸、侧弯和旋转),观察应力集中部位;测量相邻椎体相对活动度(ROM)。本模型经文献中ROM结果验证可靠;预测生理载荷下枕骨大孔前部及侧部是上颈椎的应力集中点;中下段颈椎的应力集中点大部分为活动方向对侧的椎弓根和小关节突。全颈椎(C0~T1)非线性有限元模型的建立,可以为深入了解颈椎及其相关疾病的生物力学机制提供更理想的理论研究平台。
引用本文: 王辉昊, 詹红生, 陈博, 石印玉, 李玲慧, 杜国庆. 正常人全颈椎(C0~T1)三维有限元模型的建立与验证. 生物医学工程学杂志, 2014, 31(6): 1238-1242,1249. doi: 10.7507/1001-5515.20140235 复制
引言
颈椎是整个脊柱最容易发生损伤的区域,多根重要神经和血管经颈椎穿行至颅,颈椎区域一旦受损,危害极大,甚至威胁生命,因此选择合适的脊柱生物力学模型是探索脊柱损伤诊断和治疗机制的重要途径[1-3]。不同种类的模型中,计算机模型可同时反映外部运动行为和内部力学响应,获得体内或体外实验研究难以实现或精确测量的结果,也可进一步为研究提供指导性建议[4]。
目前颈椎有限元模型大致分为动态模型和静态模型两类[5],前者主要侧重于多个椎体在载荷下总体反应的有效预测,为了减少计算时间,该类模型的几何图像、连接重建和材料分配等常被简化;后者更合适研究复杂载荷条件下内部应力、应变和其他生物力学反应,因此更为精细[6-7]。
以前的研究多集中在某个功能节段,并为其创建非线性模型,如寰枕复合体等[8-9]。本研究的目标则是通过高分辨率CT图像改进几何图形准确性和材料属性的精细分配,创建正常人非线性全颈椎(C0~T1)三维有限元模型。研究通过静态加载生理载荷,模拟正常颈椎活动(前屈、后伸、左右侧弯和左右旋转),测量相邻椎体间的相对活动度(range of motion,ROM),建立经验证后确定可靠的颈椎有限元模型,为深入了解颈椎及其相关疾病生物力学机制提供理想的理论研究平台。
1 材料与方法
1.1 样本采集
招募1例健康志愿者(女性,46岁,身高165 cm,体重63 kg),问诊无症状主诉,动静态触诊及X线检查排除颈椎明显退变、脊柱感染、骨折、肿瘤、结核、严重脊柱畸形、严重骨质疏松、强直性脊柱炎、Paget病和颈椎外创史等情况。自愿参加本研究,签署知情同意书。中国注册临床试验伦理审查委员会批准号:(ChiECRCT-2013009)。
1.2 数据获取
上海中医药大学附属曙光医院放射科提供的GE Light Speed VCT 64层螺旋CT扫描及断层图像灰度处理。受试者仰卧位,颈肩背部自然放松,保持扫描断面与身体垂直。扫描枕骨底C0上2 mm至T1下2 mm获得体层图像。条件:140 kV,200 mA,层厚0.625 mm。分辨率512×512像素,获得385幅二维断层图像,调整图像灰度级对比度。
1.3 初步固体模型建立
首先将受试者CT数据导入比利时Materialise公司的交互式医学影像控制系统(MIMICS13.1软件),基于阈值分析进行边界拾取,阈值界定为226~3 071 Hu。经区域增长、编辑蒙罩、腔隙填补、边界划分等功能进行分割、修补,以清晰显示各骨骼,计算生成颈椎3D模型,即三角网格模型,以.stl格式数据保存;将其导入逆向工程软件Geomagic8.0软件,经点阶段、多边形阶段、成形阶段,对图像进行修补、去噪、铺面处理,转化为NURBS曲面模型。
1.4 划分网格
曲面模型导入有限元前处理软件Hypermesh11.0,经拓扑分区及网格划分,将网格质量Jacobian比控制在0.6以上。皮质骨采用平均厚度为1 mm的C3D6单元,松质骨采用C3D4单元,小关节软骨和终板采用0.1 mm的C3D6单元。椎间盘(含髓核和纤维环)采用增强沙漏控制的C3D8R单元,忽略骨髓影响。所有关节(寰枕关节、寰枢关节、寰椎前弓和齿状突关节、齿状突与横韧带关节以及小关节突关节)的关节面均定义为滑动接触关系,摩擦系数为0.1。
横韧带采用M3D4R单元,厚度为0.5 mm;其余韧带采用SPRINGA弹簧单元。根据文献资料[10-11],设置12种关键韧带起止点:包括寰枕前膜(anterior atlanto-occipital membrane,AAOM)、寰枕后膜(posterior atlanto-occipital membrane,PAOM)、十字韧带垂直部(cruciate ligamentum vertical portion,CLV)、齿突尖韧带(alar ligament,AP)、翼状韧带(apical ligament,AL)、覆膜(membranae tectoria,TM)、前纵韧带(ligamenta longitudinale anterius,ALL)、后纵韧带(ligamenta longitudinale posterius,PLL)、黄韧带(ligamentum Flavum,FL)、关节囊韧带(capsular ligament,CL)、棘突间韧带(interspinous ligament,ISL)、棘上韧带(supraspinal ligament,SSL)。除横韧带(ligamentum transversum,TL)外,所有韧带模型采用非线性面面通用接触关系模拟关节间的相互作用。
1.5 材料赋值
骨骼(皮质骨、松质骨)、终板和TL采用黏弹性材料。椎间盘(髓核、纤维环)和小关节面软骨采用不可压缩的超弹材料。各材料基于应变能理论Mooney-Rivilin超弹性材料公式进行运算:W=C10(I1-3)+C01(I2-3)+1/d(J-2)2(其中C10、C01为材料常数,I1、I2为应力张量的第1、第2不变量)。本研究不涉及椎间盘内部应力应变运算,故忽略胶原纤维的影响。各组织结构的材料属性(密度、弹性模量和泊松比)参数依据参考文献设置[5, 12-17]。
韧带(除TL)采用弹塑性材料,忽略材料塑性区和破坏区,中性区采用抛物线方程拟合,直线段采用线性方程拟合,中性区终点(dn,fn)和失效点定义弹性区的斜率为(df,ff) [12]。预设中性区为失效变形公式:dn/df=1/3;中性区高度为失效载荷公式:fn/ff=1/10。韧带设置参数依据参考文献设置[18-19]。
1.6 载荷与边界条件设置
根据研究目的,约束T1下面的6个自由度作为边界条件。在颅底旋转轴上方选择一参考点,建立C0上表面所有单元节点与此参考点的Distribution Coupling(参考点上的受力情况会换算成均布载荷施加于C0所有的从节点上)。对参考点施加纯扭矩,活动方向参考x、y、z全局坐标(x-y平面为水平面、x-z平面为冠状面、y-z平面为矢状面),屈伸方向与冠状面平行,旋转方向参考颈曲切线方向,侧弯时垂直于颈曲切线方向并与矢状面平行。
2 结果
2.1 非线性全颈椎(C0~T1)有限元模型
本模型共包含171 433个节点,563 879个单元,77个组件,5种单元类型,8种材料,27种材料属性。模型清晰完整地从不同角度观察椎体解剖结构,能多彩色、透明或任意组合显示,较真实地反映出各结构空间位置关系;体网格形状规整,整体分布合理,结构清楚 (见图 1)。

(a)全颈椎(C0~T1)椎体及椎间盘模型; (b)颈椎关键韧带; (c)横韧带、小关节软骨和终板
Figure1. Non-linear FE model of cervical spine (C0-T1)(a) Features of the finite element model of the lower part of the occipital (C0) and 8 vertebrae (C1-T1) and 6 intervertebral discs (C2-T1); (b) Features of the finite element model of 12 major ligaments; (c) Features of the finite element model of TL and 16 pairs of facet joint cartilages (C0-T1) and 6 pairs of endplates (C2-T1)
2.2 全颈椎有限元模型的运算与验证
2.2.1 模型的运算结果
采用有限元软件Abaqus 6.12-1进行运算:在1.5 Nm扭矩下颈椎模型分别模拟前屈、后伸、左右侧屈和轴向旋转的6种工况。图像左侧彩虹条图为von Mises等效应力变化数值(分为24级,范围:0~1.598×102),左下方为三维坐标信息,下部为模型信息) (见图 2)。

(a)前屈活动:59.9°; (b)后伸活动:43.8°; (c)左侧弯活动:40.2°; (d)右侧弯活动:38.6°; (e)左旋活动:59.5°; (f)右旋活动:57.5°
Figure2. von Mises stress of cervical spine movements in 6 degrees of freedom unde 1.5 Nm moment (Unit: MPa)(a) flexion: 59.9°; (b) extension: 43.8°; (c) leftbending: 40.2°; (d) rightbending: 38.6°; (e) lefttorsion: 59.5°; (f) righttorsion: 57.5
2.2.2 有效性验证
通过测算模型的ROM与文献结果进行对比。ROM计算方法:屈伸和侧屈:测量各椎体上切迹斜率,计算各椎体活动角度,获得相邻椎体活动角度差;旋转:将棘突与椎体后缘中点连线,计算各椎体旋转活动角度。相同设置下模型ROM与文献[19]、[9]的数据比较,如图 3、4所示。
由图 3可见,两个模型在各活动度中C0~C2的ROM角度差异较大,C2~C3、C3~C4、C4~C5、C5~C6差异较小,C7~T1未比较。
由图 4可见,两个模型在屈伸和侧弯活动中C0~C2,旋转的C1~C2的ROM角度无明显差异,旋转时C0~C1略有差别。
3 讨论
本研究根据受试者CT数据经MIMICS13.1、Hypermesh11.0、Abaqus 6.12-1等有限元软件依次创建、预处理、运算和分析,建立了一个解剖精细、验证可靠(经文献验证)的非线性全颈椎有限元模型。模型在解剖上包括全部椎体、椎间盘、终板、小关节软骨和主要韧带,较完整地涵盖了颈椎生理结构;材料上涉及多种材料,物理属性参照人体标本,较全面地反映了颈椎的真实状态;功能上能模拟人体正常生理活动。
结果显示,枕骨髁是颈椎活动时的上颈椎的应力集中点;齿状突与其前部的C1前弓和后部的横韧带的相互支持作用提供了C1~C2主要稳定性,这与文献[9]结果一致,也符合体外实验结果[20]。中下段颈椎的应力集中点大部分为活动方向对侧的椎弓根和小关节突,其中C4左侧下关节突在后伸及左右侧弯时应力异常集中,可能提示受试者存在潜在病理状态。模型各节段的ROM与文献[19]建立的C0~C7模型相似度较高,变化趋势一致,但C0~C2节段ROM相对略小,这可能与文献采用的颈椎标本节段有关,C0~C2部位软组织结构复杂,体内试验中软组织在解冻液体作用下,产生更大活动度。因此,本模型C0~C2节段与同为非线性模型的文献[9]的结果更为接近。
有限元技术在现代骨科生物力学分析中获得了不断进展,但在中医骨伤科学领域尚有很大进步空间。例如,常用于治疗脊柱慢性筋骨病损的矫正“骨错缝”手法等,虽然具有明显的临床优势[21-22],但由于生物力学理论薄弱,其安全性及其作用机制受到一定程度的质疑[23-24]。现代生物力学原理和医学图像三维有限元分析的结合,可以较全面地反映颈椎及相关病理状态下的生物力学特性。这些基础研究数据可能为手法诊治的定位定量分析,关键技术规范、治疗方案优化和手法安全性的提高提供重要的试验依据[25]。本研究参照文献设置,分别对模型施加1.0 Nm和1.5 Nm两种低扭矩载荷模仿生理活动,基本可以满足颈椎相关疾病及中医骨伤手法诊治的生物力学研究需要。
虽然理论上的有限元模型能够模拟脊柱的多种物理状态,但实际应用还需要谨慎对待下列问题:① 根据研究目的建立模型:颈椎周围的肌肉、韧带、椎间盘和椎体理化成分复杂,组织特性非均一性,而且肌肉体外实验数据较为少见[26],完全模仿真实颈椎目前难以实现。因此,研究者需要根据不同目的进行假设、调制、简化组建模式,建立一个设计合理且符合实际需要的力学模型。② 颈部血管与血流研究:颈部血管是颈椎重要组成部分,通过骨组织、椎间盘、韧带和血管以及血流的真实动态模拟,分析外力对颈椎结构内部力学效应对血管壁与血流之间的影响,可以准确地展现颈椎血管相关疾病的生物力学特征[27]。因此,动脉血管与周围组织相互作用的流固耦合分析可能是未来有限元模型的研究方向之一。
引言
颈椎是整个脊柱最容易发生损伤的区域,多根重要神经和血管经颈椎穿行至颅,颈椎区域一旦受损,危害极大,甚至威胁生命,因此选择合适的脊柱生物力学模型是探索脊柱损伤诊断和治疗机制的重要途径[1-3]。不同种类的模型中,计算机模型可同时反映外部运动行为和内部力学响应,获得体内或体外实验研究难以实现或精确测量的结果,也可进一步为研究提供指导性建议[4]。
目前颈椎有限元模型大致分为动态模型和静态模型两类[5],前者主要侧重于多个椎体在载荷下总体反应的有效预测,为了减少计算时间,该类模型的几何图像、连接重建和材料分配等常被简化;后者更合适研究复杂载荷条件下内部应力、应变和其他生物力学反应,因此更为精细[6-7]。
以前的研究多集中在某个功能节段,并为其创建非线性模型,如寰枕复合体等[8-9]。本研究的目标则是通过高分辨率CT图像改进几何图形准确性和材料属性的精细分配,创建正常人非线性全颈椎(C0~T1)三维有限元模型。研究通过静态加载生理载荷,模拟正常颈椎活动(前屈、后伸、左右侧弯和左右旋转),测量相邻椎体间的相对活动度(range of motion,ROM),建立经验证后确定可靠的颈椎有限元模型,为深入了解颈椎及其相关疾病生物力学机制提供理想的理论研究平台。
1 材料与方法
1.1 样本采集
招募1例健康志愿者(女性,46岁,身高165 cm,体重63 kg),问诊无症状主诉,动静态触诊及X线检查排除颈椎明显退变、脊柱感染、骨折、肿瘤、结核、严重脊柱畸形、严重骨质疏松、强直性脊柱炎、Paget病和颈椎外创史等情况。自愿参加本研究,签署知情同意书。中国注册临床试验伦理审查委员会批准号:(ChiECRCT-2013009)。
1.2 数据获取
上海中医药大学附属曙光医院放射科提供的GE Light Speed VCT 64层螺旋CT扫描及断层图像灰度处理。受试者仰卧位,颈肩背部自然放松,保持扫描断面与身体垂直。扫描枕骨底C0上2 mm至T1下2 mm获得体层图像。条件:140 kV,200 mA,层厚0.625 mm。分辨率512×512像素,获得385幅二维断层图像,调整图像灰度级对比度。
1.3 初步固体模型建立
首先将受试者CT数据导入比利时Materialise公司的交互式医学影像控制系统(MIMICS13.1软件),基于阈值分析进行边界拾取,阈值界定为226~3 071 Hu。经区域增长、编辑蒙罩、腔隙填补、边界划分等功能进行分割、修补,以清晰显示各骨骼,计算生成颈椎3D模型,即三角网格模型,以.stl格式数据保存;将其导入逆向工程软件Geomagic8.0软件,经点阶段、多边形阶段、成形阶段,对图像进行修补、去噪、铺面处理,转化为NURBS曲面模型。
1.4 划分网格
曲面模型导入有限元前处理软件Hypermesh11.0,经拓扑分区及网格划分,将网格质量Jacobian比控制在0.6以上。皮质骨采用平均厚度为1 mm的C3D6单元,松质骨采用C3D4单元,小关节软骨和终板采用0.1 mm的C3D6单元。椎间盘(含髓核和纤维环)采用增强沙漏控制的C3D8R单元,忽略骨髓影响。所有关节(寰枕关节、寰枢关节、寰椎前弓和齿状突关节、齿状突与横韧带关节以及小关节突关节)的关节面均定义为滑动接触关系,摩擦系数为0.1。
横韧带采用M3D4R单元,厚度为0.5 mm;其余韧带采用SPRINGA弹簧单元。根据文献资料[10-11],设置12种关键韧带起止点:包括寰枕前膜(anterior atlanto-occipital membrane,AAOM)、寰枕后膜(posterior atlanto-occipital membrane,PAOM)、十字韧带垂直部(cruciate ligamentum vertical portion,CLV)、齿突尖韧带(alar ligament,AP)、翼状韧带(apical ligament,AL)、覆膜(membranae tectoria,TM)、前纵韧带(ligamenta longitudinale anterius,ALL)、后纵韧带(ligamenta longitudinale posterius,PLL)、黄韧带(ligamentum Flavum,FL)、关节囊韧带(capsular ligament,CL)、棘突间韧带(interspinous ligament,ISL)、棘上韧带(supraspinal ligament,SSL)。除横韧带(ligamentum transversum,TL)外,所有韧带模型采用非线性面面通用接触关系模拟关节间的相互作用。
1.5 材料赋值
骨骼(皮质骨、松质骨)、终板和TL采用黏弹性材料。椎间盘(髓核、纤维环)和小关节面软骨采用不可压缩的超弹材料。各材料基于应变能理论Mooney-Rivilin超弹性材料公式进行运算:W=C10(I1-3)+C01(I2-3)+1/d(J-2)2(其中C10、C01为材料常数,I1、I2为应力张量的第1、第2不变量)。本研究不涉及椎间盘内部应力应变运算,故忽略胶原纤维的影响。各组织结构的材料属性(密度、弹性模量和泊松比)参数依据参考文献设置[5, 12-17]。
韧带(除TL)采用弹塑性材料,忽略材料塑性区和破坏区,中性区采用抛物线方程拟合,直线段采用线性方程拟合,中性区终点(dn,fn)和失效点定义弹性区的斜率为(df,ff) [12]。预设中性区为失效变形公式:dn/df=1/3;中性区高度为失效载荷公式:fn/ff=1/10。韧带设置参数依据参考文献设置[18-19]。
1.6 载荷与边界条件设置
根据研究目的,约束T1下面的6个自由度作为边界条件。在颅底旋转轴上方选择一参考点,建立C0上表面所有单元节点与此参考点的Distribution Coupling(参考点上的受力情况会换算成均布载荷施加于C0所有的从节点上)。对参考点施加纯扭矩,活动方向参考x、y、z全局坐标(x-y平面为水平面、x-z平面为冠状面、y-z平面为矢状面),屈伸方向与冠状面平行,旋转方向参考颈曲切线方向,侧弯时垂直于颈曲切线方向并与矢状面平行。
2 结果
2.1 非线性全颈椎(C0~T1)有限元模型
本模型共包含171 433个节点,563 879个单元,77个组件,5种单元类型,8种材料,27种材料属性。模型清晰完整地从不同角度观察椎体解剖结构,能多彩色、透明或任意组合显示,较真实地反映出各结构空间位置关系;体网格形状规整,整体分布合理,结构清楚 (见图 1)。

(a)全颈椎(C0~T1)椎体及椎间盘模型; (b)颈椎关键韧带; (c)横韧带、小关节软骨和终板
Figure1. Non-linear FE model of cervical spine (C0-T1)(a) Features of the finite element model of the lower part of the occipital (C0) and 8 vertebrae (C1-T1) and 6 intervertebral discs (C2-T1); (b) Features of the finite element model of 12 major ligaments; (c) Features of the finite element model of TL and 16 pairs of facet joint cartilages (C0-T1) and 6 pairs of endplates (C2-T1)
2.2 全颈椎有限元模型的运算与验证
2.2.1 模型的运算结果
采用有限元软件Abaqus 6.12-1进行运算:在1.5 Nm扭矩下颈椎模型分别模拟前屈、后伸、左右侧屈和轴向旋转的6种工况。图像左侧彩虹条图为von Mises等效应力变化数值(分为24级,范围:0~1.598×102),左下方为三维坐标信息,下部为模型信息) (见图 2)。

(a)前屈活动:59.9°; (b)后伸活动:43.8°; (c)左侧弯活动:40.2°; (d)右侧弯活动:38.6°; (e)左旋活动:59.5°; (f)右旋活动:57.5°
Figure2. von Mises stress of cervical spine movements in 6 degrees of freedom unde 1.5 Nm moment (Unit: MPa)(a) flexion: 59.9°; (b) extension: 43.8°; (c) leftbending: 40.2°; (d) rightbending: 38.6°; (e) lefttorsion: 59.5°; (f) righttorsion: 57.5
2.2.2 有效性验证
通过测算模型的ROM与文献结果进行对比。ROM计算方法:屈伸和侧屈:测量各椎体上切迹斜率,计算各椎体活动角度,获得相邻椎体活动角度差;旋转:将棘突与椎体后缘中点连线,计算各椎体旋转活动角度。相同设置下模型ROM与文献[19]、[9]的数据比较,如图 3、4所示。
由图 3可见,两个模型在各活动度中C0~C2的ROM角度差异较大,C2~C3、C3~C4、C4~C5、C5~C6差异较小,C7~T1未比较。
由图 4可见,两个模型在屈伸和侧弯活动中C0~C2,旋转的C1~C2的ROM角度无明显差异,旋转时C0~C1略有差别。
3 讨论
本研究根据受试者CT数据经MIMICS13.1、Hypermesh11.0、Abaqus 6.12-1等有限元软件依次创建、预处理、运算和分析,建立了一个解剖精细、验证可靠(经文献验证)的非线性全颈椎有限元模型。模型在解剖上包括全部椎体、椎间盘、终板、小关节软骨和主要韧带,较完整地涵盖了颈椎生理结构;材料上涉及多种材料,物理属性参照人体标本,较全面地反映了颈椎的真实状态;功能上能模拟人体正常生理活动。
结果显示,枕骨髁是颈椎活动时的上颈椎的应力集中点;齿状突与其前部的C1前弓和后部的横韧带的相互支持作用提供了C1~C2主要稳定性,这与文献[9]结果一致,也符合体外实验结果[20]。中下段颈椎的应力集中点大部分为活动方向对侧的椎弓根和小关节突,其中C4左侧下关节突在后伸及左右侧弯时应力异常集中,可能提示受试者存在潜在病理状态。模型各节段的ROM与文献[19]建立的C0~C7模型相似度较高,变化趋势一致,但C0~C2节段ROM相对略小,这可能与文献采用的颈椎标本节段有关,C0~C2部位软组织结构复杂,体内试验中软组织在解冻液体作用下,产生更大活动度。因此,本模型C0~C2节段与同为非线性模型的文献[9]的结果更为接近。
有限元技术在现代骨科生物力学分析中获得了不断进展,但在中医骨伤科学领域尚有很大进步空间。例如,常用于治疗脊柱慢性筋骨病损的矫正“骨错缝”手法等,虽然具有明显的临床优势[21-22],但由于生物力学理论薄弱,其安全性及其作用机制受到一定程度的质疑[23-24]。现代生物力学原理和医学图像三维有限元分析的结合,可以较全面地反映颈椎及相关病理状态下的生物力学特性。这些基础研究数据可能为手法诊治的定位定量分析,关键技术规范、治疗方案优化和手法安全性的提高提供重要的试验依据[25]。本研究参照文献设置,分别对模型施加1.0 Nm和1.5 Nm两种低扭矩载荷模仿生理活动,基本可以满足颈椎相关疾病及中医骨伤手法诊治的生物力学研究需要。
虽然理论上的有限元模型能够模拟脊柱的多种物理状态,但实际应用还需要谨慎对待下列问题:① 根据研究目的建立模型:颈椎周围的肌肉、韧带、椎间盘和椎体理化成分复杂,组织特性非均一性,而且肌肉体外实验数据较为少见[26],完全模仿真实颈椎目前难以实现。因此,研究者需要根据不同目的进行假设、调制、简化组建模式,建立一个设计合理且符合实际需要的力学模型。② 颈部血管与血流研究:颈部血管是颈椎重要组成部分,通过骨组织、椎间盘、韧带和血管以及血流的真实动态模拟,分析外力对颈椎结构内部力学效应对血管壁与血流之间的影响,可以准确地展现颈椎血管相关疾病的生物力学特征[27]。因此,动脉血管与周围组织相互作用的流固耦合分析可能是未来有限元模型的研究方向之一。