鲁尔科夫映射

来自集智百科 - 复杂系统|人工智能|复杂科学|复杂网络|自组织
跳到导航 跳到搜索

鲁尔科夫映射(Rulkov map),又称鲁尔科夫模型,是用离散时间迭代描述神经元放电活动的一类低维非线性模型。映射在这里指从当前状态计算下一步状态的更新规则。模型用快变量表示膜电位样的快速变化,用慢变量表示恢复或适应过程。两者相互作用,可以产生静息、持续发放和簇放电,即一组快速锋电位与安静间隔交替出现的活动。锋电位(spike)是神经元快速产生的动作电位事件。

鲁尔科夫映射主要包括 2001 年的光滑版本和 2002 年引入显式锋电位重置的分段版本。Nikolai F. Rulkov 用前者研究不规则簇节律如何在耦合后变得规则,用后者描述簇内单次锋电位及其对输入的响应[1][2]。两个版本采用不同的快更新函数和参数定义,下面分别介绍其方程与机制。

这类模型属于现象学模型:用少量变量和更新规则保留选定的放电现象,适合研究神经动力学中的节律、同步与集体活动。变量与真实膜电位、离子电流之间的关系,需要根据所用版本另行标定。

基本信息

项目 内容
中文名称 鲁尔科夫映射、鲁尔科夫模型
英文名称 Rulkov map、Rulkov neuron model
数学形式 快变量与慢变量以不同速度更新的离散非线性映射
状态变量 快变量 [math]\displaystyle{ x_k }[/math]、慢变量 [math]\displaystyle{ y_k }[/math],[math]\displaystyle{ k }[/math] 为迭代序号
主要版本 2001 年光滑映射、2002 年分段映射及后续改进形式
建模对象 单个神经元,以及由这些单元构成的耦合网络
典型现象 静息、持续发放、周期或混沌簇放电,以及簇起止时刻的同步和慢节律趋于规则的变化

历史溯源

2001 年的研究考察耦合如何使性质不同、单独活动时簇节律不规则的神经元形成较规则的集体节律。Rulkov 用光滑二维映射表示各个单元,再加入平均场耦合,即让每个单元受到群体平均活动的影响[1]。2002 年的分段模型把快速上升和重置过程写入更新规则,使簇内的单次锋电位更容易识别,并分析了两个电耦合单元的同步。这里的电耦合以单元状态差作为相互作用的输入[2]。

2003 年,Andrey L. Shilnikov 与 Rulkov 进一步研究分段模型的混沌机制。他们分别分析快子系统和完整二维系统,说明慢变量的反馈如何改变整体轨道的稳定性[3]。后续研究又修改单元更新函数并加入化学突触,即通过神经递质传递信号的连接,以表示不同类型的皮层神经元及其网络活动[4]。这些扩展沿用了映射建模思路,具体能力由各版本的方程决定。

模型概述

光滑映射

光滑版本用两条更新规则连接快速活动与缓慢恢复过程[1]:

[math]\displaystyle{ \begin{aligned}x_{k+1}&=\frac{\alpha_{\mathrm R}}{1+x_k^2}+y_k,\\y_{k+1}&=y_k-\sigma_{\mathrm R}x_k-\beta_{\mathrm R}.\end{aligned} }[/math]

 

 

 

 

(1)

其中 [math]\displaystyle{ \alpha_{\mathrm R}\gt 0 }[/math] 控制快映射的非线性,[math]\displaystyle{ \sigma_{\mathrm R} }[/math] 和 [math]\displaystyle{ \beta_{\mathrm R} }[/math] 是小的正参数。常用示例取 [math]\displaystyle{ \sigma_{\mathrm R}=\beta_{\mathrm R}=0.001 }[/math],使慢变量 [math]\displaystyle{ y_k }[/math] 在典型轨道上比快变量 [math]\displaystyle{ x_k }[/math] 变化缓慢。这种更新速度的差异称为时间尺度分离。[math]\displaystyle{ x_k,y_k }[/math] 都是无量纲标量。

