纳秒电脉冲对肝脏组织不可逆电穿孔消融区分布的时域有限差分法仿真

2021-10-11胡亚楠包家立朱金俊朱朝阳岑泽南

电工技术学报 2021年18期

胡亚楠包家立朱金俊朱朝阳岑泽南

（浙江大学医学院浙江省生物电磁学重点实验室生物物理与医学工程研究组杭州 310058）

0 引言

纳秒电脉冲是指脉冲宽度在纳秒级的一类脉冲电场，这类电场作用在细胞膜上可以产生可逆电穿孔（Electroporation）[1]和不可逆电穿孔（Irreversible Electroporation, IRE）[2]。可逆电穿孔是脉冲电场作用于细胞，在细胞膜上产生一种可逆的瞬时性亲水孔洞的现象[3]，最早由E. Neumann和K. Rosenheck发现电脉冲引起囊膜的渗透率变化[4]，之后，U. Zimmermann等发现1～10kV/m脉冲电场可以击穿红细胞膜，形成细胞膜介质击穿（dielectric break down），血红蛋白释放，溶血机制引起细胞机械性破裂或热破裂，在37℃破裂细胞重新封闭，这种现象称为细胞膜电穿孔[5]。然而，A. J. Sale等发现，在高压脉冲电场下，红细胞和原生质体产生溶菌作用，细胞出现死亡[6]。孙才新等认为，随着脉冲电场剂量增高，细胞膜出现不可逆恢复的破裂导致细胞死亡，把这种现象称为不可逆电击穿（irreversible electrical breakdown），并观察到20～250V脉冲电压下，人卵巢腺癌SKOV3细胞出现外形留存、细胞膜完整，但线粒体肿胀、空化，表明高压脉冲电场可以致细胞死亡[7]。

2011年，美国食品药品监督管理局（Food and Drug Administration, FDA）批准了纳秒刀（nanoknife）为一种医疗器械，之后，有关纳秒电脉冲的IRE临床应用研究出现了爆发性增长[8]，尤其是恶性肿瘤的治疗[9]，如肝癌、胰腺癌、前列腺癌、肺癌、胃癌、乳腺癌、肾肿瘤、胆囊外周肿瘤、甲状腺癌、直肠癌肝转移、盆腔肿瘤，以及疏通血管、促进血脑屏障、治疗心律不齐等。一些新的IRE技术也正在研发，如高频IRE[10]、复合脉冲IRE[11]等。

不可逆电穿孔技术的关键是脉冲电场的非热效应[12]。尽管肿瘤细胞经过高压电脉冲处理后，细胞膜外表形态没有变化，但细胞核和细胞器已经受损、死亡[13]。然而，脉冲电场能量通过生物组织释放转化为热能是不可避免的，在多个电脉冲发放之后，组织的温度上升也是必然的[14]。电脉冲能量的提高增加了IRE的技术风险性[15]。实际上，电脉冲能量是可逆电穿孔和不可逆电穿孔的关键因素，也是热效应和非热效应的关键因素[16]。因此，在IRE手术前，预测术中和术后的安全有效性，指导手术按计划执行是十分必要的[17]。

组织内电场分布与电极构型、脉冲参数和组织异质性有关，使IRE治疗计划变得复杂。为了解决这一问题，施行IRE术前仿真模拟，获得最佳脉冲参数，优化治疗计划是非常重要的[18]。IRE术前仿真模拟的主要内容是预测器官和组织在各种脉压、脉宽、脉率的电脉冲作用下的电场分布、能量分布和温度分布[19]，评价该电脉冲参数的手术适用性和安全性，指导手术方案的制定。有限元分析是IRE术前仿真模拟的有力工具[20]。其中，时域有限差分法（Finite Difference Time Domain, FDTD）是电磁生物仿真的常用方法[21-22]。

FDTD是K. S. Yee于1966年建立起来的[23]，其核心思想是把带时间变量的Maxwell旋度方程转化为差分形式，模拟出电脉冲与理想导体的时域响应。本文运用FDTD方法对300ns电脉冲作用在肝脏组织的电场、能量密度和消融区域分布进行仿真，研究脉冲波形、脉冲电压、脉冲宽度、电极间距对电场、能量密度和消融区域分布的影响。

