1 引言

静脉激光治疗(EVLT)是一种静脉曲张的微创治疗方法,光纤通过导管被送至静脉曲张的血管腔内,光纤末端出射的激光直接作用于周围组织,激光的选择性光热效应使血液凝固,形成血栓,并间接损伤静脉壁,从而使得曲张静脉纤维化闭合。

1999年,西班牙静脉学家Boné首先报道了应用EVLT治疗下肢静脉曲张^[1];2002年,血管腔内激光治疗获美国食品药品管理局(FDA)正式批准^[2]。EVLT通过皮肤穿刺造成的损伤远小于传统静脉剥落手术造成的损伤,静脉关闭率高,并发症少,恢复快,已成为治疗静脉曲张最有前途的技术之一^[3]。我国于2001年在上海率先开展了EVLT^[4]。

成功的EVLT依赖于最佳治疗参数的选择,如何达到有效的血管损伤,同时避免并发症,是EVLT最值得深究的问题之一。激光辐照静脉血管腔内,光穿透区域的血管及周围组织直接吸收激光能量^[5],实现对组织的有益治疗或破坏。光穿透区域的大小受组织中散射的影响,因此波长的选择要综合考虑散射和吸收这两方面的因素。目前,810,940,980,1064,1320 nm波长的激光器已用于EVLT^[6-10],鉴于水在1470~1500 nm范围内的吸收系数较大,近年来该波段的半导体激光器开始被引入到EVLT中^[11-12];但是连续激光和脉冲激光的选择使得辐照剂量过大,从而引起疼痛或感觉异常,导致出现静脉曲张复发等诸多问题,因此仍有大量研究试图寻找最佳的辐照剂量和波长^[13-14],以期达到血管闭塞和纤维化,同时减少并发症。

目前对EVLT的大部分研究集中在临床试验上,对其数值模拟的研究很少。Mordon等^[15-16]构建了一个二维血管模型,采用有限元方法比较了810 nm和980 nm脉冲模式和连续模式下静脉血管的损伤,并采用扩散近似的方法模拟了980 nm和1320 nm激光辐照下静脉血管的损伤;龚玮^[17]结合蒙特卡罗算法和有限差分法建立了1064 nm激光光分布和热损伤的数学模型;van Ruijven等^[18]通过求解生物传热方程计算了血管周围的温度分布。现有的血管数学模型均是平面多层结构,未考虑血管的管状结构对光子分布的影响。

建立EVLT的光学模型可以更好地了解激光在血管及其周围组织中的传输过程,有助于选择最佳治疗参数。本文构建了静脉曲张血管及其周围组织的管状光学模型,采用三维蒙特卡罗算法对不同波长激光辐照静脉腔内时静脉血管及其周围组织的光分布进行了数值模拟,探讨了波长、血管管径、辐照位置、光学特性参数等对治疗效果的影响。

2 基本原理

2.1 三维蒙特卡罗算法

采用解析法难以直接求解组织内的光分布。蒙特卡罗模拟作为一种数值计算方法,尤其适用于复杂结构组织的研究^[19-20],可以将其结果视为实际物理测量的低噪声描述^[21]。

三维蒙特卡罗算法在经典马尔可夫链蒙特卡罗(MCML)算法的基础上,创建三维笛卡尔网格,同时给每个网格分配一个组织类型,每种组织类型在特定波长下具有相应的光学特性参数(包括吸收系数μ_a,散射系数μ_s,各向异性因子g)。因此对于非均匀组织,可以根据组织结构,设置各个网格相应的组织类型。

图1所示为建模所需的球坐标系示意图,移动球坐标系中z轴的方向动态跟随光子的传输方向,p₁(x₁,y₁,z₁)为光子自原点O沿光子入射方向移动抽样步长s₁后的笛卡尔坐标,当前移动球坐标系的z轴与笛卡尔坐标系中的z轴重合。之后,抽样p₁处的散射偏向角θ₁和方位角φ₁,根据方向余弦公式(见表1)确定下一步的传输方向。光子在该方向上移动抽样步长s₂到达p₂(x₂,y₂,z₂),以此类推。ξ为(0,1)之间均匀分布的随机数,上述步长s的抽样值为

s=-lnξμa+μs;(1)

偏向角θ∈[0,π]的抽样值为

θ=arccos12g1+g2-1-g21-g+2gξ2,(g≠0);(2)

方位角φ∈[0,2π]的抽样值为

φ=2πξ。(3)

图 1. 移动球坐标系示意图

Fig. 1. Moving diagram of spherical coordinate

下载图片 查看所有图片