第一条规则中的有理函数具有单峰形状,负责产生快速振荡。第二条规则让快速活动反过来改变慢变量,从而影响下一阶段的快更新。光滑版本用连续形状的函数生成活动起伏,主要保留放电节律与快慢相互作用;生理动作电位的精细波形需要更具体的描述。

分段映射

分段版本为快速上升和重置分别设置分支,并允许输入进入快、慢更新[2]:

[math]\displaystyle{ \begin{aligned}x_{k+1}&=f_\alpha(x_k,y_k+b_k),\\y_{k+1}&=y_k-\mu(x_k+1-s_k).\end{aligned} }[/math]

 

 

 

 

(2)

参数 [math]\displaystyle{ \alpha\gt 0 }[/math] 控制快函数,[math]\displaystyle{ 0\lt \mu\ll1 }[/math] 控制慢过程的更新速度,常用示例取 [math]\displaystyle{ \mu=0.001 }[/math]。[math]\displaystyle{ b_k,s_k }[/math] 分别表示进入快、慢更新的输入,原论文记为 [math]\displaystyle{ \beta_k,\sigma_k }[/math]。本词条采用 [math]\displaystyle{ b_k,s_k }[/math],以区分式(1)中的固定参数 [math]\displaystyle{ \beta_{\mathrm R},\sigma_{\mathrm R} }[/math]。在通常的 [math]\displaystyle{ \alpha+v\gt 0 }[/math] 工作区间,快函数为:

[math]\displaystyle{ f_\alpha(u,v)=\begin{cases}\displaystyle\frac{\alpha}{1-u}+v,&u\le0,\\\alpha+v,&0\lt u\lt \alpha+v,\\-1,&u\ge\alpha+v.\end{cases} }[/math]

 

 

 

 

(3)

式(3)中,[math]\displaystyle{ u }[/math] 是当前快状态,[math]\displaystyle{ v }[/math] 是慢状态与快输入的合成量。第一支描述阈下变化,即尚未形成锋电位的活动。第二支使状态快速升高,第三支把状态重置为 [math]\displaystyle{ -1 }[/math],为下一次活动提供起点。第一支为有理非线性函数,后两支为常数,所以该模型是分段非线性映射。

研究无外部输入的单元时,常取 [math]\displaystyle{ b_k=0 }[/math]、[math]\displaystyle{ s_k=s }[/math] 为常数。后续版本会改变分支形状,或加入前一步状态以限制锋电位持续时间。复现某项研究时,输入位置、分支条件与重置规则应采用该研究给出的完整定义[4]。

动力学机制

快慢反馈与簇的起止

簇放电的关键是:快速活动改变慢变量,慢变量又改变快子系统能够维持的运动。为了分开观察这两种作用,可以在分析光滑映射时暂时固定 [math]\displaystyle{ y_k=\gamma }[/math],得到一维快子系统 [math]\displaystyle{ x_{k+1}=F(x_k;\gamma)=\alpha_{\mathrm R}/(1+x_k^2)+\gamma }[/math]。“冻结慢变量”是把慢变量当作快子系统的参数;完整系统仍按式(1)更新它。

由慢更新规则可知,当 [math]\displaystyle{ x_k\lt -\beta_{\mathrm R}/\sigma_{\mathrm R} }[/math] 时,慢变量上升;当 [math]\displaystyle{ x_k\gt -\beta_{\mathrm R}/\sigma_{\mathrm R} }[/math] 时,慢变量下降,等号成立时保持不变。这种反馈使快子系统在不同运动区间之间缓慢移动。

在原论文讨论的混沌簇放电区间,静息阶段的轨道沿稳定不动点分支移动。不动点指更新后保持原状态的点。慢变量变化使稳定与不稳定不动点相遇并消失,即发生鞍结分岔,随后轨道进入快速活动[1]。

活动阶段的轨道在混沌吸引子上运动。吸引子是邻近轨道长期演化后趋近的状态集合,吸引域是会趋向该吸引子的初始状态范围。快速活动的反馈使慢变量反向漂移,混沌吸引子与吸引域边界相遇并消失,即发生边界危机,轨道随后回到静息分支[1]。这套解释对应所研究的混沌簇放电区间,其他参数或版本的转变机制需由相应方程分析。

