王智慧
(北方工业大学 理学院, 北京 100144)
癫痫是一种由神经元群异常超同步放电引起的慢性神经系统疾病[1],根据癫痫发作的特点可分为局灶性癫痫发作和全身性癫痫发作。失神发作是全身性癫痫的一种,以慢性、反复发作为特征, 表现为短暂的意识障碍,严重影响患者的健康和生活质量[2-3]。在脑电图(EEG)记录中,典型的失神发作以2~4 Hz的棘慢波(SWD)为特征[4-5],而非典型的失神发作频率约为1~2 Hz[6]。虽然有研究表明失神癫痫发作的解剖学原因包含皮质丘脑环路异常的相互作用[7],但是其发病机制不完全清楚。失神发作的病理特征为拟合生理实验建模提供了可能, 因而从动力学角度详细描述癫痫发作过程中癫痫样放电的演变过程是必要的,可以为临床提供理论指导。
为了探究癫痫的病因和发病机制,人们提出了很多基于电生理实验的数学模型,其中平均场模型是最受欢迎的模型之一。平均场模型比神经元网络模型的优势在于,其可以预测神经元集群的大规模属性,并估计神经元集群间的连接强度。此外,由于平均场模型表达简单,数值和分析结果更容易得到,更有利于揭示癫痫发作相关的内在机制。针对失神发作,ROBINSON等[8]提出了具有一定生理基础的皮质丘脑平均场模型, 可以再现临床观察到的棘慢波。MARTEN等[4]和DEEBA等[9]基于平均场模型,复现出更为复杂的多棘慢波(m-SWD)。近年来,改进的平均场模型层出不穷,关于癫痫发作动态性质的研究也在不断推进。例如,加入从丘脑网状核到丘脑中继核通路的时滞[10],将抑制性神经元集群分为两类[11]等。而揭示改进模型中的内在机理对于理解真正的癫痫发作起着重要作用。BREAKSPEAR等[12]将皮层作为电活动波传播的媒介,通过大脑平均场非线性模型的动力学分岔,结合关键的生理过程,解释了全身性癫痫发作的重要特征。CHEN等[13]证明丘脑中的前馈抑制,可能参与控制失神发作,特别是增加前馈抑制的兴奋性输入,可以显著抑制失神发作期间的放电活动。FAN等[14]改进皮质丘脑平均场模型,提出一种单脉冲交替复位刺激,通过计算神经元集群的平均放电率,给出刺激减弱癫痫发作的动力学解释。
生理学研究证实,丘脑网状核(TRN)向中继核(SRN)的抑制性投射被认为是引起失神发作的重要病理因素[15-16]。一些学者已经开始利用这个因素来进行癫痫的理论研究。例如:FAN等[17]通过改进皮质丘脑环路网络模型,针对中继核所在递归兴奋性环路,发现丘脑如果不能够有效地中继兴奋性皮层传递的信息,会导致2~4 Hz棘慢波放电的产生。 CHEN等[10]利用平均场模型,证实GABAB受体在TRN 中的慢突触动力学可以诱导失神发作,表明皮质丘脑系统内的异常改变会引发失神性癫痫。 WANG等[18]针对丘脑中继核团,探究其关联的前馈抑制和反馈抑制微环路,发现两个微环路可以同时调控失神发作,而在抑制失神发作方面,丘脑反馈抑制强于前馈抑制。大脑皮层与丘脑之间相互连接构成复杂网络,而丘脑网状核是丘脑中重要的神经核团, 其参与癫痫的发生和发展过程。分析丘脑网状核与丘脑中继核的相互关系对癫痫的影响机制,阐述癫痫发作及动力学转迁机理,是现在亟须解决的问题。
空间拓展多室耦合模型,是模拟宏观失神发作的良好模型,且应用广泛。TAYLOR等[19]和SINHA等[20]基于空间扩展神经场模型,发现棘慢波放电是由兴奋性皮质区域内的局部活动引起,并通过皮质丘脑网络快速传播。CAO等[21]通过建立单向流动的耦合神经场模型,阐述了系统随某一参数变化的动力学转迁行为,发现双室不同于单室,其动力学更丰富。WANG等[22]基于单向流动的耦合神经场模型,在病态区域丘脑网状核中加入不同的刺激脉冲,通过定性分析,找到控制失神发作的最优控制策略。AHMADIZADEH等[23]利用双室耦合的Jansen-Rit神经质量模型,通过改变系统输入、耦合强度和网络结构,得到棘慢波放电状态。双室耦合模型是最基本的皮质丘脑网络模型,结构简单却更具生理意义,并能对更大的网络提供理论支撑。以上研究虽然是基于神经场或神经质量模型,但对研究癫痫的发病机制具有重要的参考意义。因此,本文将建立单向流动的双室皮质丘脑平均场模型,通过改变GABAB受体的慢动力学特性,探究系统放电模式的动态转迁。
在ROBINSON等[8]提出的确定性模型的基础上,建立单向流动的双室皮质丘脑平均场模型。如图1所示,此模型可以描述丘脑和皮质神经元集群的动态变化,其中皮层子系统包括兴奋性锥体神经元群体(EPN)和抑制性中间神经元群体(IIN),丘脑子系统由特异性中继核(SRN)和丘脑网状核(TRN)组成。黑色直线箭头代表GABAA受体调控的快抑制性投射,黑色虚线箭头代表GABAB受体调控的慢抑制性投射,紫色箭头线代表谷氨酸受体调控的兴奋性投射。失神发作是一种典型的全身性癫痫,其动态活动被认为是在整个大脑同时发生的,因此神经元集群的动力学变量只与时间有关。
一般来说,神经元群体放电的过程会产生脉冲场,这里用模型中的平均放电率来模拟。 脉冲场通过阻尼波传播,在突触后神经元集群中,脉冲场的传入信号产生光滑的树突反应,被过滤产生膜电压,膜电压通过S形曲线转变为神经元集群的平均放电率,整个过程是一个循环的过程。
平均场模型可以简单而有效地研究神经元集群的宏观动态特征。其中,神经元集群的平均放电率Rx(t)可被平均膜电位Vx(t)表示为:
Rx(t)=F(Vx(t))=
(1)
(2)
微分算子Dαβ可以被表示为:
(3)
式(2)和式(3)中:Dαβ表示输入脉冲通过树突的滤波效果,α和β表示细胞体对输入信号的反应时间,φy(t)表示神经元集群y产生的传入脉冲率。vxy表示神经元集群y作用于x的输入脉冲的突触强度。由于皮质内突触数量与连接成正比,可忽略皮层抑制性膜电位Vi,利用Ve=Vi和Re=Ri将模型进行简化。
由于IIN、TRN和SRN 3个核团的轴突较短, 无法有效传播脉冲场,因此用以下函数近似表示:
φz(t)≈Rz(t)=F(Vz(t)) ,
(4)
式中:z=i,s,r。大脑皮层兴奋性轴突场具有有效的传播效应,因此φe可单独被表示为:
Re(t)=F(Ve(t)) ,
(5)
结合上面的描述,利用方程(1)—(3)和(5)模拟兴奋性锥体神经元集群的动力学形态,利用方程(1)—(4)模拟其他神经元集群的动力学形态。在单向传递过程中,只考虑相同子系统间的传递投射。 因此,将上述方程改写为一阶微分方程组,模型可以描述如下:
双室耦合模型足以解释癫痫的一些定性行为。 为了简单起见,在目前的工作中不考虑双室耦合间的传输延迟,仅考虑丘脑内的传输延迟τ,来描述受体GABAB诱导的慢动力学特性。φn1(t)和φn2(t)表示丘脑的外部输入。此外,耦合双室间的距离,有三种连接形式,分别为短程连接、远程连接和远距离兴奋性连接[21]。因为仅考虑相同子系统间的传递投射,故而模型采用短程连接。即:假设SRN1对TRN1的投射强度为-vsr,那么SRN1对TRN2的投射强度为-vsr/3,其他情形类似。各方程的连接形式在图1中已经详细地给出。
模型中的所有参数是根据实验数据所估计,参数也是根据已有建模研究所取值,均在电生理实验意义范围之内[10,13]。MATLAB R2013a(MathWorks, USA)仿真环境非常灵活,用于执行模型中的所有数值计算。采用dde23求解模型中的时滞微分方程组。
研究表明,TRN中GABAB受体所产生的慢动力学特性是诱导癫痫样放电的病理因素之一。这一结论在动物实验和皮质丘脑网络的生物物理模型中都有所表现[10-11]。为了探索这种机制是否也适用于耦合平均场模型,分别对TRN-SRN通路的抑制性耦合强度-vsr和延迟参数τ进行一维状态分析,来探究GABAB受体的慢动力学特性对室I中的癫痫样放电的影响机制。首先讨论不同的耦合强度-vsr下,室I放电率φe的变化情况,对应的状态分析、主频分析和典型的时间序列如图2(a)—(f)所示。当耦合强度-vsr较弱时,TRN1对SRN1的抑制性也很弱,故而无法有效抑制SRN1的活性,SRN1的活性水平迅速达到最大,SRN1的高活性驱动皮层兴奋性神经元集群放电达到饱和状态(区域I)。随着耦合强度-vsr的增加,TRN1对SRN1的抑制性开始发挥作用,当时滞足够大时,系统会经历两种不同的放电模式:第一种是棘慢波放电(区域III,简称SWD 放电),其波形的特征是每个周期内有两个极大值和两个极小值,这种波形是癫痫样放电的代表形式之一,并且已经在临床患者癫痫发作时的脑电图记录中被观察;第二种是简单振荡(区域II),其波形的特征是每个周期内有一个极大值和一个极小值。当耦合强度-vsr过强时,SRN1的活性完全被TRN1所抑制,系统进入低放电状态(区域V),即稳定状态,没有任何振荡形式。
图2 室I中耦合强度-vsr对癫痫样放电的影响(a) 状态分析;(b) 主频分析;(c) 饱和状态;(d) 棘慢波振荡;(e) 简单振荡;(f) 低放电状态Fig. 2 Effect of coupling strength -vsr on epileptic discharges in compartment I (a) State analysis; (b) Dominant frequency analysis; (c) High saturation state;(d) SWD oscillation; (e) Simple oscillation; (f) Low firing state
影响GABAB受体的慢动力学特性的另外一个因素是延迟参数τ。模型中引入时滞τ,探究时滞τ的一维状态分析(图3(a))、主频分析(图3(b))以及对应的时间序列(图3(c)—(f))。如图3(a)所示,当时滞τ较小时,TRN1产生的抑制性受体GABAA和GABAB同时作用在SRN1上,导致SRN1的活性减弱,系统处于简单振荡模式(区域II)。随着时滞τ逐渐增加,GABAB受体诱导的慢抑制性发挥作用,TRN1对SRN1的抑制性开始减弱,系统经历两种病态的放电模式:第一种是棘慢波放电,第二种是2-棘慢波放电(区域IV,简称2-SWD放电)。2-SWD波形的特征是每个周期内有3个极大值和3个极小值,这种波形也包含在癫痫样放电模式中。当时滞τ足够大时,SRN1开始仅受到抑制性受体GABAA的作用,也就是说TRN1对SRN1的抑制性作用在整体上较弱,SRN1的活性水平较高,这种较高水平的活性驱动皮层兴奋性神经元集群放电达到饱和状态(区域I)。从解剖学上来说, SRN1神经元集群通过GABAA和GABAB两种不同时间尺度下的受体介导的抑制性通路来接收TRN1神经元集群发出的信号。在一定条件下,这两种受体在不同时间瞬间产生的双重抑制性为SRN1产生多个峰放电提供了可能,SRN1进而影响皮层神经元集群放电形式,从而导致了棘慢波或多棘慢波的产生。
图3 室I中时滞τ对癫痫样放电的影响(a) 状态分析;(b) 主频分析;(c) 简单振荡;(d) 棘慢波振荡;(e) 2-棘慢波振荡;(f) 饱和状态.Fig. 3 Effect of time delay τ on epileptic discharges in compartment I (a) State analysis; (b) Dominant frequency analysis; (c) Simple oscillation; (d) SWD oscillation; (e) 2-SWD oscillation; (f) Saturation state
时滞τ和耦合强度-vsr均可以诱导癫痫样放电, 因此关于两个变量,分别探究了二维状态分析(图4(a))和对应的主频分析(图4(b))。图4(a)被分成5个区域,区域I为饱和状态,当耦合强度-vsr较弱时,增大时滞τ,此时TRN1无法有效抑制SRN1的活性,SRN1的兴奋性驱动皮层神经元集群放电达到饱和状态。
图4 室I中时滞τ和耦合强度-vsr共同诱导癫痫样放电(a) 二维状态分析;(b) 对应的主频分析Fig. 4 Time delay τ and coupling strength -vsr induce epileptic discharges jointly in compartment I(a) Two-dimensional state analysis; (b) Corresponding dominant frequency analysis
随着耦合强度-vsr的增强,当时滞τ较小时,SRN1的活化水平较弱,系统出现简单振荡状态(区域II);时滞τ逐渐增加,GABAB受体诱导的慢抑制性发挥作用,TRN1对SRN1的抑制性减弱, SRN1的活化水平开始升高,此时系统出现棘慢波放电(区域III)和2-棘慢波放电(区域IV)两种病态区域。 当耦合强度-vsr增加到一定程度,SRN1的活性完全被抑制,皮层神经元集群达到静息状态,即低放电状态(区域V)。 时滞τ和耦合强度-vsr的二维状态分析及其对应的主频分析能够充分展示GABAB受体的慢动力学特性诱导癫痫样放电。
采用单向流动的双室皮质丘脑平均场模型,从图1结构图可以发现,室II除了自身结构的运转之外还受到来自室I的兴奋性或抑制性作用。这也决定了GABAB受体的慢动力学对室II 中癫痫样放电的影响不同于室I。基于此,本节集中讨论GABAB受体的慢动力学对室II 中癫痫样放电的影响机制。首先,分别对TRN-SRN通路的抑制性耦合强度-vsr和延迟参数τ进行一维状态分析。如图5(a)所示,当耦合强度-vsr较弱时,TRN2对SRN2的抑制性较弱,SRN2的活化水平较高,此时两种额外的作用,即:SRN1对TRN2的兴奋性作用,以及TRN1对SRN2的抑制性作用促使SRN2的高活化水平被抑制,故室II没有出现饱和状态,而是低放电状态(V)。
图5 室II中耦合强度-vsr对癫痫样放电的影响(a) 状态分析;(b) 主频分析Fig. 5 Effect of coupling strength -vsr on epileptic discharges in compartment II(a) State analysis;(b) Dominant frequency analysis
随着耦合强度-vsr的增大,当时滞τ足够大时,系统出现棘慢波放电(III)。-vsr继续增大,意味着TRN2对SRN2的抑制性继续增强,系统由癫痫样放电状态转变为简单振荡状态(II)。当耦合强度-vsr足够大时,SRN2的活性被完全抑制,系统进入低放电状态(V)。比较图2(a)和5(a),不难发现,室II的癫痫样放电区域明显小于室I。
时滞τ的一维状态分析和主频分析如图6(a)和6(b)所示。当时滞τ较小时,TRN2产生的抑制性受体GABAA和GABAB几乎同时作用在SRN2上,导致SRN2的活性减弱,系统处于简单振荡模式(区域II)。 随着时滞τ逐渐增加,TRN2对SRN2的抑制性开始减弱,由于TRN1对SRN2也存在抑制作用,因此简单振荡模式存在的时间更长。继续增大时滞τ,GABAB受体诱导的慢抑制性发挥作用,系统出现两种病态放电模式,棘慢波放电模式(区域III)和2-棘慢波放电模式(区域IV)。当时滞τ足够大时,TRN2对SRN2的抑制性相对系统整体而言较弱,SRN2的活性水平升高。然而此时TRN1对SRN2的抑制性和SRN1对TRN2的兴奋性也逐渐增强,使得SRN2的高活化水平得到抑制,这种抑制性对系统产生的影响比时滞τ的更直观,作用效果更明显,因而系统出现低放电状态(区域V)。比较图3(a)和6(a),不难发现,室II的癫痫样放电区域明显小于室I。
图6 室II中时滞τ对癫痫样放电的影响(a) 状态分析;(b) 主频分析Fig. 6 Effect of time delay τ on epileptic discharges in compartment II (a) State analysis; (b) Dominant frequency analysis
关于两个变量,分别探究了二维状态分析(图7(a))和对应的主频分析(图7(b))。图7分为4部分,分别为低饱和放电区域(V)、简单振荡区域(II)、棘慢波放电区域(III)和2-棘慢波放电区域(IV),每个区域用不同的颜色所表示。如图7(a)所示,当耦合强度很小,时滞很大时,SRN2的活性被抑制,系统出现低饱和状态。随着耦合强度增大,当时滞很小时,系统出现简单振荡状态,时滞逐渐增加, TRN2对SRN2的抑制性减弱,但TRN1对SRN2的抑制性,使得简单振荡区域增大,当GABAB受体诱导的慢抑制性发挥作用,系统出现癫痫样放电模式,即:棘慢波放电模式和2-棘慢波放电模式。当耦合强度很大,时滞很小时,SRN2的活性被TRN2所抑制,系统出现低放电状态。比较图4(a)和图7(a),发现两室均出现癫痫样放电状态,区别是室I比室II 拥有更加丰富的动力学行为,且室I的棘慢波及2-棘慢波放电区域明显大于室II。
图7 室II中时滞τ和耦合强度-vsr共同诱导癫痫样放电(a) 二维状态分析;(b) 对应的主频分析Fig. 7 Time delay τand coupling strength -vsr induce epileptic discharges jointly in compartment II(a) Two-dimensional state analysis;(b) Corresponding dominant frequency analysis
癫痫发作的根源在于大脑功能网络中兴奋性和抑制性信息流的不平衡,而丘脑网状核中GABAB受体的慢动力学特性的改变正是信息流不平衡的一种表现。基于此,本文建立了单向流动的双室皮质丘脑平均场模型,来探究这种不平衡对癫痫样放电的影响。首先,TRN在一定程度上可以抑制SRN的兴奋性,通过模拟发现,当时滞达到一定水平时,SRN活性由强到弱所对应的动力学状态依次为:高饱和状态、2-棘慢波放电状态、简单振荡状态和低放电状态。且在同等参数取值下,室I出现全部的动力学状态,室II虽然也出现棘慢波和2-棘慢波放电状态,但在区域面积上很明显少于室I。通过二维状态的分析,无论在室I还是室II中,发现较长时滞τ均更容易诱导棘慢波及多棘慢波的放电状态,而由于室II除了自身结构的运转之外,还受到来自室I的兴奋性或抑制性作用,因此室I表现出更加丰富的动力学状态,且室II的癫痫样放电区域明显少于室I。与以往的模型相比,本文所建立的双室皮质丘脑平均场模型,能模拟出多棘慢波放电模式,动力学状态更加丰富。这些结果为更好地构建皮质丘脑网络模型提供了理论参考。