1 FDTD原理

1.1 介质的麦克斯韦方程

在介质中，麦克斯韦方程的微分形式为

式中，E为介质中的电场强度（V/m）；D为介质中的电位移（C/m2）；B为介质中的磁感应强度（T）；H为介质中的磁场强度（A/m）；je为介质的电流密度（A/m2）；ρ为介质的电荷密度（C/m3）。介质的仿真条件为：

（1）介质具有各向同性性质，其电导率σ、电容率ε和磁导率μ为常数，并且有D=ε E，B=μ H。

（2）介质内无电源，则有ρ=0。

（3）介质有电磁损耗，电损耗的电流密度为je=σE，磁损耗的磁电流密度为jm=sH，s为磁阻率。

这样，介质内的电场和磁场旋度微分形式为

1.2 FDTD蛙跳算法

在介质形状不对称、边界特性不一致、初始条件不确定的情况下，用式（2）计算介质内的电场Et(x,y,z)和磁场Ht(x,y,z)分布是困难的。FDTD 是解决上述问题的一种方法，其核心思路是把介质实体分割为有限个体元（Δx, Δy, Δz），体元上电场和磁场分量如图1所示，时间分割为有限个时间步Δt。所有电场分量在体元棱中点，磁场分量在面心，电场和磁场的分量正交。

图1 体元上电场和磁场分量 Fig.1 Electric and magnetic field components on a voxel

在t=nΔt时，空间点（inΔx,jnΔy,knΔz）的体元（in,jn,kn）的电磁场的可离散表达为

在该方法中，空间和时间都可以被等距离散化，每个离散间隔内的电场和磁场可认为是线性的。以时间计算为例，时域电场En和磁场Hn可以用两个离散点的平均值来表示，有

于是式（2）离散化为

可得

用蛙跳算法（leap frog algorithm），即式（6）交替计算时域的电场En和磁场Hn，如图2所示。进而，可以建立空间内各个体元的蛙跳算法，模仿电磁场在空间域的变化，直至计算出仿真全程的En(i n,jn,kn)和Hn(in,j n,kn)。

图2 蛙跳算法 Fig.2 Leap frog algorithm

1.3 能量密度分布

肝组织是一种密度均匀的电介质，暴露在电磁场中所聚集的能量密度为

根据蛙跳算法在空间各点的电场E(x,y,z)和磁场H(x,y,z)，再用式（7）可以计算空间（x,y,z）肝组织的能量密度分布。E(x,y,z)和H(x,y,z)分别为特定频率下电场强度和磁场强度分量。

1.4 消融区域分布

当脉冲电场在肝组织超过阈值，肝细胞出现IRE形成凋亡，这是纳秒电脉冲消融肝细胞的基本原理。肝细胞消融区域可根据电场的绝对值在肝组织中的分布进行预测。在电场达到最大时刻，电场阈值等值面包络空间内部即为消融区域。

2 肝脏组织电磁仿真

2.1 肝组织参数

采用人体生物电磁仿真软件（Sim4Life，IT’IS基金会，瑞士）对肝脏器官进行电磁仿真，仿真的体元大小为2mm×2mm×2mm。选择Sim4Life中人体模型Duke的肝组织，其参数见表1。肝组织为非铁磁性介质，可设s=0，μ=μ0=4π×10−7N/A2。

表1 Duke肝组织参数 Tab.1 The parameters of the liver tissue in Duke

2.2 边界条件

电磁波在体元边界上产生反射和透射，体元的边界条件有四种类型：

（1）单轴全匹配边界（Uniaxial Perfectly Matched Layers, UPML）：在肝脏模型以外适当的位置建立只有吸收、没有反射的长方体边界，将边界上电场与磁场数值强制赋值为0。如果电场磁场集中在很局限的范围内，弥散性不明显，则可以采用这种边界条件。该条件下的运算方程组同质性很强，处理速度较快，应用最为广泛。