静息不动点与局部稳定性

快子系统分析说明簇如何开始和结束,完整二维系统的局部稳定性则说明小扰动会衰减还是增长。令式(1)中的两个变量更新后均保持不变,得到不动点 [math]\displaystyle{ (x^*,y^*) }[/math]:

[math]\displaystyle{ \begin{aligned}x^*&=-\frac{\beta_{\mathrm R}}{\sigma_{\mathrm R}},\\y^*&=x^*-\frac{\alpha_{\mathrm R}}{1+(x^*)^2}.\end{aligned} }[/math]

 

 

 

 

(4)

在不动点附近,雅可比矩阵给出小扰动的一阶更新关系。该矩阵由更新函数对状态变量的一阶偏导数组成,在式(4)处为:

[math]\displaystyle{ \begin{aligned}\mathbf J&=\begin{pmatrix}a&1\\-\sigma_{\mathrm R}&1\end{pmatrix},\\a&=-\frac{2\alpha_{\mathrm R}x^*}{[1+(x^*)^2]^2}.\end{aligned} }[/math]

 

 

 

 

(5)

对离散时间映射,矩阵特征值的模表示相应扰动每一步的放大或缩小程度。所有特征值的模小于 1 时,不动点局部渐近稳定,即充分小的扰动会逐渐消失。这里 [math]\displaystyle{ a\gt 0,\sigma_{\mathrm R}\gt 0 }[/math],线性稳定条件化为 [math]\displaystyle{ a+\sigma_{\mathrm R}\lt 1 }[/math]。

例如 [math]\displaystyle{ \beta_{\mathrm R}=\sigma_{\mathrm R}=0.001 }[/math] 时,式(4)给出 [math]\displaystyle{ x^*=-1 }[/math],式(5)给出 [math]\displaystyle{ a=\alpha_{\mathrm R}/2 }[/math],因此线性稳定范围为 [math]\displaystyle{ \alpha_{\mathrm R}\lt 1.998 }[/math]。这个界限回答不动点附近的小扰动问题;远离不动点的初值以及簇放电与持续活动的转变,还需考察整体轨道。在临界边界处,稳定性要结合非线性项判断。

对分段模型,在 [math]\displaystyle{ b_k=0,s_k=s }[/math] 且不动点位于非正分支时,慢更新给出 [math]\displaystyle{ x^*=s-1 }[/math]。Shilnikov 与 Rulkov 对完整二维系统得到以下局部分岔曲线[3]:

[math]\displaystyle{ \begin{aligned}s_{\mathrm{NS}}&=2-\sqrt{\frac{\alpha}{1-\mu}}\\&\simeq2-\sqrt{\alpha}.\end{aligned} }[/math]

 

 

 

 

(6)

式(6)在原研究的参数范围内对应亚临界 Neimark–Sacker分岔:一对共轭特征值(实部相同、虚部互为相反数)穿越单位圆,即模为 1 的边界,使离散系统的不动点发生振荡失稳。这里的“亚临界”指在稳定不动点一侧存在不稳定的不变闭曲线,参数变化使闭曲线缩入不动点时,不动点失稳[3]。右侧近似式忽略了有限 [math]\displaystyle{ \mu }[/math] 的修正,适合把握慢更新很小时的阈值位置。

分析时应区分两个对象:冻结慢变量的快子系统用于解释快速运动,完整二维系统的雅可比矩阵用于判断相互反馈后的局部稳定性。

周期、混沌与多稳态

簇放电同时包含簇内快速发放和簇间慢节律。慢节律规则,仍可伴随不规则的簇内活动,因此两种时间尺度的规则性应分别判断。混沌指确定性动力学对初值具有敏感依赖等特征的一类复杂运动,需要用相应动力学证据识别。

