0 引言
近年来，脑调控技术已经在神经科学研究、脑功能疾病诊断干预等领域被广泛应用，不仅成为神经精神类疾病康复治疗的有效手段，同时也是开展神经系统工作机制研究、探索和认识大脑功能的重要工具[1-2]。以深部脑刺激、经颅直流电刺激、经颅磁刺激、经颅超声刺激等为代表的神经刺激方法是目前神经科学研究领域几种重要的神经调控方法[3]。其中，经颅直流电刺激和经颅磁刺激具有安全无创、操作简便的优点，被广泛用于脑神经活动的调控，但空间分辨率较低、刺激深度较浅等不足仍有待进一步完善[4]。经颅超声刺激是一种高空间分辨率的刺激手段，可以在大脑深部区域产生有效刺激效果[5-6]。
本文研究的经颅磁声电刺激（Transcranial Magneto-Acousto-Electrical Stimulation, TMAES）是一种新型无创式脑神经刺激和调控方法，该方法使用超声波和静态磁场共同作用于神经组织，可在神经组织周围产生感应电流来调节相应脑区或核团的神经活动。超声在生物组织中传播，聚焦精度可以达到毫米级别，TMAES感应电场分布与超声分布具有高度一致性，因而具有较高的空间分辨率。同时，相比于高频交变磁场，静磁场在生物组织中衰减大大减弱，使得刺激可以达到大脑深部区域。
通过电生理实验能够客观地观察外加物理场对神经系统活动的影响，但是受到实验条件的限制，很多实验无法应用于人体，实验结果难以复现。脑皮层神经网络包含了形态复杂的神经元胞体、树突、轴突等生物物理学细节，神经元的形态、相互之间的连接方式以及具体的电生理特性等信息都很难从实验的层次一一探究。建立皮层神经网络计算模型可以将不同层次的实验数据进行整合[7-8]，能够定量地刻画神经系统放电活动与自身生物物理特性之间的联系。因此，神经计算模型为研究电磁刺激的神经调节机制提供了一种重要方法，是目前揭示电磁场[9-11]生物效应的重要工具。
2003年，S. J. Norton提出了一种在磁场下施加超声产生感应电流，进而刺激神经组织的新方 法[12]。2015年，李慧雨等开展了基于电介质样本的聚焦超声电刺激实验，发现感应电场与超声分布存在高度一致性[13]。2015～2017年，Yuan Yi等进行磁声电刺激的理论分析，从仿真角度探索了磁声电刺激的神经作用机制[14]。2018年，刘志朋等进行小动物TMAES刺激实验，对TMAES的作用机制进行更深入的研究，进一步证明了TMAES的有效性[15]。2019年，本课题组基于H-H神经元模型，探索了TMAES下神经元放电模式的变化，为TMAES刺激参数的优选提供了一定的理论指导[16-17]。
以上研究并未分析TMAES中刺激参数变化对具有真实生物物理特性的神经网络放电活动的影响，而神经信号传递依赖于由神经元组成的神经元集群，因此，对神经元集群传递生物信息机制以及磁声电刺激的神经作用机理的分析至关重要。本文建立了可以反映脑皮层三维结构和生物物理特性的真实神经元模型，并搭建了神经网络，模拟了外加阶跃电流刺激、自身突触激活、TMAES三种方式下神经元局部场电位（Local Field Potentials, LFPs）响应特性。采用时频分析方法对不同刺激参数下阵列电极采集的神经网络放电数据进行处理，分析了感应电流密度与调制频率对神经网络LFPs信号能量分布的影响。
1 方法与模型
1.1 TMAES原理
生物体组织具有导电性，在静磁场条件下，沿着垂直于静态磁场的方向施加超声波，生物组织中的带电离子就会在超声波作用下沿着垂直于磁场的方向做往复运动。在洛伦兹力作用下，正负离子朝着相反方向运动，从而产生感应电流刺激神经组织。TMAES原理如图1所示。
图1 TMAES原理
Fig.1 Schematic diagram of TMAES
带电离子的振动速度v与超声声压p关系为
式中，ρ为生物组织密度；c为其中传播的声速。
假设静磁场的磁通密度为B，离子带电量为q，带电离子在超声波的作用下受到的洛伦兹力F为
生物组织产生的感应电流密度J为
式中，n+和n−分别为正、负带电离子的浓度；u+和u−分别为正、负带电离子的位移。将式（2）代入式（3）中可得
生物组织的电导率σ为
计算式（4）和式（5）可得
即
感应电流随超声频率振动，电流方向与超声波传播方向、磁场方向互相垂直。超声功率密度W为
引起的感应电场强度E为
式中，感应电场强度E为TMAES的刺激条件。
1.2 皮层神经元模型
本文构建了具有真实形态的皮层神经元模型，模拟TMAES对皮层刺激的单细胞和神经元集群反应。模型的建立参考了蓝脑计划所公布的相关电生理实验数据[18]，包含真实的神经元三维结构参数以及对应的生物物理学信息。建模过程中，神经元三维结构参考皮层三维追踪数字重建所获得的原始神经元形态。图2为单个皮层神经元三维结构。
图2 单个皮层神经元三维结构
Fig.2 Three dimensional structure of single cortical neuron
脑皮质内绝大多数神经元的轴突结构中有髓鞘生成，有髓轴突具有较低的电刺激激活阈值。为了更加真实地模拟皮层神经元中的轴突结构，本文在保留神经元原始几何结构的基础上对模型的轴突段添加髓鞘。在仿真过程中，将Ranvier节点作为分界点，用有髓节间段替换轴突中间段的空白分支。每一节间段的长度设置为该段纤维直径的100倍[19]。由于轴突末端分支的髓鞘相对较短，将节间段长度设置为对应段纤维直径的70倍。而轴突终端的分支纤维直径过小，不设置髓鞘。
采用基于特征的多目标优化方法[20]对神经元胞体、基底树突、顶端树突和轴突起始段的离子通道进行拟合，添加瞬时钠、持续钠、瞬时钾、持续钾、M电流、H电流、高压活化钙、低压活化钙、A型钾、D型钾、随机钾和SK钙活化钾等H-H型离子通道。模拟有髓节段的离子通道时，将钙通道移除，并使瞬时钠通道的电导加倍，其他离子通道设置与轴突起始段中相同。
1.3 神经网络的构建
为了分析TMAES对神经元集群放电活动的影响，选取36个皮层神经元模型进行神经网络构建，通过控制特定的离子通道活性，实现不同的神经元功能。所构建的神经网络共包含三种细胞类型：规则放电（Regular Spiking, RS）锥体神经元、快速尖峰中间神经元（Fast Spiking-Interneurons, FS-INs）和低阈值尖峰中间神经元（Low Threshold Spiking- Interneurons, LTS-INs），其中，FS-INs和LTS-INs各6个[21]，其余为RS锥体神经元。将所有细胞固定在一个6×6的正方形网格空间中，搭建脑皮层初级体感新皮质的神经网络模型。用LFPs信号采集模拟工具LFPsim[22]采集神经网络的电信号，图3为神经网络模型。
图3 神经网络模型
Fig.3 Schematic diagram of neural network
基于上述三种类型细胞的突触耦合实验数据[23]，在所构建的网络中神经元间连接比例设置为随机均匀选择，并在仿真过程中保持不变，见表1～表3。代表AMPA的突触兴奋性连接用0.1ms上升和3ms衰减速率时间常数来模拟，代表GABAA的突触抑制性连接用0.5ms上升和8ms衰减速率时间常数来模拟。神经元两两之间的连接方式如图4所示，包括抑制-兴奋类、抑制-抑制类、抑制-长时抑制类三种连接方式。箭头指向代表网络中不同位置神经元间的信息传递方向。
表1 兴奋性RS连接
Tab.1 Excitatory RS connections
参 数 RS→RS RS→FS RS→LTS 连接率(%) 5 44 58 所有可能连接 1 262 214 217 随机选择连接 75 92 124 最大突触电导/μS 0.001 635 0.001 316 0.000 995 12
表2 抑制性FS连接
Tab.2 Inhibitory FS connections
参 数 FS→RS FS→FS FS→LTS 连接率(%) 43 55 35 所有可能连接 215 34 35 随机选择连接 94 15 13 最大突触电导/μS 0.008 0.023 421 0.1
表3 抑制性LTS连接
Tab.3 Inhibitory LTS connections
参 数 LTS→RS LTS→FS LTS→LTS连接率(%) 36 62 3 所有可能连接 217 35 31 随机选择连接 77 21 0 最大突触电导/μS 0.105 9 0.002 714 —
图4 神经元连接方式
Fig.4 Neuronal connection mode
2 结果与分析
本文建立了可以反映脑皮层三维结构和生物物理特性的真实神经元模型，并搭建了神经网络，模拟了阶跃电流刺激、突触激活、TMAES三种刺激方式下神经元LFPs信号的响应状态，采用时频分析方法对不同刺激参数下阵列电极采集的神经网络LFPs数据进行处理，分析了TMAES感应电流强度与调制频率对LFPs信号能量分布的影响。
2.1 TMAES下神经元模型的电活动响应
2.1.1 不同刺激方式效果对比
神经元的放电模式反映了神经系统的信息编码过程，在大脑信息处理过程中发挥了重要作用。为了研究TMAES对神经元的刺激效果，采用外加阶跃电流刺激、自身突触激活、TMAES三种不同方式进行对比。在神经元胞体处施加阶跃电流刺激，如图5a所示，在突触激活的情况下相对应区域的突触被激活，如图5b所示。
图5 不同刺激方式
Fig.5 Different stimulation modes
规定所施加的阶跃电流刺激方式的电流幅值为0.5nA，并对模型施加相同幅值的TMAES感应电流；调节神经元突触设置，使得TTPC2类突触处于激活状态。不同刺激方式下的膜电压曲线如图6所示，图6a、图6c、图6e分别为阶跃电流刺激、TTPC2突触激活、TMAES三种刺激方式产生的膜电压数据，图6b、图6d、图6f为各自对应的动作电位局部放大曲线。由图6g中对比可以发现，三种方式均可以引发神经元放电活动，所产生的动作电位幅值差别不明显，相比于突触激活方式，阶跃电流刺激和TMAES刺激方式的动作电位发放频率逐渐趋于稳定。
图6 不同刺激方式下的膜电压曲线
Fig.6 Membrane potential curves under different stimulation modes
为了进一步对比外加阶跃电流刺激和TMAES两种刺激方式对神经元放电活动产生的影响，分别用幅值为0.5nA的阶跃电流和等幅值的TMAES感应电流刺激神经元，通过短时傅里叶变换对所采集的LFPs数据进行时频分析，所得结果如图7所示。对比图7a、图7b可得，TMAES刺激方式下超过2Hz的神经元放电活动能量带出现要晚于阶跃电流刺激方式，但两者所到达到的能量强度相近。
图7 不同刺激方式下LFPs信号能量时频分布
Fig.7 Time frequency distribution of LFPs signal energy under different stimulation modes
2.1.2 不同深度处刺激效果对比
神经元的放电活动受到自身形态特征和结构的影响，本节对TMAES作用下皮层真实神经元不同深度处的电活动响应特性进行分析。对模型施加TMAES感应电流，并沿神经元胞体主轴线上下移动测量点位置，采用线源记录方式模拟计算在神经元不同位置处的LFPs数值，得到的LFPs数据如图8和图9所示。
图8 不同深度LFPs
Fig.8 LFPs at different depths
图9 不同深度LFPs信号能量时频分布
Fig.9 Time frequency distribution of LFPs signal energy in different depths
如图8和图9所示，对神经元模型施加TMAES刺激后在不同深度采集到的LFPs信号能量并不相同，其中，−0.1～0.1mm处LFPs信号幅值最大、能量最强，其他深度能量较低，这与胞体处神经元分支密集、电活动最为活跃的情况相吻合。神经元胞体附近测量到的LFPs信号能量时频分布如图10所示，由图可见，在神经元胞体附近记录到的LFPs信号相对稳定，随时间变化呈现出规则放电状态，且能量相对集中于0～50Hz的频率范围内。
图10 胞体附近LFPs信号能量时频分布
Fig.10 Time frequency distribution of LFPs signal energy near cell body
2.2 TMAES下神经网络的电活动响应
2.2.1 不同深度处刺激效果对比
为了研究TMAES对神经元网络放电活动的影响，本文建立了脑皮层神经元网络模型，以模拟皮质中存在的小簇共激活神经元。本节主要对TMAES中超声调制频率和感应电流强度两个指标的改变对皮层中小簇共激活神经元放电模式的影响进行分析，设置调制频率为50Hz，加入幅值为0.5nA的感应电流，利用线源记录方法测量神经网络不同位置的LFPs数据，如图11和图12所示。
图11 神经网络不同深度处LFPs
Fig.11 LFPs at different depths of neural network
对神经元网络模型施加TMAES刺激，在不同深度采集到的LFPs信号能量分布并不相同，其中，1.2mm处能量最为集中，这与网络在该深度处神经元胞体密集，电活动相对规则的状态相吻合；其他深度处能量分布较为分散，与网络其他深度神经元结构分支简单、电信号传导随空间衰减有关。
图12 神经网络不同深度LFPs信号的能量时频分布
Fig.12 Energy time frequency distribution of LFPs signals with different depths in neural network
2.2.2 不同刺激参数下的刺激效果对比
利用虚拟阵列电极在6×6的神经网络中进行测量，模拟皮层的体外多电极阵列记录，电极间距为100μm。神经网络的多电极阵列记录如图13所示。
图13 神经网络的多电极阵列记录
Fig.13 Multi electrode array recording of neural network
为了研究TMAES不同刺激参数对神经网络电活动响应的影响，选用超声基波频率为0.5MHz，设置调制频率为50Hz，加入不同幅值的TMAES感应电流信号，所得不同感应电流强度下LFPs信号能量时频分布如图14所示。将施加的TMAES感应电流固定为0.5nA，改变调制频率，所得不同调制频率下LFPs信号能量时频分布如图15所示。
图14 不同感应电流强度下LFPs信号能量时频分布
Fig.14 Time frequency distribution of LFPs signal energy under different induced current intensity
图15 不同调制频率下LFPs信号能量时频分布
Fig.15 Time frequency distribution of LFPs signal energy under different modulation frequency
由图14和图15可见，随着感应电流强度的增加，LFPs能量分布在时域范围内经历了逐渐增强而 后减弱的过程，而在频域范围内基本稳定在0～4Hz的范围内；随着调制频率的增加，LFPs能量分布在时域和频域范围内均经历了逐渐增强而后减弱的过程。这表明TMAES可对神经网络放电活动产生促进与抑制两种效果，改变刺激参数可以实现对生物体神经活动的调节。
3 结论
为了分析TMAES中刺激参数对神经元及神经元集群放电活动的影响，本文建立了具有真实形态结构的皮层单神经元和神经网络的计算模型，模拟了不同刺激方式和TMAES方式下不同感应电流密度、不同调制频率对神经元网络电特性的影响。通过时频联合分析方法对LFPs仿真数据进行处理，结果表明，与神经元自身突触激活方式和外加电流刺激方式相比，TMAES可以实现对神经元放电活动的有效激活。进一步分析不同TMAES刺激参数下神经网络的局部场电位能量分布，发现随着调制频率和感应电流强度的增加，LFPs信号在时域的能量分布呈现先增大后减小的趋势。这表明TMAES可以实现促进和抑制神经元放电活动两种调节效果。
本文建立的神经网络计算模型为研究电磁场对生物体影响提供了一种有效的补充工具，通过模拟具有真实形态的单神经元及小簇共激活神经元集群的放电活动，可以更为详细地刻画脑皮层在外加电磁场影响下的电特性变化。对于TMAES下神经元网络放电活动的仿真分析，可以将多个层次的实验数据进行统一，有助于深入理解TMAES的神经作用机理，推动无创式脑调控技术的进一步发展。
[1] 伊国胜, 王江, 魏熙乐, 等. 无创式脑调制的神经效应研究进展[J]. 科学通报, 2016, 61(8): 819-834.
Yi Guosheng, Wang Jiang, Wei Xile, et al. Deve- lopments of neural effects induced by noninvasive brain modulation[J]. Chinese Science Bulletin, 2016, 61(8): 819-834.
[2] 刘婧, 刘国强. 电粒子成像方法及其正问题数值研究[J]. 电工技术学报, 2020, 35(22): 4621-4626.
Liu Jing, Liu Guoqiang. Electrical particles imaging method and numerical study on the forward pro- blem[J]. Transactions of China Electrotechnical Society, 2020, 35(22): 4621-4626.
[3] 黎国锋, 邱维宝, 钱明, 等. 超声神经调控技术与科学仪器[J]. 生命科学仪器, 2017, 15(1): 3-8.
Li Guofeng, Qiu Weibao, Qian Ming, et al. Ultrasonic neuromodulation technology and scientific equipment[J]. Life Science Instruments, 2017, 15(1): 3-8.
[4] 熊慧, 王玉领, 付浩, 等. 一种应用于经颅磁刺激脉冲宽度可调的节能型激励源[J]. 电工技术学报, 2020, 35(4): 679-686.
Xiong Hui, Wang Yuling, Fu Hao, et al. An energy efficient excitation source for transcranial magnetic stimulation with controllable pulse width[J]. Transa- ctions of China Electrotechnical Society, 2020, 35(4): 679-686.
[5] Rezayat E, Toostani I G. A review on brain stimu- lation using low intensity focused ultrasound[J]. Basic & Clinical Neuroscience, 2016, 7(3): 187-194.
[6] 刘素贞, 魏建, 张闯, 等. 基于FPGA的超声信号自适应滤波与特征提取[J].电工技术学报, 2020, 35(13): 2870-2878.
Liu Suzhen, Wei Jian, Zhang Chuang, et al. Adaptive filtering and feature extraction of ultrasonic signal based on FPGA[J]. Transactions of China Electro- technical Society, 2020, 35(13): 2870-2878.
[7] 郭磊, 刘东钊, 黄凤荣, 等. 基于突触可塑性的自适应脉冲神经网络在高斯白噪声刺激下的抗扰功能研究(英文)[J]. 电工技术学报, 2020, 35(2): 225-235.
Guo Lei, Liu Dongzhao, Huang Fengrong, et al. Research on disturbance rejection of adaptive spiking neural network based on synapticplasticity under white gaussian noise[J]. Transactions of China Elec- trotechnical Society, 2020, 35(2): 225-235.
[8] 刘运华, 张波, 谢帆, 等. 多尺度和多物理场的电力电子变换器建模方法初探[J]. 电力系统自动化, 2020, 44(16): 61-69.
Liu Yunhua, Zhang Bo, Xie Fan, et al. Preliminary study on modeling methods with multiscale and multiphysics for power electronic converters[J]. Automation of Electric Power Systems, 2020, 44(16): 61-69.
[9] 赵禹, 杨仕友. 无线电能传输系统的频率跟踪技术与控制方法[J]. 电机与控制学报, 2020, 24(9): 22-29.
Zhao Yu, Yang Shiyou. Frequency tracking and controlling of wireless power transfer system[J]. Electric Machines and Control, 2020, 24(9): 22-29.
[10] Li Dandan, Qiao Zhenyang, Yang Na, et al. Study on vector magnetic properties of magnetic materials using hybrid hysteresis model[J]. CES Transactions on Electrical Machines and Systems, 2019, 3(3): 292-296.
[11] 赵靖英, 周思诺, 高佳雄. 智能电动机保护器的电磁辐射定量评估方法研究[J]. 电机与控制学报, 2017, 21(6): 34-43.
Zhao Jingying, Zhou Sinuo, Gao Jiaxiong. Research on quantitative assessment method of electromagnetic radiation of intelligent motor protector[J]. Electric Machines and Control, 2017, 21(6): 34-43.
[12] Norton S J. Can ultrasound be used to stimulate nerve tissue?[J]. BioMedical Engineering OnLine, 2003, 2(1): 6-14.
[13] 李慧雨, 周晓青, 张顺起, 等. 基于磁声耦合效应的聚焦电刺激方法的初探[J]. 生物医学工程研究, 2015, 34(4): 201-206.
Li Huiyu, Zhou Xiaoqing, Zhang Shunqi, et al. Approach for focused electric stimulation based on the magneto-acoustic effect[J]. Journal of Biomedical Engineering Research, 2015, 34(4): 201-206.
[14] Yuan Yi, Li Xiaoli. Theoretical analysis of trans- cranial hall-effect stimulation based on passive cable model[J]. Chinese Physics B, 2015, 24(12): 373-378.
[15] 刘世坤, 张鑫山, 刘志朋, 等. 经颅磁声耦合电刺激技术应用于小鼠的实验研究[J]. 生物医学工程研究, 2018, 37(1): 11-15.
Liu Shikun, Zhang Xinshan, Liu Zhipeng, et al. Experimental study in mice on the technology of transcranial magneto-acoustic coupling electrical stimulation[J]. Journal of Biomedical Engineering Research, 2018, 37(1): 11-15.
[16] 张帅, 崔琨, 史勋, 等. 经颅磁声电刺激参数对神经元放电模式的影响分析[J]. 电工技术学报, 2019, 34(18): 3741-3749.
Zhang Shuai, Cui Kun, Shi Xun, et al. Effect analysis of transcranial magnetic-acousticao-electrical stimu- lation parameters on neural firing patterns[J]. Transactions of China Electrotechnical Society, 2019, 34(18): 3741-3749.
[17] 张帅, 史勋, 尹宁, 等. 基于H-H神经元模型的经颅磁声刺激对神经元放电活动的影响[J]. 高电压技术, 2019, 45(4): 122-128.
Zhang Shuai, Shi Xun, Yin Ning, et al. Effects of transcranial magneto-acoustical stimulation on neuronal firing activities based on H-H neuron model[J]. High Voltage Engineering, 2019, 45(4): 122-128.
[18] Michael W R, Costas A, Rodrigo P, et al. A biophysically detailed model of neocortical local field potentials predicts the critical role of active membrane currents[J]. Neuron, 2013, 79(2): 375-390.
[19] Pinto D J, Jones S R, Kaper T J, et al. Analysis of state-dependent transitions in frequency and long- distance coordination in a model oscillatory cortical circuit[J]. Journal of Computational Neuroence, 2003, 15(2): 283-298.
[20] Shaul D. A novel multiple objective optimization framework for constraining conductance-based neuron models by experimental data[J]. Frontiers in Neuroscience, 2007, 1(1): 7-18.
[21] Toledo R M, Wang Yun, Gupta A, et al. Interneurons of the neocortical inhibitory system[J]. Nature Reviews Neuroscience, 2004, 5(10): 793-807.
[22] Einevoll G T. Modeling the spatial reach of the LFP[J]. Neuron, 2011, 72(5): 859-872.
[23] Gibson J R, Beierlein M, Connors B W. Functional properties of electrical synapses between inhibitory interneurons of neocortical layer 4[J]. Neurophysiol, 2005, 93(10): 467-480.