（2）解析吸收边界条件（Analytical-Absorbing Boundary Condition, A-ABC）[24]：定义了吸收系数的边界条件。这种条件的模拟更加符合实际情况。该条件要求电磁场有较高频率，波长需短于边界尺度，并且运算量受入射角度等多种因素影响。异质性和损耗性的介质不适宜采用该边界条件。

（3）完全电导体边界条件（Perfectly Conductive boundary Condition, PEC）：将边界上的电场强度切向分量赋值为0。金属外壳屏蔽条件下的模型可采用该边界条件。

（4）完全磁导体边界条件（Perfectly Magnetic boundary Condition, PMC）：将边界（图1体元边界外半边长位置）上的磁场强度切向分量赋值为0。高导磁性材料屏蔽下的模型可采用该边界条件。

2.3 仿真的电脉冲

电极附近的电流密度最高，这是组织受热最强的位置。建模所选择的电脉冲参数应当能够反映放电速率、病变大小、组织温度和宿主反应等各种因素的实际总和。在选定参数时，考虑肌肉收缩为最小，不需要与心跳同步的脉冲。仿真的电脉冲主要指标是脉冲参数和电极构型，脉冲参数包括脉冲波形、脉冲电压、脉冲宽度、脉冲率、脉冲数等。脉冲波形有梯形波、指数波、方波三种，总脉宽为300ns，仿真的脉冲波形及其时间参数见表2。脉冲电压分别为30kV、10kV、3kV。

表2 仿真的脉冲波形及其时间参数 Tab.2 The parameters of the waveform and time for simulation

缺省电脉冲为梯形波TW50-50、脉冲电压10kV、电极间距10mm。

2.4 观测指标

仿真的脉冲电极插入肝脏位置的三视图如图3所示。电极构型为长度50mm，半径1mm，间距分别为5mm、10mm、20mm的平行双针，从肝的正前部位插入至肝中部（见图3c），电极中线距左叶尖端102mm，距右叶顶端51mm，插入部分长46mm。

（1）电场分布。两根电极导入到肝组织产生脉冲电场，电极与肝组织形成电导体边界条件，并利用FDTD方法计算出Et(x,y,z)。观察各种电脉冲参数，以及电极间距对电场分布影响，有利于对电脉冲参数和电极间距进行优化选择。

（2）能量密度分布。电磁能量是不可逆电穿孔肿瘤消融的关键物理因素，观察能量密度在肝组织内分布，有利于评价个体化的纳秒电脉冲治疗效果。

（3）消融区域。利用纳秒电脉冲治疗肿瘤的技术中，术前预测术后疗效是准确选择电脉冲参数的关键，消融区域是术后疗效的关键观察指标。从消融区域的仿真结果中，可以分析和判断各种纳秒电脉冲参数下肿瘤的消融效果，优选适宜于个体化的纳秒电脉冲参数。消融区域是根据FDTD仿真的电场强度进行阈值判断，低于阈值时，为非消融区域，反之，为消融区域。

图3 电极插入肝脏位置的三视图 Fig.3 The three-view drawing of electrode-inserted liver

2.5 电场阈值

消融区域主要由组织中的电场分布决定，这构成治疗计划算法的基础。电场阈值是确定消融区域仿真的依据，当体元内的电场强度超出阈值，模拟了组织细胞坏死，达到组织消融的目的。从细胞悬浮液的生存能力中很难确定细胞死亡的电场阈值。细胞死亡是以生存能力判断，悬浮细胞的电致生存能力估计小于5%[25]，实验的悬浮细胞致死电场阈值约为1 500V/cm，但组织体内电场强度明显低于悬浮细胞，约为500V/cm[26]。一般认为电场阈值在500V/cm[20]。本仿真采用的电场阈值为500V/cm。

3 仿真结果