分段模型的快子系统可以具有稳定不动点和锋电位周期轨道,加入慢反馈后,完整系统仍可能产生混沌。2003 年的研究分析了鸭解(canard),即先沿吸引的慢结构运动、再沿排斥的慢结构停留一段时间的特殊轨道,以及相关混沌机制。该研究还展示了静息吸引子与混沌吸引子共存的参数区间[3]。这种多稳态意味着同一组参数下,不同初值可能趋向不同的长期活动。研究混沌时应结合轨道敏感性等证据,研究多稳态时应比较不同初值。

典型数值示意

图1和图2均展示式(1)的光滑映射,取 [math]\displaystyle{ \sigma_{\mathrm R}=\beta_{\mathrm R}=0.001 }[/math]。横轴是迭代步数,纵轴是无量纲模型变量。图1用于观察快变量与慢变量如何共同形成簇,图2用于比较三个参数示例的活动差异。

图1:光滑映射中快变量与慢变量的配合。取 [math]\displaystyle{ \alpha_{\mathrm R}=4.2 }[/math]。上图红线为快变量 [math]\displaystyle{ x_k }[/math],下图蓝线为慢变量 [math]\displaystyle{ y_k }[/math],快速活动与安静间隔交替出现,慢变量在两种阶段呈不同的变化方向。曲线开头含初始暂态。横轴为迭代步数,变量无量纲。图源:原词条自绘数值示意,上传记录见文件页。模型依据:[1]。
图2:光滑映射的三个参数示例。从上到下 [math]\displaystyle{ \alpha_{\mathrm R} }[/math] 分别为 1.8、4.2、4.7,展示收敛到静息不动点、簇放电和较密集的快速振荡。原图标签“紧张式发放”指第三个示例中的持续快速活动。横轴为迭代步数 [math]\displaystyle{ k }[/math],纵轴为无量纲快变量 [math]\displaystyle{ x_k }[/math]。这些参数标记所示示例,具体转变位置需结合参数扫描判断。图源:原词条自绘数值示意,上传记录见文件页。

神经科学意义与典型现象

簇同步与节律规则化

簇同步指不同单元的簇起止时刻或慢节律保持协调,锋电位同步指快速发放事件在更短的时间尺度上协调。在 2001 年的平均场网络中,耦合使簇节律变得近似周期,而簇内快速混沌振荡仍然不同步[1]。这种节律规则化主要发生在慢节律上,说明集体活动可以在保留单元快速不规则运动的同时形成稳定的时间组织。

2002 年的分段模型研究进一步比较了两个耦合单元的相对簇时序。同相簇同步指两者的簇大致同时出现,反相簇同步指两者交替出现。原研究通过改变耦合符号与参数观察到这两种情形[2]。反相结果采用模型中的负耦合参数,解释其生理对应时需结合电耦合机制与参数约定。

阈下振荡与噪声诱发放电

阈下振荡是尚未触发锋电位的小幅状态起伏。2004 年,Shilnikov 与 Rulkov 修改快函数的阈下分支,使模型先出现稳定的小幅阈下振荡,再进入锋电位发放。加入噪声后,该改进模型可以产生偶发锋电位,其锋电位间隔分布与阈下节律有关[5]。这个实例说明,调整快函数的局部形状,可以为映射模型增加特定的放电机制;相应结论以修改后的方程为依据。

放电时刻可靠性

放电时刻可靠性(spike timing reliability)衡量重复施加同类刺激时,神经元是否在相近时刻产生相似的锋电位序列。它考察跨次刺激的重复性,簇同步与锋电位同步则考察不同单元在同一次试验中的时序协调。

Bazhenov 等人于 2005 年用改进映射表示规则发放神经元、快速发放中间神经元及其耦合网络。在所考察的刺激和网络条件下,兴奋性连接有助于传播较弱输入,抑制反馈限制过多放电,两者共同改变可靠响应的频率范围[6]。研究把单元对输入节律的响应与网络相互作用联系起来,表明可靠响应的范围会随连接强度改变。

皮层网络的时空活动

映射模型也可把单元的放电事件与空间连接结合起来。2004 年的皮层网络研究用改进映射表示规则发放神经元、快速发放神经元和内在簇放电神经元,并通过兴奋性、抑制性化学突触连接这些单元。在所构造的一维和二维网络中,模型产生传播波、旋转螺旋波以及不同网络活动状态[4]。传播波表示活动依次到达不同位置,螺旋波则形成旋转的空间传播结构。