模拟过程中,定义三维数组W(x,y,z)储存各网格的权重,光子的初始权重为1,当光子与组织相互作用时,光子损失一部分权重,损失的权重储存在与光子发生相互作用的网格中,计入W(x,y,z)。光子每一步权重、位置、方向余弦的变化如表1所示,其中ω为光子当前的权重,(x,y,z)为光子当前的笛卡尔坐标,(μ_x,μ_y,μ_z)为光子当前的方向余弦,ω'为光子移动步长s后的权重,(x',y',z')为光子移动步长s后的笛卡尔坐标,(μ_x',μ_y',μ_z')为光子移动步长s后的方向余弦。光子到达该步的最后一个网格后,会根据当前网格的光学特性参数确定下一步的散射方向。当光子权重小于阈值权重时,由轮盘赌决定是否终止光子。

2.2 静脉血管建模

静脉血管及其周围组织的三维模型如图2(a)所示,从中心向外依次为血液、血管壁、静脉周围组织,它们的光学特性参数^[3]如表2所示,g取0.9^[22],所选波长均已在文献[ 3,6-18]中用于EVLT。本研究中长度单位均采用cm,μ_a和μ_s的单位均采用cm^-1。静脉曲张血管壁内径d₁=0.253 cm,血管壁厚度d₂=0.14 cm^[23],组织体模型大小为0.6 cm×0.6 cm×0.6 cm,分割网格个数为300×300×300,每个网格大小为0.002 cm×0.002 cm×0.002 cm,光束的总能量均为1 J,发射光子数为10⁷。

光源采用360°径向辐射式光源^[24],如图2(b)所示,在仿真过程中,将其简化为各向同性点光源,光源坐标为(0,0,0.3),简化模型如图3所示。

图 2. (a)静脉血管及其周围组织的三维模型;(b) 360°辐射式光源

Fig. 2. (a) 3D model of vein and perivenous tissues; (b) 360° radial light source

下载图片 查看所有图片

图 3. 静脉血管简化模型。(a) xz平面(y=0);(b) yz平面(x=0)

Fig. 3. Simplified model of vein. (a) xz plane(y=0); (b) yz plane(x=0)

下载图片 查看所有图片

3 结果与讨论

3.1 波长对光能流率的影响

采用蒙特卡罗方法分别模拟波长为810,940,980,1320,1500 nm的激光辐照静脉腔内时血管及其周围组织的光分布,即在上述不同波长辐照下,xz平面(y=0)的光分布,结果如图4所示,其中Φ为光能流率。图5所示为相应的纵向光能流率分布,x=0,y=0,z∈[0.3,0.6],见图3(a)中虚线处。因为简化模型中的血液、血管壁同轴,且光源位于轴上,光能流率呈柱状对称分布,所以图5仅绘出一侧的纵向光能流率分布。

图 4. 不同波长下xz平面(y=0)的光分布。(a) 810 nm;(b) 940 nm;(c) 980 nm;(d) 1320 nm;(e) 1500 nm

Fig. 4. Light distributions of xz plane (y=0) at different wavelengths. (a) 810 nm; (b) 940 nm; (c) 980 nm; (d) 1320 nm; (e) 1500 nm

下载图片 查看所有图片

图 5. 不同波长时纵向光能流率分布

Fig. 5. Longitudinal distribution of fluence with different wavelengths

下载图片 查看所有图片

波长为810,940,980 nm的激光辐照静脉腔内时,血管及其周围组织的光分布相似,纵向光能流率曲线几乎重合。在这3个波长下,血液及血管壁的吸收系数较小,部分光子会进入静脉周围组织,在实际手术中会造成对静脉周围组织的损伤。而1320 nm的光子分布相对集中,小部分光子进入静脉周围组织;1500 nm的光子几乎全部被血液及血管壁吸收,在同等剂量的情况下,对周围组织的损伤最小。各波长在不同组织中的光学特性参数不同,导致光的分布不同,在选择激光治疗剂量时,需要考虑光学特性参数的影响,应选择合适的波长,这样有利于增加靶组织吸收,减轻周围组织的损伤,减少并发症。

3.2 血管直径对光能流率的影响

发生静脉曲张的部位不同,病情轻重不同,血管的直径差异很大,而管径对EVLT的光分布影响很大。在3.1节比较同一管径下各波长光分布差异的基础上,选择对周围组织损伤较小的1500,1320 nm激光进行研究,比较其辐照血管壁内径变为0.6d₁、1.2d₁、1.8d₁、2.4d₁以及血管壁厚度变为0.6d₂、1.2d₂、1.8d₂、2.4d₂时的畸形血管,分析它们的光分布差异,结果如图6所示。由于血管管径差异较大,因此模拟参数修改如下:组织体大小为1.4 cm×1.4 cm×1.4 cm,分割网格个数为350×350×350,每个网格大小为0.004 cm×0.004 cm×0.004 cm,其余参数不变。

图6(a)~(d)所示为选择1500 nm激光时,不同管径xz平面(y=0)的光分布,相应的纵向光能流率分布如图7(a)所示。其中0.6d₁和0.6d₂时的纵向光能流率与其他管径的差异较为明显。由于血管较细,因此只有少量光子能到达静脉周围组织。在图7(a)中,当z>0.86 cm(z=0.86 cm为血管壁与周围组织交界处)后,0.6d₁和0.6d₂时的光能流率衰减梯度明显减小,这是由于血管周围组织的μ_a和μ_s均小于血管壁的μ_a和μ_b,因此衰减系数较小,而1.2、1.8、2.4倍管径的光子仅分布在血液中,衰减系数未出现明显改变,不能对血管壁造成足够损伤,无法起到闭合静脉的作用。

图6(e)~(h)所示为选择1320 nm激光,不同管径时xz平面(y=0)的光分布,相应的纵向光能流率分布如图7(b)所示。在1320 nm激光的辐照下,这四种管径的血管周围组织中均有光子进入,管径越小的血管,其周围组织中的光能流率越大,损伤越大。2.4倍管径时,光子基本分布在血液和血管壁中,只有少量光子能进入周围组织中。

图 6. 不同波长、血管直径时xz平面(y=0)上的光分布。(a) 1500 nm,0.6d1,0.6d2;(b) 1500 nm,1.2d1,1.2d2;(c) 1500 nm,1.8d1,1.8d2;(d) 1500 nm,2.4d1,2.4d2;(e) 1320 nm,0.6d1,0.6d2;(f) 1320 nm,1.2d1,1.2d2;(g) 1320 nm,1.8d1,1.8d2;(h) 1320 nm,2.4d1,2.4d2

Fig. 6. Light distribution of xz plane (y=0) with different wavelengths and vessel diameters. (a) 1500 nm, 0.6d1,0.6d2; (b) 1500 nm, 1.2d1, 1.2d2; (c) 1500 nm, 1.8d1, 1.8d2; (d) 1500 nm, 2.4d1, 2.4d2; (e) 1320 nm,0.6d1, 0.6d2; (f) 1320 nm, 1.2d1, 1.2d2; (g) 1320 nm, 1.8d1, 1.8d2; (h) 1320 nm, 2.4d1, 2.4d2

下载图片 查看所有图片

图 7. 不同波长、血管直径时的纵向光能流率分布。(a) 1500 nm;(b) 1320 nm

Fig. 7. Longitudinal distribution of fluence with different wavelengths and vessel diameters. (a) 1500 nm; (b) 1320 nm

下载图片 查看所有图片

由此可知,在进行静脉腔内激光治疗时,应根据畸形血管的粗细选择合适的波长,光能流率分布可以作为参考。较细的血管(1倍管径以下)适合选择1500 nm的激光,而较粗的血管适合选择1320 nm的激光,以减轻对周围组织的损伤。

3.3 光源位置对光能流率的影响

在进行激光腔内治疗时,光源位置无法始终保持在静脉中轴线上,因此,比较光源位置在z轴上分别上移0.02,0.04,0.06,0.08,0.10 cm时静脉血管及其周围组织的光分布差异,模拟采用的参数同3.1节,结果如图8所示,相应的轴向光能流率如图9(a)所示。比较光源上移0.02 cm和0.10 cm时的纵向光能流率可知:在近光源端血管壁与周围组织交界处(x=0,z=0.0335 cm),后者的光能流率是前者的50多倍;而在远光源端血管壁与周围组织交界处(x=0,z=0.5665 cm),后者的纵向光能流率已减至0 J/cm²。若光源偏离中轴线较远,靠近一侧的血管壁时,靠近光源的这一侧血管壁的光能流率会远大于另一侧的,从而导致靠近光源一侧血管壁出现溃疡和穿孔;而远离光源一侧的血管壁光能流率较小,无法达到损伤血管的目的。为了避免纤维末端与血管壁直接接触,Vuylsteke等^[25]提出在光纤外套一个超弹性材料制成的具备形状记忆的导管,如图9(b)所示。

图 8. 不同光源位置下xz平面(y=0)的光分布。(a) (0,0,0.30);(b) (0,0,0.28);(c) (0,0,0.26);(d) (0,0,0.24);(e) (0,0,0.22);(f) (0,0,0.20)

Fig. 8. Light distributions of xz plane (y=0) at different light source positions. (a) (0,0,0.30); (b) (0,0,0.28);(c) (0,0,0.26); (d) (0,0,0.24); (e) (0,0,0.22); (f) (0,0,0.20)

下载图片 查看所有图片

图 9. (a)不同光源位置的纵向光能流率分布;(b) Vuylsteke等[25]研制的导管

Fig. 9. (a) Longitudinal distribution of fluence at different light source positions; (b) catheter developed by Vuylsteke et al[25]

下载图片 查看所有图片

3.4 血液光学特性参数对光能流率的影响

对于身体不同部位静脉血管中的血液,其血氧饱和度、血细胞比容、静脉周围组织结构及生色团均存在差异,因此血液、血管壁、静脉周围组织的光学特性参数也不相同。以血细胞比容为例^[26],在同一波长、剂量的激光辐照下,血细胞比容分别为0.84%和42.1%的血液,其μ_a和μ_s相差数十倍。由于大部分光子在血液中被吸收,且血液的光学特性参数容易受其他因素影响,因此比较1320 nm激光辐照下,血液光学特性参数μ_a和μ_s分别改变为原参数的1/20时,静脉血管及其周围组织的光分布。模拟采用的血液分别编号为血液1、血液2和血液3,它们的μ_a分别为3.8,0.19,3.8 cm^-1,μ_s分别为54,54,2.7 cm^-1。血液1采用表2中1320 nm时的原始参数;血液2的其余各光学特性参数保持不变,μ_a改变为原参数的1/20;血液3的μ_s改变为原参数的1/20,其余各光学特性参数保持不变。管径改变为原参数的2.4倍,模拟参数如下:组织体大小为1.4 cm×1.4 cm×1.4 cm,分割网格个数为350×350×350,每个网格大小为0.004 cm×0.004 cm×0.004 cm。

图10所示为1320 nm激光辐照下,改变血液光学特性参数时xz平面(y=0)的光分布,相应的纵向光能流率分布如图11所示。与血液1的光分布相比,因血液2的μ_a改变为原参数的1/20,故而更多光子进入静脉及其周围组织,导致血液区域(即z∈[0.7,1.04])的纵向光能流率衰减梯度减小,光的穿透深度增大。血液3中的μ_s改变为原参数的1/20,其光分布与血液1差异较小,光子在静脉周围组织中的分布略多于血液1,两者在血液区域的纵向光能流率分布基本重合;在血液区之外(即z∈(1.04,1.4]),血液3的纵向光能流率略大于血液1。血液的散射系数减小,更多的光子发生前向散射而进入更深层的组织中。

图 10. 不同光学特性参数血液在xz平面(y=0)的光分布。(a)血液1;(b)血液2;(c)血液3

Fig. 10. Light distributions of xz plane (y=0) for different optical parameters of blood. (a) Blood 1; (b) blood 2; (c) blood 3

下载图片 查看所有图片

图 11. 不同光学特性参数血液的纵向光能流率分布

Fig. 11. Longitudinal distribution of fluence of different optical parameters of blood

下载图片 查看所有图片

4 结论

在构建静脉血管及其周围组织管状模型的基础上,将激光腔内治疗静脉曲张的360°径向辐射式光源简化为各向同性点光源,采用三维蒙特卡罗算法模拟其在组织中的光分布情况。比较了波长为810,940,980,1320,1500 nm的激光辐照静脉腔内时光能流率的分布差异:在相同剂量的情况下,1500 nm激光辐照1倍管径血管时,光子在血液及血管壁中几乎全部被吸收,对血管周围组织的损伤最小。比较了在同等剂量的1500 nm和1320 nm波长激光辐照下,血管粗细对光能流率的影响:1倍管径以下的较细血管适合选择1500 nm激光治疗,1500 nm激光对1倍管径以上的血管壁损伤较小,无法达到手术的预期目的;1倍管径以上的较粗血管适合选择1320 nm的激光,相比于较短的波长,1320 nm激光对周围血管组织的损伤较小。分析了光源位置对激光腔内治疗的影响:光源偏离血管中轴线越远,靠近光源一侧血管壁的光能流率越大,越容易引起该侧血管壁发生溃疡和穿孔;而另一侧光能流率很小,无法达到良好的治疗效果。比较了1320 nm激光辐照下,血液光学特性参数对光能流率的影响:减小血液吸收系数和散射系数均有利于更多的光子进入血管周围组织,并损伤周围组织。

不同于已有研究采用平面多层结构模型构建管状三维组织模型,本研究采用蒙特卡罗算法模拟静脉腔内治疗过程中的光分布,在综合考虑畸形血管管径以及光能流率分布的基础上,选择可以凝固血液、闭合血管且对周围组织损伤最小的波长。激光腔内治疗过程中光纤的回撤速度对光分布的影响,以及在静脉血管光热响应基础上确定安全的治疗剂量还有待进一步研究。