肝脏器官中纳秒电脉冲电场分布、能量密度分布、消融区域的仿真结果是三维的，为了展示分布特征，采取在三维分布中最具特征的一个视角平面图。每张图左上角有强度标尺色图，对数分度，浅色最强，深色最弱。电场强度的范围是31.6V/m(−50dB)～ 1×104V/m(0dB)。能量密度取1MHz中的分量，范围是5×10−17J/m3(−50dB)～5×10−12J/m3(0dB)。消融区域以电场阈值为界，高于电场阈值为消融区域，标尺与电场强度相同。表3～表5中的例图是所示电脉冲期中电场强度或能量密度最强时刻（方波和梯形波的平台期，指数波的峰值），以及三维分布中电场强度或能量密度最强（电极平行双针所在平面），并从人体右前方，俯视45°方向（θ=−45°,ϕ=45°）的一个视角图。

3.1 脉冲波型

在脉冲电压10kV、电极间距10mm时，各种脉冲波型产生的仿真电场分布、能量密度分布和消融区域见表3。仿真结果显示，在肝组织及旁组织中，电极中心电场强度较高（中心点1.9×105V/m，内侧边缘2.7×105V/m），向四周扩散，至左叶尖端（302V/m），右下叶尖端（136V/m）。在肝内部，电场分布呈椭圆状，层次分明、连续、光滑，随电极距离变大，电场强度逐渐降低。在肝边缘凹进的部分，电场强度有所增加。在肝外部，由于胃、胰、脾、肾、脊液、皮肤、皮下脂肪等组织器官的电学参数各异、电场分布呈复杂的扭曲状。及至体外，电场分布重新呈现规则椭圆，并且较临近的皮肤内电场强度高，这是体外空气不导电且电容率低的缘故。在脉冲电压10kV、电极间距10mm时，肝组织对各种波型电脉冲的电场分布均没有显著改变。肝组织中能量密度分布基本集中在电极周围，且与脉冲波形有关，指数脉冲EW20电极能量密度中心点2.2×10−15J/m3，内侧边缘4.3×10−15J/m3，其他脉冲均中心点2.9×10−13J/m3，内侧边缘5.5×10−13J/m3。非电极周围肝组织能量密度低于电极区4～6个数量级（左叶尖端6.3×10−19J/m3，右下叶尖端3.6× 10−19J/m3）。消融区域是一个被电极包裹，截面长22mm、宽16mm的柱形区域，在电极尖端略有膨大。表3结果显示，电脉冲波型对消融区域影响不大。

3.2 脉冲电压

在50-50梯形波、电极间距10mm时，脉冲电压分别为3kV、10kV、30kV，产生的肝组织仿真电场分布、能量密度分布和消融区域见表4。仿真结果显示，3kV电脉冲的肝右叶下段月牙状电场分布明显缩小（尖端41V/m），肝组织中的电场强度（中心点5.8×104V/m，内侧边缘8.1×104V/m）明显小于10kV，肝左叶上段部位的电场分布不光滑、层次不分明，局部区域明显减小（尖端91V/m）。30kV的肝组织电场分布（中心点5.8×105V/m，内侧边缘8.1×105V/m）明显大于10kV，肝右叶下段及其旁组织的月牙状电场分布（尖端409V/m）明显大于10kV。以缺省电脉冲肝组织能量密度分布区域中心2.9×10−13J/m3，内侧边缘5.5×10−13J/m3为对照，3kV电脉冲能量密度分布区域（中心点2.6×10−14J/m3，内侧边缘4.9×10−14J/m3）明显小于10kV电脉冲。30kV中心2.6×10−12J/m3，内侧边缘4.9×10−12J/m3，明显高于10kV，但仍局限在电极周围。以缺省电脉冲消融区域22mm×16mm柱形域为对照，3kV电脉冲的消融区域为12mm×6mm，小于10kV，并上下向内凹陷呈哑铃状。30kV电脉冲的消融区域为36mm×30mm，明显广于10kV。表4结果显示，脉冲电压愈高，电场强度、能量密度、消融区域均愈高。

表3 脉冲波型对电场分布、能量密度、消融区域的影响 Tab.3 The influence of pulse waveform on electric field, energy density and ablation zone