这些时空模式由单元更新与网络连接共同形成。持续波动是否出现,取决于网络大小、连接范围、边界条件和突触参数[4]。解释或复现这些结果,需要同时给出单元模型和网络设置。

耦合与网络扩展

平均场耦合

平均场耦合用群体均值作为各单元共同接收的输入。2001 年论文的全局耦合形式为[1]:

[math]\displaystyle{ \begin{aligned}x_{i,k+1}&=\frac{\alpha_i}{1+x_{i,k}^2}+y_{i,k}+\varepsilon\bar x_k,\\y_{i,k+1}&=y_{i,k}-\sigma_i x_{i,k}-\beta_i,\\\bar x_k&=\frac1N\sum_{j=1}^{N}x_{j,k}.\end{aligned} }[/math]

 

 

 

 

(7)

其中 [math]\displaystyle{ N }[/math] 是单元数,[math]\displaystyle{ i,j }[/math] 是单元编号,[math]\displaystyle{ \varepsilon }[/math] 是耦合强度,[math]\displaystyle{ \bar x_k }[/math] 是快变量的群体均值。各单元可以具有不同参数。式(7)直接加入均值,使一个单元的活动通过群体平均量影响其他单元;扩散耦合则根据单元之间的状态差产生输入。

电耦合与化学突触

对分段单元,可以用状态差构造电耦合输入。下面将原论文的两单元形式推广为多个邻接单元的加权输入[2]:

[math]\displaystyle{ \begin{aligned}C_{i,k}&=\sum_j g_{ij}(x_{j,k}-x_{i,k}),\\b_{i,k}&=c_b C_{i,k},\\s_{i,k}&=s_i+c_s C_{i,k}.\end{aligned} }[/math]

 

 

 

 

(8)

[math]\displaystyle{ g_{ij} }[/math] 是从单元 [math]\displaystyle{ j }[/math] 到 [math]\displaystyle{ i }[/math] 的耦合权重,[math]\displaystyle{ C_{i,k} }[/math] 是合成输入,[math]\displaystyle{ s_i }[/math] 是单元的基础参数。[math]\displaystyle{ c_b,c_s }[/math] 决定该输入进入快、慢更新的相对强度。若赋予无量纲输入物理电流的含义,需要给出变量尺度与电导之间的对应关系。

化学突触采用另一种输入机制:突触前神经元的锋电位触发电流,电流随后按一定时间尺度衰减,大小还与突触后神经元的状态和反转电位有关。反转电位是突触电流改变方向时对应的膜电位。2004 年皮层网络的改进模型用差分方程描述这类突触电流[4]。因此,平均场均值输入、电耦合的状态差输入和化学突触的事件触发输入,应按各自机制实现。

数值仿真与结果判读

复现鲁尔科夫映射,需要记录具体版本、完整参数、初值、输入和网络更新方式。结果判读围绕两个问题展开:更新规则是否正确实现,测量指标是否对应所关心的放电现象。

  1. 同步更新状态:先用第 [math]\displaystyle{ k }[/math] 步的全部旧状态计算新状态,再统一写入。式(1)的慢更新使用 [math]\displaystyle{ x_k }[/math]。这种方式保留原映射规定的同时更新关系,逐个原地覆盖则会增加更新顺序的影响。
  2. 区分暂态与长期行为:暂态是从初始状态过渡到长期运动的阶段。改变初值重复计算,并按需要沿参数增大、减小方向扫描,以检查多稳态及扫描方向对结果的影响。
  3. 定义放电事件:说明检测依据是阈值穿越、局部峰值还是分段映射的重置事件。用于生理比较时,还需解释所定义事件与锋电位之间的对应关系。
  4. 分别测量快、慢节律:锋电位间隔(interspike interval,ISI)是相邻锋电位的时间间隔,簇间隔则测量相邻簇之间的间隔。同步分析也应说明关注单次锋电位还是簇的起止。
  5. 核查混沌证据:最大李雅普诺夫指数描述相近轨道的平均指数分离率,返回映射比较相邻事件的状态关系。这些指标可与参数扫描共同分析混沌。分段不连续映射的指数估计还需考虑分支切换。
  6. 明确时间与电压标定:比较实测频率或膜电位前,给出迭代间隔和变量的转换关系。时间标定与所用版本及参数相联系。

