乳腺是由皮肤、纤维组织、乳腺腺体和脂肪组成，乳腺癌是一种发生在乳腺腺上皮组织的恶性肿瘤．乳腺癌的病因尚未完全清楚，但有研究发现乳腺癌的发病存在一定的规律性，有乳腺癌高危因素的女性易患乳腺癌[1]．由于不同类型的乳腺癌往往对现有的癌症治疗方法表现出不同的反应，因此治疗乳腺癌的一般治疗手段很难有普遍性成效．消除癌变的理想策略是将光与热相结合，实现非侵入性治疗，降低死亡率，提高患者的存活率[2]．目前，文献[3]的研究表明复合肽两亲性纳米胶束与光敏剂可以通过光的作用共同治疗，证明了其在体内治疗乳腺肿瘤的有效性．这种非侵入性的策略成为乳腺癌联合治疗一种极具吸引力的方法．文献[4]研究了光动力治疗(PDT)技术在影像引导下治疗原发性乳腺癌的有效性，并成功证实了PDT在早期乳腺癌治疗中的潜在作用．因此研究激光在乳房组织和乳腺肿瘤中的传输及之间的相互作用对乳腺癌的治疗具有重要的理论和现实意义．基于光粒子的光子传输理论应用最广，辐射传递方程(RTE)准确地描述了光在生物组织中的传播[5]，蒙特卡罗(MC)方法为光在弱吸收和高散射介质传播提供了准确的解决方案[6]．扩散近似(DA)是将RTE作为正向光子传播模型的一种典型方法[7]．DA由于其计算速度快，精度高而被广泛使用．目前有限元软件很多，大型通用有限元分析软件ANSYS能够很好地对非线性问题进行分析，主要应用于机械和核工业领域[8]，但使用十分复杂；有限元分析软件ALGOR操作简单，主要应用于工业产品设计中[9]；COSMOS软件影响力较小[10]．COMSOL Multiphysics是一款以有限元法为基础，通过数学方法求解偏微分方程组来仿真真实物理现象的软件[11]，主要优点是操作界面简洁，能够通过任意模块的组合实现优异高效的多物理场耦合，计算性能突出，仿真结果精确，对各种物理现象都能进行较好的计算．本研究选用COMSOL Multiphysics对乳腺癌光热成像进行可视化模拟．利用COMSOL构建乳房的简化模型并进行网格划分，随后对COMSOL中的系数形式偏微分方程(在数学模块中)和生物传热模块进行耦合．将系数形式偏微分方程简化近似为亥姆霍兹方程，实现对光在乳房以及肿瘤中的传播进行可视化模拟；用生物传热模块对温度变化进行可视化模拟．1 模型与方法1.1　几何模型利用COMSOL构建了一个二维半圆形模型来模拟乳房的局部截面．该模型对乳房实际结构进行了简化，如图1(a)所示．乳房半径为40 mm，乳腺肿瘤半径为4 mm，肿瘤的深度为37 mm，其横向位置偏离中线3 mm，乳房上方填充水层，水层总深度为45 mm，脉冲激光点源放置在水层中间以辐照乳房组织[12-13]．10.13245/j.hust.220310.F001图1模型的建立采用自由三角形网格对模型进行划分，设置最大单元尺寸为870 μm，最小单元尺寸为20 nm．最大单元增长率为1.4，曲率因子设定为0.2，狭窄区域分辨率为20．划分网格后模型如图1(b)所示．1.2　光沉积光在乳房组织中的传输使用COMSOL中系数形式偏微分方程进行模拟．系数形式偏微分方程通过亥姆霍兹方程[14]得到光通量率φ，具体为nc∂φ∂t+∇⋅(-D∇φ)+μa=f；(1)f=14π1.5Wpτpexp-4(t-τc)2τp2δ(r-r0)，(2)式中：c为真空中的光速；D=1/[3(μa+ μs')]为扩散系数；r表示空间坐标；t为时间；μs为散射系数；g为各项异性因子；μs'=μs(1-g)为约化散射系数；μa为吸收系数；f表示具有高斯分布的位于点r0处的激光源；激光脉冲能量Wp =4 mJ/cm2；激光脉宽τp=20 ns，在τc=30 ns处达到峰值．r0位于水层下部中央位置，模型的光学参数列于表1．10.13245/j.hust.220310.T001表1乳房模型光学参数对象μa/mm-1μs/mm-1肿瘤乳房水0.7410.0500.0060.900.890.101.3　热传递乳房组织吸收了光能量会产生温度T的变化，这种变化可利用COMSOL中的生物传热模块通过求解下面方程来模拟，即ρC∂T/∂t-∇⋅(k∇T)=μaφ，(3)式中：组织密度ρ=998 kg/m3；比热容C=840 J/kg∙m3；导热系数k=0.7 W/(m∙K)，组织的初始温度设置为310.15 K．为简化起见，本模型仅仅模拟乳房组织与乳腺肿瘤的温度变化．由于激光作用时间非常短，不考虑体温调节机制．2 结果与讨论2.1　光能量沉积仿真模型的剩余光能量分布如图2所示．图2(a)~(d)分别展示了在四个不同时刻乳房组织中的归一化光能量分布．颜色越蓝表示剩余的光能量越少，光能量吸收的越多．其中图2(a)为当t=26 ns时剩余光能量分布，由于此时激发光脉冲能量较小，光能量沉积较少；图2(b)为当t=28 ns时，激光源已经发射激光脉冲，但是激光脉冲的强度还未达到峰值，光能量沉积略有提高；图2(c)为当t=30 ns时，激光源发射的激光脉冲的强度达到峰值，此时整个乳房组织中光能量沉积继续增大；图2(d)为t=32 ns时，激光源发射的激光脉冲强度减弱，乳房组织中的剩余光能量在减小．10.13245/j.hust.220310.F002图2不同时刻的归一化光能量分布对光能量沉积进行定量分析，分别在横向和纵向进行分析．图3给出了z=37 mm处在t=26，28，30，32 ns这4个时刻的光能量分布，从图3可以看出：光能量在乳房组织与肿瘤之间呈近似指数下降，且在肿瘤中心处达到最小值．光能量在乳房组织和肿瘤处有明显的突变，这是因为这二者的光学参数不同，乳房肿瘤对光能量的吸收远大于乳房组织对光能量的吸收．图4给出了纵向穿过肿瘤即x=37 mm位置处的光能量分布，当t=26，28，30，32 ns时，肿瘤剩余的光能量低于周围乳房组织剩余的光能量，这是因为乳腺肿瘤的吸收系数显著高于乳房组织的吸收系数，肿瘤吸收了更多的光能量；在激光脉冲发射时刻(t=30 ns)，肿瘤内部剩余的光能量约为周围乳房组织剩余光能量的1.524%；同时肿瘤中心剩余的光能量略低于肿瘤边缘剩余的光能量．10.13245/j.hust.220310.F003图3z=37 mm处t=26，28，30，32 ns时刻的光能量10.13245/j.hust.220310.F004图4x=37 mm处t=26，28，30，32 ns时刻的光能量2.2　温度分布图5为获取的温度分布．可以看出：肿瘤温度高于周围乳房组织温度，肿瘤内部温度分布不均匀，而周围乳房组织的温度几乎均匀分布．进行定量分析，分别绘制横向和纵向的温度分布图．图6为横向穿过乳房组织和乳腺肿瘤中心处(z=37 mm)的温度分布图，图7为纵向穿过乳房组织和乳腺肿瘤中心处(x=37 mm)的温度分布图．从图6可以看出：肿瘤的温度显著高于肿瘤周围乳房组织的温度，温度差约为0.06 K，这是因为肿瘤的吸收系数大于乳房组织的吸收系数，肿瘤吸收了大量光子导致肿瘤升温显著．从图7可以看出：乳腺肿瘤中心的温度低于乳腺肿瘤边缘的温度，同时肿瘤靠近光源一侧的温度高于远离光源一侧的温度，这与光能量分布是符合的．10.13245/j.hust.220310.F005图5乳房和肿瘤温度分布的局部放大图10.13245/j.hust.220310.F006图6z=37 mm处横向穿过乳房和肿瘤的温度分布10.13245/j.hust.220310.F007图7x=37 mm处纵向穿过乳房和肿瘤的温度分布3 结论提出了一种基于COMSOL的脉冲激光与乳腺癌相互作用的光热成像可视化模型．首先，在COMSOL中构建了一个二维半圆形简化模型来模拟乳房的局部截面．采用COMSOL中系数形式偏微分方程和生物传热模块来实现乳腺癌的光热成像可视化．将系数形式偏微分方程简化近似为亥姆霍兹方程，实现对光在乳房以及肿瘤中的传播进行了可视化模拟；生物传热模块对光能量吸收导致的温度变化进行了可视化．结果显示肿瘤处剩余的光能量低于周围乳房组织剩余的光能量，其剩余的光能量约为周围乳房组织剩余光能量的1.524%．由于肿瘤吸收了更多光能量导致肿瘤的温度显著高于周围乳房组织的温度，温度差约为0.06 K．肿瘤中心的温度低于肿瘤边缘的温度，同时肿瘤靠近光源一侧的温度高于远离光源一侧的温度，这与光能量分布是契合的．可以看出：在入射光脉冲能量为4 mJ/cm2下乳腺肿瘤的温度上升对于乳腺癌的光热治疗还显得不足．更大的入射光能量可以获得更高的肿瘤温度上升，但容易导致光损伤．采用光热转换效率高的材料则是一种非常好的途径来增强光热治疗效果．未来的研究中，将进一步研究在有较高光热转换效率的材料下乳腺癌光热成像的可视化．该可视化模型只须要获得对应的结构特征、光学参数和热学参数，便可应用于其他类型肿瘤．本研究对乳腺癌的光热治疗具有一定的理论指导意义．