3.3 电极间距

在50-50梯形波、脉冲电压10kV时，电极间距分别为5mm、10mm、20mm，产生的肝组织仿真电场分布、能量密度分布和消融区域见表5。仿真结果显示，5mm电极间距的肝右叶下段月牙状电场分布明显缩小（66V/m），肝组织中的电场强度（中心点3.8×105V/m，内侧边缘4.4×105V/m）明显强于10mm，肝左叶上段部位的电场分布不光滑、层次不分明区域明显减小（145V/m）。20mm的肝组织电场分布（中心点7.2×104V/m，内侧边缘2.3×105V/m）明显弱于10mm，肝右叶下段及其旁组织的月牙状电场分布明显较大（227V/m）。以缺省电脉冲肝组织能量密度分布区域中心2.9×10−13J/m3，内侧边缘5.5× 10−13J/m3为对照，5mm电极间距的肝组织能量密度分布区域中心1.1×10−12J/m3，内侧边缘1.5×10−12J/m3，明显强于10mm电脉冲。20mm中心4.4×10−14J/m3，内侧边缘4.0×10−13J/m3，明显弱于10kV，但仍局限在电极周围。以缺省电脉冲消融区域22mm×16mm为对照，5mm电极间距消融区域14mm×14mm，明显小于10mm。20mm电极间距消融区域33mm× 17mm，明显广于10kV，同时也上下向内少许凹陷呈哑铃状，仍局限在电极周围。表明电极间距越大，电场强度、能量密度越低，消融区域越大。

表4 脉冲电压对电场分布、能量密度、消融区域的影响 Tab.4 The influence of pulse voltage on electric field, energy density and ablation zone

表5 电极间距对电场分布、能量密度、消融区域的影响 Tab.5 The influence of electrode spacing on electric field, energy density and ablation zone

4 讨论

4.1 仿真模型

除了FDTD模型之外，还有一种准静态方程，或称Laplace方程，其特点是无时域变量，电场仅是空间函数[18-20]，有

式中，ϕ为电动势；σ为组织电导率。电场是电势的梯度，E=−∇ϕ。FDTD所需的算力要求远高于准静态，但准静态一般采用预设标准电压，这样，场源波形无法在仿真中得以体现，无法用不同激励源的脉冲电场得出一般规律。而FDTD可以通过记录全程的波形变化，得出一个同一标准的仿真结果。

在300ns电脉冲，用准静态模型仿真可以得到，发出脉冲到检测到脉冲有1～2ns的延迟，表明这种纳秒电脉冲产生的电磁场波动性不强。在波动性不强的电磁场中，采用准静态模型仿真是可行的。但脉冲宽度变窄，频谱变宽，波动性增强，则准静态模型有较大偏差，需要用FDTD模型仿真。在消融效应上，可以根据FDTD仿真结果对不同电场的肝脏组织坏死/损伤区域进行分析，获得经验规律，进而定量研究脉冲电场对细胞/组织的效应。

网络传输网络模型是一种模拟传输的有限体积方法，可以仿真在纳秒脉冲下单细胞电穿孔的动态变化过程，仿真出孔道半径有一定的累积效应[27]。

4.2 软件平台

Laplace方程仿真软件通常是Comsol多物理场软件。这个软件是一种借助数值仿真理解、预测和优化工程设计的通用软件平台，涵盖力学、电磁场、流体、传热、化工、微电子、声学等，可用于任何基于物理系统的建模和模拟，也可以将模型封装为仿真应用软件，提供给设计、制造、实验测试。大部分不可逆电穿孔肝癌消融仿真都是在这个软件平台上进行[26,28-30]，其优点是可以在网上进行仿真工作。

Sim4Life是专业的人体生物电磁仿真平台，该软件提供的人体和动物模型来自于实体检测的数据。尽管个体差异使人体和动物模型的数据会有一点差异，但是，这种仿真应该更接近于实际人体的仿真。这些特定的人体和动物对研究个体化生物效应具有意义。

4.3 能量密度与SAR