优势与局限

鲁尔科夫映射的状态维度低,单步更新运算简单。时间尺度分离有助于分析快速发放与慢节律的关系,映射方法也可以用较粗的采样保留目标放电模式和较慢的网络过程,适合大量单元的模拟[4]。计算效率的比较还取决于所用版本、网络规模和所需时间分辨率。

模型参数主要服务于现象学描述,其生理解释需要额外建立映射关系。分段模型把锋电位上升与重置写入规则,便于表示放电事件;研究真实动作电位的精细波形或离子通道变化时,需要相应的生物物理模型。离散步长、重置规则和输入位置都会影响结果,加入噪声、时延或突触可塑性时,应给出这些扩展的方程、参数和依据。

与其他模型的关系

以下按时间描述方式、主要变量及解释目标比较常见神经元模型:

模型 时间与主要变量 解释侧重点及与鲁尔科夫映射的联系
Hodgkin–Huxley模型(HH) 连续时间;膜电位与表示通道开放状态的门控变量 用电导和门控机制描述动作电位,适合研究较具体的生物物理过程
FitzHugh–Nagumo模型(FHN) 连续时间;快兴奋变量与慢恢复变量 与鲁尔科夫映射共享快慢建模思路,具体放电模式与分岔结构由各自方程决定
Hindmarsh–Rose模型(HR) 连续时间;膜电位样变量及恢复、慢适应变量 用于研究簇放电,Rulkov 早期工作比较了两种模型中簇起止的机制[1]
鲁尔科夫映射 离散时间;快变量与慢变量 用低维更新规则保留选定的放电、节律与耦合现象

模型选择应围绕研究问题和所需观测量展开。通道机制与精细波形需要较具体的生物物理描述,节律、同步和较大网络的集体活动则可采用经过相应验证的映射模型。

参考文献

  1. ↑ 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 Rulkov, Nikolai F. (2001). "Regularization of Synchronized Chaotic Bursts". Physical Review Letters. 86 (1): 183–186. doi:10.1103/PhysRevLett.86.183.
  2. ↑ 2.0 2.1 2.2 2.3 2.4 Rulkov, Nikolai F. (2002). "Modeling of spiking-bursting neural behavior using two-dimensional map". Physical Review E. 65 (4): 041922. doi:10.1103/PhysRevE.65.041922.
  3. ↑ 3.0 3.1 3.2 3.3 Shilnikov, Andrey L.; Rulkov, Nikolai F. (2003). "Origin of chaos in a two-dimensional map modeling spiking-bursting neural activity" (PDF). International Journal of Bifurcation and Chaos. 13 (11): 3325–3340. doi:10.1142/S0218127403008521.
  4. ↑ 4.0 4.1 4.2 4.3 4.4 4.5 Rulkov, N. F.; Timofeev, I.; Bazhenov, M. (2004). "Oscillations in Large-Scale Cortical Networks: Map-Based Model" (PDF). Journal of Computational Neuroscience. 17 (2): 203–223. doi:10.1023/B:JCNS.0000037683.55688.7e.
  5. ↑ Shilnikov, Andrey L.; Rulkov, Nikolai F. (2004). "Subthreshold oscillations in a map-based neuron model". Physics Letters A. 328 (2–3): 177–184. doi:10.1016/j.physleta.2004.05.062.
  6. ↑ Bazhenov, Maxim; Rulkov, Nikolai F.; Fellous, Jean-Marc; Timofeev, Igor (2005). "Role of network dynamics in shaping spike timing reliability" (PDF). Physical Review E. 72 (4): 041903. doi:10.1103/PhysRevE.72.041903.

署名

贡献者:杨明哲。

原词条与配图的贡献记录见页面历史及相应文件页。