近年来，磁纳米热疗技术已成为一种新兴的肿瘤热消融技术，相比传统肿瘤治疗它具有安全性高、定位准确、副作用小等优点[1-2]．含有磁纳米粒子(magnetic nanoparticles，MNPs)的磁流体在磁热疗前被靶向注入肿瘤区域，随后MNPs在交变磁场的作用下产热并使肿瘤区温度提升至损伤癌细胞所需理想治疗温度(通常是42°C~46°C)[3]．精准控制肿瘤区域的治疗温度大小和分布是决定磁热疗治疗效果的关键，然而其在治疗过程中受诸多因素影响，如外加磁场属性、MNPs属性、生物组织属性和磁流体注射策略等．其中，治疗温度的大小主要受外磁场和MNPs属性影响，而肿瘤区域内的温度分布则与生物组织属性和磁流体注射策略关联较大[4]．因此，鉴于磁热疗治疗的复杂性，肿瘤区域治疗温度的控制应可根据实际状态动态调整至理想温度．尽管以往的报道已对磁热疗的方方面面做了较多的贡献，但同时仍存在一些亟需解决的科学问题．在以往的研究中，MNPs一般被假定为尺寸恒定的球体且均匀分布在治疗区域内[5]，或者为简化研究将MNP设定为一个或多个大尺寸的MNP核[6]．此外，Iacob等[7]对长短径比不同的椭圆形状MNPs产热功率进行比较，并得到了最佳传输功率下的长短径比；Seehra等[8]对不同尺寸MNPs的尺寸依赖特性进行了测量，并证明了MNPs磁性和弛豫参数的尺寸依赖性．事实上，技术的局限性导致制备过程中并非所有的MNPs都具有相同的尺寸，文献[9]指出铁磁流体的粒径分布与对数正态分布较为相符．值得一提的是，MNPs的多分散性也会明确影响其产热并最终影响治疗温度[10]，然而现有的报道并未对这个特点进行有效关注．此外，磁热疗温度控制系统方面也亟需进一步深入研究，因为癌细胞只有在特定的治疗温度范围才会因其高热敏感特性受到损伤甚至凋亡．Haemmerich等[11]创建了一个由有限元模型和控制施加电压的软件组成的热疗闭环系统，该系统通过优化控制器的控制参数，允许以最小的输入实现温度控制射频消融．以上研究均将热疗的控制系统近似等效为一个传递函数，然而实际磁热疗系统的非线性特点决定了这些等效并不合理．这里基于实际文献的多分散MNPs标准差构建了肿瘤区域内MNPs分布模型，并分析单分散性下均匀和随机MNPs分布对磁热疗温度分布的影响．此外，针对磁热疗系统的非线性特点，提出一种基于反馈的模糊自适应PID控制系统，用于实时调控磁热疗期间的温度使其始终保持在最佳治疗温度46℃，同时可以有效补偿MNPs分布不均匀性和多分散性造成的局部温差，从而获得最佳治疗效果．1 原理及方法1.1　数学模型1.1.1　生物传热模型磁热疗过程中，肿瘤区域的MNPs在外部交变磁场的作用下产生治疗所需的热量，将MNPs产热量作为Pennes生物传热(Pennes biological heat transfer，PBHT)方程的输入，应用有限元方法便可预测生物组织的温度分布情况．MNPs靶向注射区域为肿瘤组织，因此可分别用两个不同的PBHT方程描述肿瘤组织和健康组织的温度分布[12]，具体为ρTcT∂TT/∂t=QmT+ωbρbcb(Tb-TT)+kT∇2TT+PT; (1)ρHcH∂TH/∂t=QmH+ωbρbcb(Tb-TH)+kH∇2TH，(2)式中：下标T和H分别表示肿瘤和健康组织；参数ρ，c，k和T分别为密度、比热容、导热系数和绝对温度；参数ωb，ρb，cb和Tb分别为血液灌注率、血液密度、血液比热容和血液温度；Qm和PT分别为生物组织代谢热量和MNPs在交变磁场下的产热．应用有限元方法求解偏微分方程(1)和(2)时，设定肿瘤和健康组织的初始温度应为正常人体温度，即37 °C，生物组织治疗温度应存在上限[13]，TTt=0=THt=0=37 °C;TT(t)Tlim,TH(t)Tlim,Tb(t)Tlim,式中Tlim为导致健康组织坏死的治疗温度．此外，两个不同组织界面应满足Newman条件，具体为TT(t)|γ=TH(t)|γ;kT∇TT(t)|γ=kH∇TH(t)|γ,式中γ表示肿瘤和健康组织间的边界．1.1.2　磁纳米粒子产热模型磁纳米粒子在交变磁场下的产热与其粒径大小直接相关，假定MNPs粒径满足超顺磁特性，因此其产热机制主要由Brownian弛豫和Néel弛豫主导．这样，单个MNP在交变磁场作用下的功率耗散可用Rosensweig理论进行描述[10]，具体为PT=πμ0χ0fH02{(2πfτ)/[1+(2πfτ)2]}，式中：μ0为真空中的磁导率；χ0为平衡磁化率；f 为磁场的频率；H0为磁场的强度；τ为有效弛豫时间．1.1.3　磁纳米粒子的粒径分布MNPs的粒径分布可分为尺寸匀一的单分散性，以及尺寸按一定规律分布的多分散性．通过TEM显微镜可对实际制备的磁流体中某一区域的MNPs的尺寸和数量进行统计，得到MNPs的尺寸分布，而统计数据符合对数正态分布[14]，即采用对数正态分布函数对MNPs粒径分布进行建模，g(D)=12πσDexp-(ln{D}mm-μ)22σ2;∫0∞g(D)dD=1,式中：D为MNP直径；μ和σ具体为μ=ln[E(D)]-2-1lnVar(D)/E2(D)+1,σ=lnVar(D)/E2(D)+11/2,ED为D的期望，VarD为D的标准差．1.2　磁纳米粒子建模方法所提出的具有多分散性MNPs模型的构建具体流程如图1所示．首先基于Matlab软件环境确定MNPs对数正态分布的期望和标准差[9,15]，随机生成符合对数正态分布MNPs的直径；随后生成MNPs的位置坐标，同时应确保MNPs不发生重叠且不能超出肿瘤组织区域；重复以上步骤，直至肿瘤组织区域内的MNPs超过指定数量M=1 000；单个MNP的功率耗散为Qi，热疗总功率为所有MNPs的功率耗散之和，即Q=∑i=1MQi．10.13245/j.hust.240472.F001图1MNPs建模方法流程图1.3　控制策略1.3.1　比例积分微分控制器MNPs的多分散性增加了磁热疗治疗温度预测的不确定性，比例积分微分(PID)控制方法则可在一定程度上实现对目标区域温度的控制并最终改善治疗温度．PID控制器由比例环节、积分环节和微分3个不同环节构成，通过反馈回路对被测参数进行调控，以保证被测参数的稳定性[16]，即u(t)=Kpe(t)+Ki∫0te(t)dt+Kd∂e(t)/∂t，式中：e(t)=yd(t)-y(t)为误差；Kp为比例系数；Ki为积分时间常数；Kd为微分时间常数．1.3.2　模糊自适应PID控制器模糊自适应控制器是一种非线性控制器，控制过程中无须建立精确的传递函数即可实现实时控制．模糊逻辑控制器的输入为偏差e(t)和偏差的变化率ce(t)，输出为Kp，Ki和Kd[17]．模糊逻辑控制器首先通过模糊语言对偏差e(t)和偏差的变化率ce(t)进行模糊化得到模糊量，然后根据模糊规则进行模糊推理得到对应的隶属度值，最后将隶属度值解模糊化得到精确值[18]．模糊化过程采用七段式的模糊方法，将e(t)和ce(t)转化为NB，NM，NS，ZE，PS，PM和PB，表示的含义分别为负大、负中、负小、零、正小、正中和正大．模糊规则如表1所示．采用重心法求解模糊输出为[19]10.13245/j.hust.240472.T001表1模糊规则(∆Kp/∆Ki/∆Kd)e(t)ce(t)NBNMNSZEPSPMPBNBPB/NB/PSPB/NB/NSPM/NM/NBPM/NM/NBPS/NS/NBZE/ZE/NMZE/ZE/PSNMPB/NB/PSPB/NB/NSPM/NM/NBPS/NS/NMPS/NS/NMZE/ZE/NSNS/ZE/ZENSPM/NB/ZEPM/NM/NSPM/NS/NMPS/NS/NMZE/ZE/NSNS/PS/NSNS/PS/ZEZEPM/NM/ZEPM/NM/NSPS/NS/NSZE/ZE/NSNS/PS/NSNM/PM/NSNM/PM/ZEPSPS/NM/ZEPS/NS/ZEZE/ZE/ZENS/PS/ZENS/PS/ZENM/PM/ZENM/PB/ZEPMPS/ZE/PBZE/ZE/NSNS/PS/PSNM/PS/PSNM/PM/PSNM/PB/PSNB/PB/PBPBZE/ZE/PBZE/ZE/PMNM/PS/PMNM/PM/PMNM/PM/PSNB/PB/PSNB/PB/PBz0=∑i=0nμc(zi)⋅zi/∑i=0nμc(zi)，式中：zi为论域内的值；μc(zi)为zi隶属度值．根据表1所示的模糊规则可实现对∆Kp，∆Ki和∆Kd的实时调整．模糊逻辑控制器输出∆Kp，∆Ki和∆Kd按照如下方式组合作为PID控制器的输入[20]，即Kp=Kp0+ΔKp;Ki=Ki0+ΔKi;Kd=Kd0+ΔKd,式中：Kp0，Ki0，Kd0为PID控制参数的初始值；∆Kp，∆Ki和∆Kd均为模糊逻辑控制器输出值．2 物理特性与参数根据参考文献总结了用于模拟研究的常数材料特性．表2列出了用于求解PBHT方程的材料参数 [21]．磁场和MNPs相关参数分别：Hm=15 kA/m，f=300 kHz，τ0=1×10-9 s，δ=2.35×10-3 kg∙m-1∙s-1，Md=446 kA/m，Keff=23 kJ/m3，μ0=4π×10-7 N/A2，kb=1.38×10-23 J/K[22-23]．10.13245/j.hust.240472.T002表2用于生物传热方程的材料参数参数肝脏肿瘤血液Fe3O4ρ/(kg∙m-3)1 0791 0501 0505 180k/(W∙m-1∙K-1)0.520.590.5940c/(J∙kg-1∙K-1)3 5403 7503 750670Qm/(W∙m-3)14 4146 842——3 结果与讨论3.1　物理模型为了验证所提出的方法，基于正常人类二维肝脏模型构建所需的几何模型．图2中的几何模型包括健康肝脏组织、被健康组织包围的肿瘤组织，以及分布在肿瘤组织内的MNPs模型．肝脏模型左右径约为250 mm，上下径约为160 mm．模型的坐标原点设定为肿瘤组织中心位置，肿瘤组织被设计为一个半径15 mm的圆．为提高可读性，后续将默认采用局部放大视图对肿瘤组织温度预测结果进行观察，局部放大视图为长度为20 mm的正方形区域．10.13245/j.hust.240472.F002图2二维几何模型3.2　单分散性MNPs下的治疗温度分布磁热疗期间MNPs粒径大小一样时，其在肿瘤区域内的分布也可以分成理想均匀分布和随机分布．如图3所示，这两种不同MNPs分布在相同耗散功率下(1.51×106 W)在肿瘤区域呈现出治疗温度分布上的差异．此时，理想MNPs分布对应的情况其最高温度刚好位于安全临界值46°C，而随机MNPs分布的最大值比前者小0.4 °C．此外，随机MNPs分布情况下的42 °C等温线所覆盖面积相比前者范围也更大，故其对健康组织的损伤影响也相对较大．相反地，均匀MNPs分布情况对应的42 ℃等温线则比较均匀，紧紧包裹在肿瘤区域外．10.13245/j.hust.240472.F003图3不同单分散MNPs分布下的治疗温度分布(色标单位：°C)3.3　多分散性MNPs下治疗温度分布多分散性MNPs粒径分布情况与其对应的标准差高度相关，选定了4个不同的标准差用于分析其对MNPs粒径分布的影响，分别为0.1，0.194，0.209，和0.5，其中，0.194和0.209来源于参考文献的实际数据[12,19]，而0.1和0.5则是为了表达差异化而设定的．由图4可见：期望值相同的情况下，MNPs粒径分布标准差越大，MNPs粒径分布偏离中心粒径的范围越大，即MNPs粒径分布越分散，分布曲线更加矮胖．相反地，MNPs粒径分布标准差越小，MNPs粒径分布偏离中心粒径的范围越小，即MNPs粒径分布越集中，分布曲线更加瘦高．10.13245/j.hust.240472.F004图4不同标准差下多分散性MNPs的粒径分布接下来将就两种真实多分散MNPs分布的标准差(0.194和0.209)对治疗温度分布影响进行进一步讨论．图5描述了相同功耗作用下不同标准差多分散性MNPs分布对肿瘤区域治疗温度分布的影响．可见MNPs方差较小情况(σ=0.194)在相同治疗条件下其可获得的最大治疗温度为45.7 °C，较另一种情况(σ=0.202)高0.3 °C．此外，前者在肿瘤区域内的3条等温线覆盖的面积也较后者大．因此，相同条件下，多分散MNPs分布的标准差越小，粒径分布的范围越窄，治疗温度的最大值反而越大，同时相同等温线覆盖的有效治疗面积也越大．或者也可总结为MNPs大小差异越小，即越接近均匀分10.13245/j.hust.240472.F005图5相同功耗不同多分散MNPs分布下的治疗温度分布(色标单位：°C)布，其所对应的磁热疗治疗温度效果越优．3.4 磁热疗温度分布评估与控制系统优化基于上述分析，对单分散性均匀分布(MU)、单分散性随机分布(MR)和两种不同标准差下的多分散性随机分布(PRn)等4种不同MNPs分布对磁热疗治疗效果影响进行评估，这里变量n取值为1和2分别表示标准差为0.194和0.209两种情况．表3给出了通过两种不同指标用于评估4种不同情况下磁热疗治疗效果，即超出治疗区域比例(proportion beyond treatment area，PBA)和最佳治疗区域比例(proportion for optimal treatment area，POA)．其中，前者表示健康组织区有效治疗温度与肿瘤面积的百分比，用于描述治疗过程中对健康组织的损伤程度．后者表示肿瘤区内治疗温度在45°C ~ 46°C范围面积与肿瘤面积的百分比，此范围下的温度治疗效果最优．由表3可以看出：均匀MNPs分布的PBA最小仅为2.14%，而随机MNPs分布对应所有情况下的PBA均高于均匀分布情况．由此可知：MNPs分布均匀性的增加可有效减少健康组织的伤害．此外，具有单分散性且均匀分布的MNPs模型可获得最大的POA(40.37%)，即治疗效果最佳．在所有随机分布情况中，单分散优于多分散，而多分散中标准差越小POA越大．10.13245/j.hust.240472.T003表3MNPs不同分布下不同模型的最佳治疗面积所占百分比评价指标MUMRPR1PR2PBA2.147.834.715.19POA40.3726.5422.8114.99%为了验证所提出的方法，基于Matlab软件中的可视化仿真工具Simulink搭建了如图6所示的基于模糊自适应PID(fuzzy adaptive PID，FAPID)磁热疗控制系统框图．实际上，磁热疗加热机制相对复杂，生物组织特性容易受到外部环境的影响，而MNPs分布的不均匀性以及分散性也同样会影响治疗区域的温度分布情况．所提出的模糊自适应PID磁热疗控制系统通过动态控制克服了难以对生物传热过程建立精确传递函数的困难；同时，该系统可对具有多分散性且随机分布MNPs所导致的复杂肿瘤区治疗温度分布情况进行有效地实时调控．在不进行干预的情况下，磁热疗的肿瘤区域温度难以精确地控制在有效治疗范围，过高的温度会组织造成不可逆的损伤而过低的温度则无法对癌细胞造成损伤．因此，磁热疗控制系统的最高目标温度可设定为安全温度上限值46 °C，在满足此条件下使治疗区域温度尽可能高．10.13245/j.hust.240472.F006图6模糊自适应PID磁热疗控制系统磁热疗期间保证肿瘤区域最大治疗温度始终保持或者小于上限值46 °C，则其他位置的温度一定也是安全的．本控制系统接下来主要实现对该位置进行实施调控．鉴于此，将基于搭建的FAPID磁热疗控制系统对系统的性能进行了评估，其中系统初始功率设为经验值1×105 W．图7描述了相同的条件下，PID和FAPID对考虑MNPs多分散且随机分布情况下磁热疗系统最高温度随时间变化的控制曲线，图中目标温度(θ)为46 °C．由图7可知：两种控制器均可最终实现稳态治疗温度最大值控制在目标温度．然而，PID控制系统的稳定性较差，温度调控曲线出现震荡，稳定的时间长，温度超调的情况较为严重．FAPID控制系统则可实现极小的超调量和较快的稳定速度，因而其对健康组织的损伤也极小．10.13245/j.hust.240472.F007图7不同控制器下多分散MNPs系统的实时温度控制为了更准确说明所提FAPID磁热疗控制系统性能的优越性，通过调节时间和超调量对传统PID控制系统和FAPID控制系统的性能进行了评估．由表4可以看出：FAPID控制系统的调节时间为62 s，而PID控制系统调节时间为110 s．由此可见：相比于PID控制系统，采用FAPID控制系统能更快速调节温度，使得治疗区域温度维持相对稳定．除此之外，FAPID控制系统的超调量为0.89%也远小于传统PID控制系统超调量(8.87%)，可说明FAPID控制系统具有较好的稳定性，能较有效地防止治疗局部温度过高．10.13245/j.hust.240472.T004表4PID控制系统与FAPID控制系统性能对比控制系统类型调节时间/s超调量/%PID1108.87FAPID620.89相同条件下，MNPs分布的均匀性与磁热疗治疗温度分布及治疗效果成正相关关系，而MNPs粒径分布的多分散性程度也会同样影响治疗温度．MNPs粒径标准差越大，粒径分布范围越宽，粒径分布越分散，肿瘤区域内治疗温度最大值和有效面积均越小．因此，提高MNPs分布的均匀性和MNPs粒径分布的均一性都可以有效提升治疗效果，改善温度分布的均匀性．此外，所提出的模糊自适应PID磁热疗系统能对因MNPs分布的不均匀性和MNPs粒径分布的不均一性造成的局部温度异常进行相应的补偿，从而使得温度最大值稳定在最佳安全温度46 ℃．同时，与传统PID控制系统相比，模糊自适应PID磁热疗控制系统能迅速调节治疗温度并具有显著的稳定性，因此也适用于其他组织的肿瘤模型．磁热疗过程中治疗温度受诸多因素影响，导致其很难被精确预测，而所提出的控制系统可一定程度上弥补温度预测模型的不足．尽管如此，磁热疗期间仍存在值得进一步深入研究和探讨的问题，如组织和MNPs的温度依赖性问题、MNPs通过不同方式进入组织后的实际分布行为、组织多孔率对传热预测的影响、模型在微纳尺度下的有效性问题，以及磁热疗微观机制问题等．