生物组织中的电磁能量将转化为热量。电磁波量转化为辐射热，符合坡印廷定律，电磁场能量转化为焦耳热，符合欧姆定律。生物体本身的生化反应产生热量，称为自产热，血流的对热换热构成了自然生态环境。电磁热不论是辐射热，还是焦耳热，都是外源性热量，改变了生物体的自然生态环境。因此，电磁热或者说电磁能量是影响生物效应的因素。能量密度可以观察电磁生物效应。生物体温度与热量用Pennes生物热方程描述[31]为

式中，K为体元热传导系数；V为通过组织的血流容积；s为血液热量；Qm为组织代谢产热速率；Ta为动脉温；k为平衡常数0≤k≤1。C. K. Sung等的研究结果表明，当能量密度超过5.9×105J/m3，电极周边出现碳化迹象[32]。肝组织利用表1的肝组织参数和式（7）可以仿真出肝组织的能量密度均在10−13J/m3数量级左右，表明单个纳秒电脉冲在肝脏组织中非常低，可以认为不可逆电穿孔是非热效应。采用能量密度分布分析肝组织温度分布是合理的。

在脉冲电压足够低，组织电导率σ是常数，脉冲电压在可逆电穿孔阈值之上，式（9）是线性方程，仿真比较简单。如果脉冲电压超过组织局部电场阈值，组织电导率会增加，并取决于电场强度E，为σ(E)，式（9）为非线性，仿真将变得复杂[33]。另外，组织电导率σ也是温度T的函数：σ(T)=σ0[1+a(T−T0)][20,26]。其中，σ0为初始温度T0的初始电导率；a为系数。因此，组织电导率σ也是电场强度、温度的函数。

电磁生物效应与电磁场的波长有关。当电磁波长与人体或器官组织的尺度相比拟或更小，则电磁波效应占主导，器官组织获得的是电磁波能量，可用比吸收率（Specific Absorption Rate, SAR）计算。如30MHz电磁场的波长是10m，与人体尺度相近。当电磁波长比人体或器官组织的尺度大很多，则电场和磁场效应占主导，器官组织获得的是电场能和磁场能，要用能量密度计算。如1MHz电磁场的波长约300m，比人体尺度大2个数量级。因此，Sim4Life用式（7）计算肝组织吸收的电磁能量密度是合理的。

根据IEC/IEEE 62704标准，利用FDTD计算肝组织的SAR，需满足的频率范围为30MHz～6GHz电磁波[21-22]。然而，用Fourier频谱分析表2各种波形电脉冲的频谱是0～1MHz，远不足30MHz。因此，不能利用FDTD方法计算表2纳秒脉冲的肝组织SAR。实践上，在Sim4Life中计算也提示了结果警示。

4.4 肝组织电容率和电导率

能量密度式（7）显示，生物组织中的能量密度不仅与电场强度E和磁感应强度B有关，还与组织的电容率ε和磁导率μ有关，因为D=εE，B=μH。电容率是综合反映电介质畸变极化、位移极化、转向极化的宏观物理量，是外电场频率的函数，复数形式[34]为

式中，ε′(ω)为介质储能特性；ε′(ω)为介质损耗特性。纳秒电脉冲的频谱为0～1MHz，肝组织的电容率会受到电脉冲频谱的影响。因此，在仿真中应当采用肝组织的电容率色散特性。肝组织电容率色散谱具有个体差异，仿真中应采用肝组织电容率色散谱的实测值最为合理。

4.5 术前预测可行性

消融区域的预测可以帮助医师术前设计电脉冲计划，确定脉冲电压、脉冲时间、脉冲数和脉冲率，通过仿真预测术后效果。消融区域的确定方法有电场阈值[26]和温度[20]两种。有研究采用选择电脉冲参数预测术后效应[19]。在本仿真中，把脉冲电压倍增，电极间距减半，消融区域出现明显变化，而改变脉冲波形变化不大。因此，当选定了一种波形的脉冲后，主要是选择脉冲电压和电极间距来仿真术后治疗效果。