力学进展, 2020, 50(1): 202010-202010 DOI: 10.6052/1000-0992-18-021

基于大偏差理论的随机动力学研究

朱金杰,1,2,†, 陈朕2,3, 孔琛2,4, 刘先斌,2,*

1 南京理工大学机械工程学院, 南京 210094

2 南京航空航天大学机械结构力学及控制国家重点实验室, 南京210016

3 罗切斯特大学大脑与认知科学系, 美国罗切斯特 14627

4 远景能源(江苏)有限公司数字化产品开发与应用部, 上海 200051

The researches on the stochastic dynamics based on the large deviation theory

ZHU Jinjie,1,2,†, CHEN Zhen2,3, KONG Chen2,4, LIU Xianbin,2,*

1 School of Mechanical Engineering, Nanjing University of Science and Technology, Nanjing 210094, China

2 State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China

3 Department of Brain and Cognitive Sciences, University of Rochester, Rochester 14627, United States

4 Department Digital Method and Product, Envision Energy, Shanghai 200051, China

通讯作者: †E-mail:zhujinjie95@126.com;*E-mail:xbliu@nuaa.edu.cn

责任编辑: 丁千

收稿日期: 2018-09-5   接受日期: 2019-03-6   网络出版日期: 2019-03-19

Corresponding authors: †E-mail:zhujinjie95@126.com;*E-mail:xbliu@nuaa.edu.cn

Received: 2018-09-5   Accepted: 2019-03-6   Online: 2019-03-19

作者简介 About authors

朱金杰,男,1989年出生,南京理工大学机械工程学院讲师.东京工业大学2020年JSPS博士后奖学金获得者,合作教授HiroyaNakao.本科毕业于南京航空航天大学工程力学钱伟长班.研究生期间主要从事非线性随机动力学研究.参与国家自然基金项目2项.目前已在PhysicalReviewE,Chaos,ProceedingsoftheRoyalSocietyLondonA,PhysicsLettersA等期刊上发表SCI论文8篇,其中第一作者6篇.担任InternationalJournalofModernPhysicsB等期刊审稿人.主要研究兴趣包括:复杂网络、神经元动力学、振动共振、大偏差理论等.

刘先斌,男,1965年出生,南京航空航天大学航空宇航学院教授,博士生导师.主要从事非线性随机动力学:局部和全局随机分岔、噪声诱发混沌系统随机动力学行为、随机颤振、离出问题和随机共振等方面的研究.目前主持国家自然科学基金面上项目2项.已主持完成国家自然科学基金项目3项、教育部博士点基金和国家留学归国基金各1项、国家重点实验室基金项目2项.此外,作为主要研究人员,参与完成国家自然科学基金重点项目、重大研究计划和面上项目4项.已发表学术论文80余篇、论著1部.论文发表在ProceedingsoftheRoyalSocietyLondonA,ASME-JournalofAppliedMechanics,PhysicalReviewE,Chaos,NonlinearDynamics,《力学学报》《中国科学》等国际、国内学术期刊,其中SCI检索期刊论文70余篇,论文论著他引500篇次.

摘要

本文介绍了大偏差理论的基本思想、基本概念以及大偏差理论在离出问题研究中的应用.本文评述了有关离出问题的三个重要指标:平均首次离出时间、离出位置分布和最优离出路径相关研究的思路和方法,而其中对最优离出路径的刻化是结构性的难题. 针对平均首次离出时间,本文介绍了它与拟势的关系,并应用平均首次离出时间的结论分析了随机共振以及自诱导随机共振中的时间匹配机制.对于离出位置分布, 本文介绍了提高蒙特卡罗模拟速度的相关算法,并重点评述了其中的概率演化算法和相关的算例. 最后,对于最优离出路径的研究, 本文讨论了几类计算方法,分析了最优路径满足的辅助哈密尔顿系统轨线由于非线性多值性形成的拉格朗日流形拓扑结构的奇异性及其动力学含义,并进一步给出了有限噪声强度激励条件下的作用量修正方法. 最后,给出了大偏差理论应用发展的一些开放性问题的展望.

关键词： 大偏差理论 ; 离出问题 ; 最优离出路径 ; 平均首次离出时间 ; 离出位置分布

Abstract

This paper introduces the basic ideas and concepts of large deviation theory and its application in the study of exit problems. Three critical indicators of exit problems are reviewed: mean first passage time, exit location distribution and most probable escape path. Among them, the characterization of the most probable escape path is a structural conundrum. For the mean first passage time, its relationship with quasi-potential is introduced, and it is applied to analyze the time matching mechanism in stochastic resonance and self-induced stochastic resonance. For the exit location distribution, the relevant algorithms to accelerate the Monte Carlo numerical simulation are discussed, and the probability evolution method is specially clarified with some interesting examples. For the study of the most probable escape path, several calculation methods are discussed, and the singularity of the topological structure and its dynamical implications of the Lagrangian manifold formed by the auxiliary Hamilton system trajectories are analyzed. Furthermore, the corrected action method under the condition of finite noise intensity is given. Finally, the prospect of some open problems for the application and development of large deviation theory is discussed.

Keywords： large deviation theory ; exit problem ; most probable escape path ; mean first passage time ; exit location distribution

PDF (12739KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

朱金杰, 陈朕, 孔琛, 刘先斌. 基于大偏差理论的随机动力学研究. 力学进展[J], 2020, 50(1): 202010-202010 DOI:10.6052/1000-0992-18-021

ZHU Jinjie, CHEN Zhen, KONG Chen, LIU Xianbin. The researches on the stochastic dynamics based on the large deviation theory. Advances in Mechanics[J], 2020, 50(1): 202010-202010 DOI:10.6052/1000-0992-18-021

1 引言

1905年, 爱因斯坦发表了关于布朗运动的第一篇重要理论分析论文(Einstein 1905). 此后, 随机动力学逐渐发展成一门重要的学科, 在土木、机械、航空航天、海洋等工程领域, 在物理、化学、生物、生态、气象等自然科学领域, 以及在经济与金融等社会科学领域都得到了广泛的应用 (朱位秋和蔡国强 2017).

噪声通常都会被认为是有害的, 然而20世纪80年代随机共振现象的发现 (Benzi et al. 1981, Gammaitoni et al. 1998), 让研究者们看到了噪声对于动力系统建设性的积极作用. 随机共振(stochastic resonance, SR), 表现为在谐和激励和噪声共同作用的系统中, 中等强度的噪声可以提高系统响应的信噪比(signal-to-noise ratio, SNR). 而即使系统中没有谐和激励, 一定强度的噪声依然可以使得系统获得较好的响应特性, 这一现象最早被命名为自治随机共振 (Longtin 1997)或相干共振(coherence resonance, CR) (Pikovsky & Kurths 1997). 其他噪声产生的非平凡作用包括随机反共振 (Gutkin et al. 2009, Tuckwell & Jost 2012)和噪声诱导的转移 (Horsthemke & Lefever 1984)等. 另一方面, 20世纪40年代, Kolmogorov注意到了噪声对于动力系统长期行为的影响, 并认识到即使确定性系统具有渐近稳定的平衡点, 但噪声的影响可以使得系统的样本轨线具有从吸引子的吸引域离出的可能性. 换句话说, 即使是充分小的随机扰动, 在其长期作用的累积效应下, 轨线具有概率1的离出几率 (刘先斌 等 1996). 因此, 在数学文献中, 计算离出率和给定域边界上离出点概率分布的问题也被称为Kolmogorov离出问题(exit problem) (Naeh et al. 1990). 离出问题又被称为首次穿越问题, 后者更为从事动力学研究的学者所熟知.

正是噪声长期作用的特点, 确定性系统中的稳定性概念已经不再适用. 事实上, 在后面的讨论中 可以看到, 噪声的加入可以将原确定性系统的相空间"拓展一倍", 具体来讲, 如果原二维系统具有一个稳定的平衡点, 则噪声可以诱导出两个不稳定的方向, 产生两个不稳定流形; 相反地, 如果原系统具有一个不稳定的平衡点, 噪声将诱导出两个稳定的方向, 产生两个稳定流形. 因此, 随机动力系统中的稳定性概念需要重新考量. Ludwig (1975)建议在随机动力系统中用"持续性(persistence)"代替稳定性.

为了求解随机动力系统的响应, 最直接的方法就是求解与Itô随机微分方程相对应的Fokker-Planck-Kolmogorov (FPK)方程, 它描述了粒子速度的概率密度函数的时间演化信息. 然而, 通常情况下, FPK方程很难求解, 只有维数较低、满足一些特殊条件(比如, 细致平衡条件, 具有广义平稳势等)的情形才能解析求解 (Lin & Cai 1995). 为此, 一部分研究者通过采用随机平均法对FPK方程进行降维, 从而简化求解难度 (朱位秋 1987). 随机平均法优点十分显著, 因而近年来一直是随机动力学领域的热点方法(Deng & Zhu 2009, Huang & Liu 2012, Xu et al. 2013, Gu & Zhu 2014, Kong et al. 2016, Kong & Liu 2017, Pei et al. 2017, 许勇 等 2017). 经典的随机平均法分为标准和能量包络两种, 后者适用于弱白噪声和宽带实噪声激励的强非线性系统. 对于经过能量包络平均法, Hamilton函数为一维扩散过程的系统, 朱位秋 (2003)给出了一个详细的分类. 此外, 他对于随机平均方法的思想、理论和应用, 进行了详尽的评述 (朱位秋 1992, 2003). 与确定性的平均法一样, 通过对快变量在其周期内的平均, 随机平均法可以降低所研究系统的维数. 而除此之外, 随机平均法最本质、也是随机模型独有的特征, 是有关Markov极限过程的弱收敛定律—在满足强混合条件的弱噪声(量级为$\varepsilon $)激励下, 平均之后的慢变量, 如振幅、能量包络或相位(Kasminskii变换之后)在$(1/{\varepsilon ^2})$的时间尺度上弱收敛于Markov过程. 除了少数解析方法外, 对于随机动力系统, 更多地研究者通过数值计算方法获得近似的结果, 最常见的是直接通过蒙特卡罗模拟(Monte Carlo simulation)方法求解 (Higham 2001). 此外, 胞映射方法, 作为一种有效的全局动力学行为分析工具, 在非线性随机动力学中也有着不可忽视的作用 (徐伟 等 2013, 2017; 孙建桥和熊夫睿 2017). 对于离出问题, 主要的研究对象包含最有可能离出路径(most probable escape path, MPEP)、平均首次离出时间(mean first passage time, MFPT)、离出位置分布(exit location distribution, ELD)等. 其中第一项研究是结构性的, 其难度也最大. 由于噪声强度一般足够小, 因而离出问题中的平均首次离出时间都具有很高的量级, 这给直接的蒙特卡洛模拟带来了很大的麻烦. 尽管如此, 由于大自然中普遍存在的噪声作用, 化学、物理、工程、生物等各领域都会不可避免地遇到离出问题, 因而研究离出问题具有重要的现实意义. 针对离出问题中的平均首次离出时间, Kramers (1940)在20世纪40年代就给出了开创性的结果, 指出了平均首次离出时间与噪声强度之间的指数率关系. 因而研究随机动力系统中的平均首次离出时间的问题后来也被称为Kramers离出问题 (Spivak & Schuss 2002a, 2003).

直到20世纪60年代末, 有关离出问题的首个完整结果由Wentzell和Freidlin应用大偏差理论(large deviation theory, LDT)给出 (Wentzell & Freidlin 1970, Naeh et al. 1990, Freidlin & Wentzell 2012). 大偏差理论, 顾名思义, 研究随机动力系统产生的大偏差行为. 它是大数定律和中心极限定理的一个拓展和补充, 是一个处理随机过程中产生大偏差的概率的指数衰减的理论 (Touchette 2009). 与传统概率理论研究大概率事件不同, 大偏差理论关注动力系统中的小概率事件或者称为罕见事件(rare events)或极端事件(extreme events). 简单来说, 大偏差理论的核心思想为: 当小概率事件发生时, 它将以极大的概率从众多"几乎不可能(almost improbable)"的方式中, 选择"最有可能(least improbable)"的那一个. Wentzell和Freidlin (2012)通过引入作用量泛函(action functional), 刻画每一条小概率事件的样本轨迹产生的难易程度. 该作用量泛函的极值不仅刻画了小概率事件的最有可能发生方式, 也定量地给出了该事件的发生率和平均发生时间. 因而求解上述离出问题的平均首次离出时间、最有可能离出路径的问题, 就转变成了求解作用量泛函的极值问题. Freidlin-Wentzell的大偏差理论从理论上严格证明了: 对于一个即使是"结构稳定"的非线性动力系统, "弱噪声"激励使得系统长期的响应幅值越界不再是小概率事件, 此时系统产生"大偏差"的概率为1. 大偏差理论旨在从根本上揭示由噪声诱导的诸多复杂的离出现象背后的数学机理. 由于大偏差理论的广泛性以及强大的应用价值, 该理论已经在各个领域获得了大量的关注, 近期的一些应用领域包括流体力学、化学、生物学、经济学、社会学、人口学、政治学 (Wang 2015, Young 2015, Grafke et al. 2015, Agazzi et al. 2017)等.

本文主要着眼于大偏差理论的基本思想及其在离出问题上的一些近期进展, 并在最后提出了一些开放性问题.

2 大偏差理论构建

在给出大偏差理论的主体内容之前, 考虑三种稳定性的定义 (Nolting & Abbott 2016), 具体示意图见图1. 第一种稳定性为图中虚线所示, 如考虑$A$点的稳定性, 根据线性稳定性分析, 取$A$点处的特征值或者最大李雅普诺夫指数(largest Lyapunov exponent)或者Floquet指数, 它们反映了蓝色小人离开$A$点处局部的难易程度(渐近稳定性), 如果把图中的势能曲线比作一座山, 那么它们对于小人来说就是"爬得多陡". 这种稳定性是局部线性稳定性分析最常见的方法. 第二种稳定性如图中沿着水平轴从$A$点到$B$点的绿条, 对于一维系统, 取吸引域半径, 对于二维系统, 可以取吸引域面积, 以此类推. 它从几何上反映了$A$点处的平衡点的相对稳定性. 近年来Kurths团队提出的域稳定性(basin stability) (Menck et al. 2013, 2014; Schultz et al. 2017; Ji et al. 2018)的概念隶属于这一类. 对于图中小人来说, 这类稳定性反映了"走得多远". 第三类稳定性, 也是大偏差理论中采用的, 如图中双箭头所示, $\Delta V$反映了小人从$A$点到$B$点需要跨越的势能差, 化学反应中体现为活化能(activation energy), 量子隧穿和半导体中体现为位势垒(potential energy barrier). 它从能量的角度反映了小人从$A$点到$B$点所需要克服的能量大小, 这类稳定性对于小人来说是"登得多高". 三种稳定性: "爬得多陡"、"走得多远"、"登得多高"从不同角度刻画了系统的稳定性. 其中前者为局部的, 后两者为全局的. 而当 考虑系统的长期受扰行为时, "爬得多陡"显然不能准确反映系统的稳定性. 图中虽然$C$点比$A$点具有更强的局部线性稳定性, 但是由于其山谷更浅, 显然更小的能量就能驱使粒子离开$C$点越过$B$点. 而对"走得多远"来说, 尽管它从一定程度上反映了平衡点的稳定范围, 但是如果 令图中平衡点的相对势差保持不变, 将$B$和$C$之间的曲线展平, 这时候仅仅考虑平衡点的吸引域半径也是不够的. 既然外部的扰动输入可以看作对于系统的一种能量输入(类比化学反应中的活化能), 那么考虑系统的长期受扰行为最合适的还是"登得多高", 大偏差理论中这个高度差由拟势(quasipotential)来刻画 (Freidlin & Wentzell 2012).

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1   三种稳定性示意图. 其中蓝色粗实线表示势能曲线, $A$, $B$, $C$分别表示三个平衡点, $A$和$C$为稳定平衡点, 而$B$为不稳定的平衡点. $A$和$C$点处的虚线表示它们的局部曲率. 水平方向由$A$点到$B$点的绿条表示$A$点的域半径. 双箭头示意的 $\Delta V$ 表示$B$点相对于$A$点的势差


2.1 大偏差理论基本思想

传统概率论关注的都是大概率事件, 这与大数定律(law of large numbers, LLN)和中心极限定理(central limit theorem, CLT)的研究主体一致. 大数定律描述了大量重复的随机试验获得结果的平均值将会趋于数学期望, 而中心极限定理则说明具有有限均值和方差的独立同分布随机变量(对于随机变量的分布没有要求)的和依分布收敛于正态分布. 这两者都反映了大量重复的随机事件会产生大概率的近乎确定的事件. 长久以来, 人们似乎习惯了对于大概率事件的敏感性, 而容易忽视小概率事件. 举个例子, 中彩票大奖是个小概率事件, 大部分人都很难有机会赢得大奖, 它对于大部分人没有很大影响, 但对于中大奖的个体却可能改变其一生的命运. 这样看起来, 小概率事件似乎仅仅影响小部分样本, 而对于大部分样本没有影响. 但是大偏差理论却否认了这一观点. 大偏差理论关注罕见事件的发生, 如图2 所示. 尽管罕见事件的发生是小概率的, 但是如果 关注系统的长期行为, 即使是弱随机扰动, 系统产生大偏差的响应(甚至越过边界)将不再是小概率事件. 根据Freidlin和Wentzell (2012)的大偏差理论, 即使是渐近稳定的平衡点, 同时让噪声强度趋于0, 随机扰动在长时间尺度上将会使得系统状态达到任何有界的区域, 这件事情概率甚至是1的. 这是因为, 噪声强度小仅仅是指其均方值(方差)很小, 但却有小概率的情形产生大波动, 而在足够长的时间尺度上, 弱噪声的大波动将不再是小概率事件.

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2   罕见事件(rare events)示意图(朱金杰 2018)


Freidlin和Wentzell (2012)通过引入作用量泛函$S\left( \varphi \right)$给出了样本轨迹$X^\varepsilon $, 通过给定轨迹$\varphi $附近的难易程度, 其概率式具有如下的渐近形式

$ P\left\{ {\rho _{\left[ {T_1 ,T_2 } \right]} \left( {X^\varepsilon ,\varphi } \right) < \delta } \right\} \approx \exp \left\{ { - \varepsilon ^{ - 2}S\left( \varphi \right)} \right\}$

式中, $\varepsilon $表示弱高斯扰动的强度, $X^\varepsilon $为样本轨迹, $\varphi $为给定轨迹(时间区间$\left[ {T_1 ,T_2 } \right]$内的绝对连续函数, 即$\varphi \in \bar {C}\left[ {T_1 ,T_2 } \right]$), $\rho $表示$\bar {C}\left[ {T_1 ,T_2 } \right]$函数空间上的距离, $\delta $为小正数. 当噪声强度趋于0的时候, 即$\varepsilon \to {0}$, 根据式(1), 只有使得作用量泛函$S\left( \varphi \right)$取得最小值的那些轨迹具有压倒性(overwhelming)的概率, 而其他轨迹出现的概率将会被指数级地削弱.

2.2 拟势

为刻画式(1)中作用量泛函的最小值, Freidlin和Wentzell (2012)定义了拟势或准势(quasipotential)的概念. 对于时间区间$\left[ {0,T} \right]$, 拟势具有如下形式

$ W\left( {x_1 ,x_2 } \right) = \mathop {\rm inf }\limits_{T > 0} \mathop {\rm inf }\limits_{\varphi \in \bar {C}_{x_1 }^{x_2 } \left( {0,T} \right)} S_{0T} \left( \varphi \right)$

即拟势是连接起始位置$(x_{1})$与最终位置$(x_{2})$之间所有的绝对连续样本轨线的作用量的最小值, 且要求同时对路径与时间取下界. 拟势满足Hamilton-Jacobi方程, 可以通过下述Wentzel-Kramers-Brillouin (WKB)近似获得.

2.3 WKB近似

考察如下随机动力系统

$ \dot {\pmb x}_i = F_i \left(\pmb x \right)+\xi _i (t)$

式中, $\left\langle {\xi _i (t)} \right\rangle = 0$, $\left\langle {\xi _i (t)\xi _j \left( {t - \tau } \right)} \right\rangle = \varepsilon \delta _{ij} \delta \left( \tau \right)$为独立同分布高斯白噪声, $\varepsilon $为噪声强度. 与其对应的FPK方程为

$ \dfrac{\partial p\left( {\pmb x,t} \right)}{\partial t} = - \nabla \left( {F\left(\pmb x \right)p\left( {\pmb x,t} \right)} \right) + \dfrac{\varepsilon }{2}\nabla ^2p\left( {\pmb x,t} \right)$

假设系统具有唯一的稳态概率密度$p_0 \left(\pmb x \right) = \mathop {\lim }\limits_{t \to \infty } p\left( {\pmb x,t} \right)$. 当噪声强度$\varepsilon \to 0$时, 即考虑弱噪声极限的情形, 可以采用如下WKB近似(或称程函近似) (Maier & Stein 1996, Luchinsky & McClintock 1997)

$ p_0 \left(\pmb x \right)\sim z\left(\pmb x \right)\exp \left( { - \dfrac{W\left(\pmb x \right)}{\varepsilon }} \right)$

式中, $z\left(\pmb x \right)$为前因子项, $W\left(\pmb x \right)$为满足Hamilton-Jacobi方程的经典作用量, 也就是式(2)定义的拟势. 将式(5)带入式(4), 并保留至$\varepsilon $的最低阶项, 可以得到如下Hamilton-Jacobi方程

$ H\left( {\pmb x,\pmb \lambda } \right) = F\left(\pmb x \right)\lambda + \dfrac{1}{2}\pmb \lambda ^{\rm T}\pmb \lambda = {0},\qquad \pmb \lambda \equiv \dfrac{\partial W}{\partial \pmb x}$

式中$\pmb \lambda $表示动量. 注意, 由于式(3)中噪声为各向同性的(isotropic), 对于各向异性(anisotropic)噪声, 需要在式(6)中纳入扩散矩阵进行修正 (Maier & Stein 1996; Chen et al. 2016, 2017; Chen & Liu 2017b). 为求解Hamilton-Jacobi方程, 可以采用特征线方法(method of characteristics), 也叫射线方法(ray method) (Cohen & Lewis 1967, Ludwig 1975). 通过特征线方法, 可以将式(6)所示的偏微分方程转化为特征线上的常微分方程来求解. 对式(6)采用特征线法, 可以得到

$ \dot {\pmb x} = \dfrac{\partial H}{\partial \pmb \lambda } = F\left(\pmb x \right) +\pmb \lambda ,\qquad \dot {\pmb \lambda } = - \dfrac{\partial H}{\partial\pmb x} = - \left[ {\dfrac{\partial F\left(\pmb x \right)}{\partial\pmb x}} \right]^{\rm T}\pmb \lambda $

方程(7)为辅助哈密顿系统(auxiliary Hamiltonian system) (Beri et al. 2005), 它的解轨线被称为Hamilton轨迹或瞬子(instanton) (Bray & McKane 1989, Bouchet et al. 2014, Grafke et al. 2015).

对比方程(7)和原方程(3), 可以揭示出噪声扰动对于确定性系统的影响. 对于原确定性系统的平衡点$\pmb x_\ast $ (即$F\left( {\pmb x_\ast } \right) =\pmb 0)$, 相应的扩展系统(7)的平衡点为$\left( {\pmb x_\ast ,\pmb 0} \right)$ (即$F\left( {\pmb x_\ast } \right) = \pmb 0$, $\pmb \lambda = \pmb 0)$, 则很容易可以看出, 方程(7)线性化得到的特征值比原确定性系统多了一倍, 多出来的特征值正好是原系统的特征值的相反数. 因而, 对于原始确定性的系统, 噪声扮演着"搅局者"的角色, 如果原始平衡点为稳定平衡点, 噪声就给它添加不稳定的方向, 而如果原始平衡点为不稳定平衡点, 噪声就给它引出稳定的流形. 噪声似乎使得系统"鞍化"了. 所以对于随机动力系统, 稳定性的概念需要重新考虑 (Ludwig 1975).

事实上, 还可以通过拉格朗日乘子法得到方程(7). 根据Freidlin和Wentzell (2012)定义的作用量泛函, 对于弱噪声极限, 随机动力系统总是选择使得作用量泛函取得极小值的路径. 对于如下作用量泛函

$ S\left(\pmb x \right) = {\int}_{t_0 }^{t_{\rm f} } {\rm d}t L\left( {\dot {\pmb x},\pmb x} \right) = \dfrac{1}{2}{\int}_{t_0 }^{t_{\rm f} } \pmb \xi ^{\rm T}\pmb \xi {\rm d}t$
$ L\left( {\dot {\pmb x},\pmb x} \right) = \dfrac{1}{2}\left( {\dot {\pmb x} - F\left(\pmb x \right)} \right)^{\rm T}\left( {\dot {\pmb x} - F\left(\pmb x \right)} \right)$

式中, $L\left( {\dot {\pmb x},\pmb x} \right)$为经典力学中的拉格朗日函数. 它与式(6)中的哈密顿函数$H\left( {\pmb x,\pmb \lambda } \right)$满足Legendre变换 (Heymann & Vanden-Eijnden 2008). 给出如下成本函数(cost function)

$ \tilde {S}\left( {\dot {\pmb x},\pmb x,\pmb \lambda } \right) = \dfrac{1}{2}{\int}_{t_0 }^{t_{\rm f} } {\pmb \xi ^{\rm T}\pmb \xi {\rm d}t } + {\int}_{t_0 }^{t_{\rm f} } {\pmb \lambda \left[ {\dot {\pmb x} - F\left(\pmb x \right)} \right]} {\rm d}t$

求式(10)的极值问题也就是满足

$ \dfrac{\partial \tilde {S}}{\partial \dot {\pmb x}} = \dfrac{\partial \tilde {S}}{\partial \pmb x} = \dfrac{\partial \tilde {S}}{\partial \pmb \lambda } =\pmb 0$

由此, 同样可以很容易得到式(7). 相关内容可参考文献 (Bray & McKane 1989, Dykman et al. 1992, Beri et al. 2005), 在此不再赘述.

2.4 拟势的可微性

拟势是对确定性系统中的势能函数在随机动力系统中的推广, 与传统意义的势能函数相关但又有一定的区别. 以一维系统为例, 若未受扰确定性系统为梯度系统(事实上, 一维系统总是梯度系统), 即式(3)中$F\left(\pmb x \right) = - {V}'\left(\pmb x \right)$, 则Hamilton-Jacobi方程(6)化为

$ - {V}'\left(\pmb x \right){W}'\left(\pmb x \right) + \dfrac{1}{2}{W}'\left(\pmb x \right)^2 = 0$

则拟势$W\left(\pmb x \right)$可以很容易求解(Cameron 2012)

$ {W}'\left(\pmb x \right) = 0\quad \mbox{或}\quad {W}'\left(\pmb x \right) = {2}{V}'\left(\pmb x \right)$

对于任何非常数的势能函数$V\left(\pmb x \right)$, 可以求得拟势的两个解, 即

$ W\left(\pmb x \right) \equiv 0\quad \mbox{或}\quad W\left(\pmb x \right) = 2\left[ {V\left(\pmb x \right) - V\left( {\pmb x_0 } \right)} \right]$

这种情况下的拟势被称为经典解(Cameron 2012). 注意式(14)拟势与势能函数之间的两倍关系.

但是拟势并不总是处处可微的, 原系统的向量场可能不是梯度场, 这时, 称拟势为Hamilton-Jacobi方程的一个黏性解 (Cameron 2012). 有关黏性解的一些性质可以参看文献 (Crandall & Lions 1983, Crandall et al. 1984). 尽管拟势并不总是处处可微, 会在相空间的某些集合上退化至连续, 将分段可微的若干个解沿着这些集合"粘合"起来, 但是这并不影响由作用量泛函的最小值所定义的拟势的唯一性. 拟势的可微性(光滑性)对于 后面讨论的拉格朗日流形拓扑结构的奇异性至关重要.

3 平均首次离出时间

平均首次离出时间(mean first passage time, MFPT)是离出问题的一个重要研究内容和指标, 它与随机微分方程理论 (Bernt Øksendal 2010)中停时(stopping time)的概念具有密切联系, 对应于首次离出时间的随机变量在随机微分方程理论中被称为首次碰撞时间(first hitting time). 它对于衡量吸引子之间的相对稳定性具有重要的意义 (孔琛和刘先斌 2014). 从应用的角度来讲, 平均首次离出时间在广泛的领域都有它实际的应用价值. 比如, 扩散限制增长, 神经元跨阈值放电, 股票期权的触发等 (Sidney Redner 2001). 尤其是在工程领域, 平均首次离出时间对应于首穿失效(first passage failure)问题, 它与工程系统中的可靠性分析密切相关.

平均首次离出时间的最早结果可以追溯到20世纪40年代, Kramers (1940)通过对Langevin方程离出问题的研究, 给出了关于系统平均首次离出时间的一个广泛接受的近似结果

$ \bar {\tau }\sim {\exp}\left( {\dfrac{2\Delta V}{\varepsilon }} \right)$

式中, $\bar {\tau }$为平均首次离出时间, $\Delta V$表示势差, $\varepsilon $为噪声强度. 它与著名的Kramers率(Kramers rate)是互为倒数的关系 (Reimann et al. 1999). 这一结果与Freidlin和Wentzell (2012)后来建立的大偏差理论的结果一致. 他们根据流的形式的不同, 将离出问题分为3类, 分别是顺流(with a flow)、跨流(across a flow)和逆流(against a flow) (Wentzell & Freidlin 1970). 而后, Ludwig (1975)首次将数学物理方法中的射线方法或者特征线方法引入离出问题, 从而将偏微分方程转化为特征线上的常微分方程来求解. 由此, Ludwig求解了生态系统模型中的离出问题, 得到了各个种群的平均寿命. 与此同时, Schuss与他的合作者们也使用了特征线方法研究了各种化学反应模型的离出问题. 他们使用奇异摄动方法、WKB近似, 系统地分析了特征边界(characteristic boundary)与非特征边界(non-characteristic boundary)条件下的各种非线性系统的离出行为 (Matkowsky & Schuss 1977, 1982; Schuss & Matkowsky 1979; Bobrovsky & Schuss 1982; Ben-Jacob et al. 1982; Matkowsky et al. 1982, 1984; Knessl et al. 1985; Spivak & Schuss 2002b), 得到了平均首次离出时间与离出点分布的近似解析结果 (Schuss & Spivak 1998, Klosek-Dygas et al. 2006). 20世纪90年代以后, Roy (1994a, 1994b, 1995)结合奇异摄动法和射线方法计算了平均首次通过时间, 并与Monte Carlo模拟结果进行比较, 得到了很好的结果.

进入21世纪以后, 随着计算机技术的迅猛发展, 研究者们对于离出问题的研究不再受限于一般的非线性系统, 越来越多的研究工作涉及到了含有混沌现象的非线性系统的离出行为. 这些工作不但涉及具有吸引性的混沌吸引子, 还涉及到不具有吸引性的混沌鞍, 这些复杂的非线性系统在噪声的激励下会产生很多有趣的现象, 如噪声诱导的混沌(noise-induced chaos) (Tél et al. 2008, Tél & Lai 2010), 噪声诱导的间歇行为(noise-induced intermittent switching behaviors) (Kong & Liu 2017), 离出激活能量下降(lowering of the activation energy) (Kraut & Feudel 2003a, 2003b). 孔琛等通过van der Pol变换、随机平均法、 奇异摄动法以及射线方法计算了一类分段线性系统的平均首次离出时间 (孔琛 等 2014), 发现了碎片化的吸引域边界会导致理论与数值结果的偏差, 而足够光滑的吸引域边界不存在上述问题. 之后, Kong等 (2016, 2017)通过广义图胞映射方法进行全局动力学分析以及结合随机平均法降维, 研究了受高斯白噪声激励的分段线性系统的噪声诱导的混沌现象, 发现平均首次离出时间可以作为判断"噪声诱导的混沌"的一个定量指标.

这里需要特别指出的是, 以随机平均法为代表的随机动力学研究方法对于求解平均首次穿越时间具有重要贡献. 近年来, 朱位秋与其团队将随机平均法引入到离出问题的研究中. 应用随机平均法、隐式有限差分法等方法针对首穿时间和首穿失效等问题做了大量重要的工作, 得到了很多有意义的结果 (Zhu & Wu 2003; Chen et al. 2009, 2013; Chen & Zhu 2010a, 2010b, 2010c). Kougioumtzoglou和Spanos (2013)基于统计线性化和随机平均法提出了非线性振子受到宽带噪声激励情形响应和首穿统计的一种数值路径积分方法. Li等 (2015)将随机平均法推广到具有分数阶导数系统的首次穿越问题. Xu (2018)应用随机平均法研究了具有非经典非弹性碰撞的线性振子的首穿失效问题. 此外, 胞映射本身也是计算平均首次离出时间的一种有效数值工具 (Sun & Hsu 1988, Han et al. 2015).

3.1 平均首次离出时间与拟势

首次离出时间是一个随机变量, 描述为 (Freidlin & Wentzell 2012)

$ \tau ={\rm inf}\left\{ {t > 0,\pmb x_t \notin D} \right\}$

式中, $\pmb x_t $表示状态变量, $D$表示待离出区域. 平均首次离出时间是其均值, 满足如下Pontryagin方程 (Matkowsky et al. 1982, Naeh et al. 1990, Zhu et al. 2002) (狄利克雷边值问题)

$ {\rm L}^\ast \bar {\tau } = - 1,\quad \pmb x \in D; \qquad \bar {\tau } = 0,\quad\pmb x \in \partial D$

式中${\rm L}^\ast $为Kolmogorov后向算子$({\rm L}^\ast = a_{ij} \left( {\pmb x_0 ,t_0 } \right)\dfrac{\partial }{\partial \pmb x_i } + \dfrac{1}{2}b_{ij} \left( {x_0 ,t_0 } \right)\dfrac{\partial ^2}{\partial x_i \partial x_j }$, $a_{ij}$和$b_{ij}$分别为漂移项和扩散项), 它是FPK算子的伴随算子, $\partial D$为区域$D$的边界. 求解方程(17)的方法有很多, 包括奇异摄动法、有限差分法、路径积分法等. 这些分析方法都不可避免地讨论区域$D$内的稳态概率密度$p_0 \left(\pmb x \right)$, 而它与拟势具有下面的关系 (Freidlin & Wentzell 2012)

$ \mathop {\lim }\limits_{\varepsilon \to 0} - \varepsilon \ln p_0 \left(\pmb x \right) = W\left(\pmb x \right)$

式中, $\varepsilon $为噪声强度, $p_0 \left(\pmb x \right)$为稳态概率密度, $W\left(\pmb x \right)$为拟势. 或者更加直接的有 (Ludwig 1975)

$ \ln \bar {\tau }\sim \dfrac{1}{\varepsilon }W\left(\pmb x \right)$

这一结果与式(15)是一致的. 只是对于梯度系统, $W\left(\pmb x \right) = 2\Delta V\left(\pmb x \right)$; 而对于非梯度系统, 式(15)中的$2\Delta V$必须用拟势来代替. 从式(19) 可以看出, 对于固定的离出点(即固定拟势), 平均首次离出时间是噪声强度倒数的指数函数; 而对于固定的弱噪声强度, 拟势的大小对于平均首次离出时间的影响也是指数的.

3.2 时间尺度匹配

平均首次离出时间对于许多随机动力学现象具有重要意义, 当它与系统内禀的或者外在的时间尺度相匹配的时候, 系统会产生类似于共振的行为. 这边 分别以随机共振和自诱导随机共振(self-induced stochastic resonance, SISR)为例.

3.2.1 随机共振

随机共振现象在20世纪80年代初被发现 (Benzi et al. 1981), 并首先被应用于大气科学 (Benzi et al. 1982). 对于受到随机扰动以及简谐激励的双稳态系统, 随机共振表现为中等强度的噪声可以增强系统响应, 在功率谱图上表现为显著的峰值. 近年来随机共振在多个领域表现出强大的应用前景. 机械系统中可以通过应用随机共振来增大弱故障特征, 从而作为一种故障诊断的手段 (Lu et al. 2017, Qiao et al. 2017, Huang et al. 2018). 在医学领域, 随机共振也在许多方面展现了它重要的价值, 比如随机共振可以促进前庭脊髓反射 (Wuehr et al. 2017), 它在人工耳蜗方面也具有重要应用 (Chatterjee & Robert 2001, Stocks et al. 2002), 随机共振甚至是癌细胞对抗治疗背后一个重要机制 (Rodrigo & Stocks 2018). 总之, 随机共振的反直觉特征以及其广泛的应用前景, 一直是各个领域的研究者们的研究热点.

随机共振背后的机制与平均首次离出时间具有密切联系 (Gammaitoni et al. 1998). 以具有对称势阱的双稳态系统随机共振为例, 由于周期激励的作用, 系统的势阱深度会被周期性地调制, 如图3 所示. 图3 中小球从势阱跨越势垒的时间尺度由式(15)的平均首次离出时间给出, 这一时间尺度与势阱深度$\Delta V$有关, 显然图中$t = \dfrac{\pi }{2\varOmega }$时, 左侧势阱深度更浅, 因而平均首次离出时间更小, 更易从势阱中逃逸, 从而转移到另一侧. 同样的, 对于$t = \dfrac{ {3}\pi }{2\varOmega }$的情况, 右侧势阱深度更浅, 因而更易从势阱中逃逸进入左侧. 理论上, 随机共振发生时, 有如下时间匹配条件 (Gammaitoni et al. 1998)

$ \bar {\tau } = \dfrac{1}{2}T_\varOmega$

式中, $\bar {\tau }$表示粒子从一侧势阱离出的平均首次离出时间, $T_\varOmega $为外激励的周期. 也就是说, 当粒子的平均首次离出时间为外激励周期的一半时, 随机共振产生, 期间粒子经历了从一侧转移至另一侧然后再转移回来的两次离出过程. 因而对于强度小的噪声, 即使是最浅时候的势阱也无法使得平均首次离出时间满足式(20). 而对于强度太大的噪声, 即使是最深时候的势阱也因为平均首次离出时间太小而无法与外激励周期产生时间尺度的匹配. 只有合适强度的噪声使得式(20)满足, 随机共振现象才能发生.

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3   对称双势阱随机共振示意图. 其中双势阱势能函数为$V\left(\pmb x \right) = x^4/4 - x^2/2$, 由于周期激励的作用, 势能函数将会变为$V\left(\pmb x \right) = V\left(\pmb x \right) - Ax\sin \left( {\varOmega t} \right)$, $A$为周期激励幅值, $\varOmega $为其频率, $\Delta V$为势阱深度


3.2.2 自诱导随机共振

当随机动力系统中没有周期激励的作用, 噪声依然可以让系统产生相干的行为, 这一现象被命名为相干共振 (Pikovsky & Kurths 1997). 而后, Vanden-Eijnden团队 (Lee DeVille et al. 2005, Muratov et al. 2005)发现噪声作用在不同时间尺度上产生的相干行为的机制不同, 在慢时间尺度上的相干即上述相干共振, 而在快时间尺度上的噪声产生的相干行为被称为自诱导随机共振. 前面讨论的随机共振是平均首次离出时间与外激励周期之间的时间尺度匹配, 而自诱导随机共振是平均首次离出时间与系统内禀的时间尺度之间的匹配. 以Brusselator系统为例, 系统方程如下 (Muratov et al. 2005)

$ \dot {u} = 1 + u^2v - \left( {1 + A} \right)u + \sqrt \varepsilon \xi ,\qquad \dot {v} = \alpha \left( {Au - u^2v} \right)$

式中, $A$为控制参数, 取$A=0.7$, $\alpha = 0.01$保证$u$和$v$变量之间的时间尺度分离. $\xi (t)$为高斯白噪声, $\varepsilon $为其噪声强度, 固定其强度为0.1. 对于确定性的Brusselator系统, 图4(a)为其确定性系统流线图, 注意当变量$v$的值取在0.747左右时, 相轨线对于初值的取值十分敏感, 有的轨线会迅速返回平衡点$(1, A)$, 而有的轨线会经历大的振荡, 这样的系统被称为可激系统(excitable system), 具有图中所示的叶状(foliation)结构. 对于上述参数, 数值模拟结果展现了惊人的相干性, 如图4(b)插图所示, 时间序列具有近乎确定的周期性(相图中随机轨线的剧烈毛刺是由于横坐标的对数形式导致的). 这样的相干性是由平均首次离出时间与系统内禀的时间尺度相匹配导致的, 这点与随机共振具有类似的地方, 下面来具体分析.

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4   自诱导随机共振现象 (Muratov et al. 2005). (a) Brusselator系统相图, 其中红色与蓝色粗实线分别为$u$零倾线及$v$零倾线, 黑色带箭头细实线为确定性系统的相轨线;(b) Brusselator系统自诱导随机共振现象, 灰色细实线为噪声激励下的样本轨线, 黑色粗实线与虚线为理论上的随机极限环, 插图为样本轨线时间序列. 注意横轴的对数尺度, 数值模拟采用Euler算法


由于系统(21)为快慢动力系统, 当轨线运行到$u$零倾线的左侧分支$u_ - \left( v \right)$时, 可以冻结慢变量$v$, 而研究快子系统的动力学行为. 由系统(21)的第一个式子, 可以求得对应于快子系统的势能函数 (Muratov et al. 2005)

$ V\left( {u,v} \right) = - u - \dfrac{1}{3}u^3v + \dfrac{1}{2}\left( {1 + A} \right)u^2$

注意, 式(22)中的$v$沿着$u_ - \left( v \right)$可以作为参数, 当$v$值固定时, 式(22)中的势能函数为$u$的函数, 见图5. 当样本轨线沿着$u_ - \left( v \right)$运动时, 对于每一个$v$值, 随机扰动都有可能使得轨线越过势垒而发生一次大幅值振荡. 因而, 类比随机共振中的讨论, 自诱导随机共振发生要求系统平均首次离出时间与沿着慢流形$u_ - \left( v \right)$运动的弛豫时间必须相匹配. 具体表达式为

$ \alpha ^{ - 1}\sim \exp \left( {\dfrac{2\Delta V}{\varepsilon }} \right)$

或者

$ \dfrac{1}{2}\varepsilon \ln \alpha ^{ - 1}\sim \Delta V$

式中, $\Delta V$为$u$零倾线右分支与左分支的势能差

$ \Delta V = V\left( {u_ + \left( v \right),v} \right) - V\left( {u_ - \left( v \right),v} \right) = \dfrac{\left( {\left( {1 + A} \right)^2 - 4v} \right)^{3 / 2}}{6v^2}$

由式(24)可以求得的临界的离出位置$v_\ast $, 对应于图4(b)中的黑色粗实线. 当样本轨线没有运动到$v_\ast $位置时, 平均首次离出时间远大于沿着慢流形运动的弛豫时间; 而当样本轨线越过$v_\ast $时, 平均首次离出时间远小于样本轨线沿着慢流形运动的弛豫时间, 离出大概率发生. 需要注意的是, 式(23)或式(24)理论上预测的临界离出点要求噪声强度$\varepsilon $和时间尺度分离参数$\alpha $都趋于0, 因而对于有限值的情形会存在一定的误差.

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5   Brusselator系统自诱导随机共振时间尺度匹配示意图. 其中, 红色和蓝色粗实线分别为$u$零倾线及$v$零倾线, 黑色细实线为不同$v$值对应的势能曲线, $T_{\rm s}$和$T_{\rm e}$分别为慢子系统时间尺度和不同$v$值对应的平均首次离出时间


相比于相干共振, 自诱导随机共振可以在很大的参数范围内产生相干行为(Muratov et al. 2005). 这是因为相干共振要求系统参数接近分岔值, 而自诱导随机共振没有这个限制, 因而具有更强的参数鲁棒性, 可以预见其拥有更广泛的应用前景.

4 最优离出路径

平均首次离出时间解决了噪声导致的离出"何时发生"的问题, 对于动力系统而言, 离出"怎样发生"也很重要, 它反映了相轨线在相空间中的运动轨迹, 特别地, 它能给出离出过程噪声扰动的大小, 这对于认识噪声在离出过程中的作用, 以及相应的最优控制问题至关重要. 离出问题中描述离出"怎样发生"的是最优离出路径(或者最有可能离出路径), 在大偏差理论中对应着最小的作用量泛函轨迹, 由于其满足Hamilton-Jacobi方程, 因而也被称为零能轨迹(zero-energy trajectories) (Maier & Stein 1993a, Luchinsky et al. 1997, Smelyanskiy et al. 1997).

事实上, 对于满足细致平衡条件的平衡态系统(equilibrium system), 涨落耗散定理(fluctuation-dissipation theorem) (Weber 1956, Kubo 1966)早就给出了最优离出路径的答案: 涨落轨迹是耗散(弛豫)轨迹的逆过程. 以如下受白噪声扰动的一维有势系统为例 (Dykman et al. 1992)

$ \dot {x} = - {V}'\left( x \right) + \xi (t)$

式中$V(t)$为系统势能函数, $\xi (t)$为高斯白噪声, 通过式(12)和式(13), 并结合辅助Hamilton系统(7) 可以很容易得到最优离出路径满足的方程 (Dykman et al. 1992)

$ \dot {x} = - {V}'\left( x \right) + 2{V}'\left( x \right) = {V}'\left( x \right)$

由式(27)很容易可以看出, 最优路径的系统向量场与原系统的向量场大小相等, 方向相反, 因而涨落是耗散的逆过程. 如果对于从稳定平衡点离出至鞍点的过程, 最优离出路径就是从鞍点弛豫回平衡点的相反路径. 尽管涨落耗散定理给出了一些系统涨落与耗散的深刻关系, 但是对于不满足细致平衡条件的情形(或者说不具有时间可逆性的系统), 这一定理还无法给出答案. 大偏差理论并不需要系统满足细致平衡条件, 它从能量的角度, 以作用量泛函的极值给出了最优路径的结果.

有趣的是, 尽管大偏差理论计算得出的最优离出路径针对的是噪声扰动的随机动力系统, 如果考虑确定性系统的最优控制, 最优离出路径所对应的动量作为一种控制力将会成为最优控制力(optimal force) (Smelyanskiy & Dykman 1997, Luchinsky et al. 2002, Khovanov et al. 2003, Chen & Liu 2017b). 此外, 由作用量泛函的最小值所确定的最优离出路径, 不仅是理论意义上的结果, 当噪声强度足够小时, 在实验中也是可以被观测的 (Dykman et al. 1996; Luchinsky et al. 1997, 1998; Chan et al. 2008).

4.1 最优离出路径的计算

理论上, 最优离出路径需要求解相应的Hamilton-Jacobi方程, 它是一个偏微分方程, 一般情形下很难解析求解, 在这里 不做过多讨论. 下面 给出一些典型的数值算法, 作为求解最优离出路径的强有力的工具.

4.1.1 历程概率密度

求解最优离出路径, 最直接的方法是采用蒙特卡罗模拟(Monte Carlo simulation). 但是由于首次离出时间是个随机变量, 无法确定某条样本轨迹的离出时间(只能求得它的平均值), 因而无法直接得到准确的离出样本. 科学家们总能找到办法, Dykman 等 (1992)通过反向思考引入了历程概率密度(prehistory probability density)方法解决了这一问题.

历程概率密度$p_{\rm h} \left( {\pmb x,t;\pmb x_{\rm f} , t_{\rm f} } \right)$描述了在$ t_{\rm f} $时刻到达$\pmb x_{\rm f} $的条件下, 在$t < t_{\rm f} $时刻系统处于$x$处的概率. 事实上, 它并不是标准意义下的两时刻转移概率, 它是下面三时刻转移概率与两时刻转移概率的比值 (Dykman et al. 1992)

$ p_{\rm h} \left( {\pmb x_{\rm i} , t_{\rm i} ; \pmb x,t;\pmb x_{\rm f} , t_{\rm f} } \right) = \dfrac{w\left( {\pmb x_{\rm f} , t_{\rm f} ; \pmb x,t;\pmb x_{\rm i} , t_{\rm i} } \right)}{w\left( {\pmb x_{\rm f} , t_{\rm f} ;\pmb x_{\rm i} , t_{\rm i} } \right)}$

式中, $w\left( {\pmb x_{\rm f} , t_{\rm f} ; \pmb x,t;\pmb x_{\rm i} , t_{\rm i} } \right)$表示状态从$\pmb x_{\rm i} \to\pmb x \to \pmb x_{\rm f} $的转移概率, $w\left( {\pmb x_{\rm f} , t_{\rm f} ;\pmb x_{\rm i} , t_{\rm i} } \right)$同理. 由于当终点$\pmb x_{\rm f} $远离吸引子时, 连续通过终点$\pmb x_{\rm f} $的时间间隔将会相当大(噪声强度足够小时到达的时刻近似符合泊松统计), 这一时间间隔远大于系统的弛豫时间以及噪声的相关时间(如果是色噪声的话), 则 可以令初始时刻$ t_{\rm i} \to - \infty $, 连同初始位置$\pmb x_{\rm i} $可以一同舍弃, 最终得到历程概率密度$p_{\rm h} \left( {\pmb x,t;\pmb x_{\rm f} , t_{\rm f} } \right)$. 显然, 可以对历程概率密度做时间平移: $p_{\rm h} \left( {\pmb x,t;\pmb x_{\rm f} , t_{\rm f} } \right)=p_{\rm h} \left( {\pmb x,t - t_{\rm f} ;\pmb x_{\rm f} ,0} \right)$.

历程概率密度方法具体操作如下: 首先 进行蒙特卡罗模拟, 当样本到达边界时, 模拟停止, 并记录下从离出时刻向前递推固定时间长度的样本路径(一般真正的离出过程时间很短, 模拟在大部分时间都在吸引子附近徘徊). 得到大量的这样的离出样本后统计不同时刻的概率密度, 最终概率密度的山脊(峰值)给出了最优离出路径.

图6给出了系统(26)中势能函数$V\left(x \right) = - x^2/2 + x^4/4$的情形的历程概率密度图, 离出点定为$x_{\rm f} = - 0.3$, 图中红色粗实线为其山脊(峰值), 对应着最优离出路径. 事实上, 历程概率密度方法除了可以从数值上求解最优离出路径, 还可以求解离出样本的弥散(dispersion)程度, 从图6 中还可以看出, 离出路径经历了弥散程度增大再缩小的过程(图中曲线变缓的情形), 存在一个最大的弥散位置.

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6   历程概率密度方法求解最优离出路径 (Dykman et al. 1992). 其中, 黑色细实线为不同时刻离出路径的概率密度统计, 红色粗实线为概率密度的峰值, 它给出了最优离出路径4.1.2 Action plot方法求解最优离出路径最准确的方法是求解Hamilton-Jacobi方程(6), 通过特征线方法, 可以转化为求解辅助Hamilton系统(7). Beri等 (2005)引入了一种拓扑方法, 名为Action plot方法, 通过枚举所有满足一定初始条件或边界条件的离出轨迹, 从中选出具有最小作用量的那一条轨迹, 即为最优离出路径. 更准确地说, 该方法可以计算出辅助Hamilton系统所有从吸引子出发、并张成整个拉格朗日流形(Lagrangian manifold)的解轨线在离出吸引域时所取得的作用量, 而这些作用量被参数化为辅助Hamilton系统的解轨线的初始条件的函数.


下面 以inverted van der Pol (IVDP)振子为例, 详述其求解过程.其动力学模型如下 (Beri et al. 2005)

$ \ddot {x} = - 2\eta \left( {1 - x^2} \right)\dot {x} - x + \sqrt {4\eta T} \xi (t)$

或者以状态变量的形式改写为

$ \dot {x} = y,\qquad \dot {y} = - 2\eta \left( {1 - x^2} \right)y - x + \sqrt {4\eta T} \xi (t)$

式中, $\eta = 0.5$为参数, $\xi (t)$为高斯白噪声, 满足$\left\langle {\xi (t)\xi \left( {t - \tau } \right)} \right\rangle = \delta \left( \tau \right)$, $\sqrt {4\eta T} $为其强度. IVDP振子具有一个稳定的平衡点和一个不稳定的极限环. 相应的Hamilton-Jacobi方程为(注意这边为各向异性噪声)

$ H\left( {\pmb x,\pmb \lambda } \right)=y\lambda _1 + \left[ { - 2\eta \left( {1 - x^2} \right)y - x} \right]\lambda _2 + \dfrac{1}{2}\lambda _2^2$

通过特征线方法得到的辅助Hamilton方程为

$ \begin{array}{l} \dot {x} = y,\qquad \dot {y} = - 2\eta \left( {1 - x^2} \right)y - x + \lambda _2 \\ \dot {\lambda }_1 = \left( {1 - 4\eta xy} \right)\lambda _2 ,\qquad \dot {\lambda }_2 = - \lambda _1 + 2\eta \left( {1 - x^2} \right)\lambda _2 \\ \end{array}$

接下来, 需要确定辅助Hamilton系统的初值选取. 首先, 根据文献 (Beri et al. 2005)的方法, 选取以平衡点为中心半径极小$( \ll \sqrt {4\eta T} )$的一个单位圆, 初始坐标取值为$\left( {\delta \cos \left( \varPhi \right) , \delta \sin \left( \varPhi \right)} \right)$, $\delta $为小量, $\varPhi $取$\left[ {0,2\pi } \right]$之间足够多的点, 然后根据拟势在初始点处的高斯近似 (Roy 1993), 或者采用平衡点处附近的线性化近似 (Beri et al. 2005), 可以选取动量的初值. 图7(a)给出了这些离出轨迹的作用量结果, 其中$a$和$b$处为两个作用量最小值点(两个最小值点是由于系统的中心对称导致的), $d$处为作用量的局部极小值点, $c$处为局部极大值点. 黑色矩形框的局部放大细节展示了作用量图的自相似性. 图7(b)给出了与图7(a)中相对应的离出轨迹, 从图上可以看出, $a$和$b$曲线作为中心对称的最优离出路径是连接平衡点与极限环的异宿轨迹(heteroclinic trajectories).

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7   Action plot方法求解IVDP系统最优离出路径 (Beri et al. 2005). (a)辅助Hamilton系统参数化轨线到达不稳定极限环处的作用量图, 其中$a$和$b$对应两个作用量全局最小值点, $d$对应局部作用量极小值点, $c$对应局部作用量极大值点, 黑色矩形框局部放大图反映了作用量图的自相似性; (b)相应的离出轨迹


4.1.3 几何最小作用量方法

从前面的讨论可以看到, 作用量是轨线的泛函, 如果 对于给定轨迹重新参数化, 那么式(1)或者式(8)中定义的Freidlin-Wentzell作用量将变成几何作用量 (Heymann & Vanden-Eijnden 2008, Freidlin & Wentzell 2012, Cameron 2012)

$ S_g \left( \psi \right) = {\int}_0^{\alpha _{\max } } {\left\{ {\left| {F\left( {\psi \left( \alpha \right)} \right)} \right|\left| {\psi _\alpha \left( \alpha \right)} \right| - F\left( {\psi \left( \alpha \right)} \right) \cdot \psi _\alpha } \right\}{\rm d}\alpha }$

该几何作用量与前面定义的Freidlin-Wentzell作用量由下面不等式给出

$ S\left( \varphi \right) = \dfrac{1}{2}{\int}_0^T {\left| {\dot {\varphi } - F\left( {\varphi (t)} \right)} \right|^2} {\rm d}t = \dfrac{1}{2}{\int}_0^T {\left\{ {\left| \dot {\varphi } \right|^2 + \left| {F\left( {\varphi (t)} \right)} \right|^2 - 2\dot {\varphi } \cdot F\left( {\varphi (t)} \right)} \right\}} {\rm d}t \\ \geq {\int}_0^T {\left\{ {\left| \dot {\varphi } \right|\left| {F\left( {\varphi (t)} \right)} \right| - \dot {\varphi } \cdot F\left( {\varphi (t)} \right)} \right\}} {\rm d}t \\ = {\int}_0^{\alpha _{\max } } {\left\{ {\left| {F\left( {\varphi \left( \alpha \right)} \right)} \right|\left| {\varphi _\alpha \left( \alpha \right)} \right| - F\left( {\varphi \left( \alpha \right)} \right) \cdot \varphi _\alpha } \right\}{\rm d}\alpha } = S_g \left( \varphi \right)$

式中, 不等式成立的条件为$\left| \dot {\varphi } \right| = \left| {F\left( \varphi \right)} \right| $, 而这一等式对于最优离出路径总是成立的 (Heymann & Vanden-Eijnden 2008). 在这样的情况下, 几何最小作用量与Freidlin-Wentzell作用量将会等价. 运用这样的思想, Heymann & Vanden-Eijnden (2008)构建了几何最小作用量方法(geometric minimum action method, gMAM)来求解最优离出路径, 方法的核心是预定义一条空间曲线, 然后按照一定的收敛规则(差分方程)让预定义的空间曲线收敛至最优路径, 方法的详细建立过程参见文献 (Heymann & Vanden-Eijnden 2008).

这里 以Maier-Stein系统为例, 运用gMAM来求解两种参数条件下的最优离出路径, 系统方程如下 (Maier & Stein 1993a, 1993b, 1996)

$ \dot {x} = x - x^3 - \alpha xy^2 + \xi _1 (t),\qquad \dot {y} = - \mu y\left( {1 + x^2} \right) + \xi _2 (t)$

其中$\left\langle {\xi _i (t)} \right\rangle = 0$, $\left\langle {\xi _i (t)\xi _j \left( {t - \tau } \right)} \right\rangle = D\delta _{ij} \delta \left( \tau \right)$. 分别选取$\alpha = 6.67$, $\mu = 2$以及$\alpha = 10$, $\mu = 0.67$这两种参数条件, 所得的最优离出路径分别如图8(a)和图8(b) 所示. 图中也给出了初始曲线在演化过程中的收敛情形.

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8   几何最小作用量方法求解Maier-Stein系统最优离出路径 (Heymann & Vanden-Eijnden 2008). (a) $\alpha = 6.67$, $\mu = 2$; (b) $\alpha = 10$, $\mu = 0.67$. 图中灰色细线为确定性系统的流线, 蓝色粗实线为演化的初始曲线, 红色粗虚线为最终的最优离出路径, 中间的彩色细线展示了在演化过程中初始曲线的收敛情况


需要注意的是, 通过gMAM得到的最优离出路径可能收敛于局部最优解, 因而在实际使用过程中需要尝试不同初始曲线并与其他方法的结果进行对比以保证正确性.

4.1.4 有序逆风算法

在介绍这一算法之前, 要强调一下拟势的重要性. 一旦随机动力系统的拟势确定了, 它在相空间的能量结构也就确定了. 式(19)表明, 通过拟势 可以计算系统的平均首次离出时间(当然, 理论上还需要计算指数前因子项). 而有序逆风算法有助于 数值地计算拟势, 进而通过拟势确定最优离出路径.

有序逆风算法(ordered upwind method, OUM)最早是由Sethian & Vladimirsky (2001, 2003)建立的, Cameron (2012)通过改进它来计算随机动力系统的拟势, 具体算法详见文献. 这里 假设已经通过OUM算法计算得到了系统的拟势$\nabla U$, 那么假设原始向量场为$\pmb F\left(\pmb x \right)$, 那么最优离出路径满足的向量场为 (Cameron 2012, Chen et al. 2017)

$ \pmb G = \dfrac{\pmb F + \nabla U}{\left| \pmb F \right|}$

根据式(36)可以从离出点反向射出流线来得出最优离出路径.

在给出算例之前, 为了更好地理解拟势与原确定性向量场的关系, 可以给出向量场分解图, 如图9 所示. 从图上可以看出, 当系统向量场的旋转分量为零时, 系统为有势系统, 向量场满足关系: $\pmb F(\pmb x ) = - \dfrac{1}{2}\nabla U = -\pmb G_F \left(\pmb x \right)$, 此时最优离出路径为确定性轨迹的时间逆过程, 即涨落耗散定理所满足的关系.

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9   向量场的分解. 图中$\pmb F(\pmb x)$ 为确定性系统的向量场, $\nabla U$为拟势的梯度, $\pmb G_F \left(\pmb x \right) = \nabla U + \pmb F$为最优离出路径满足的向量场, $\pmb l(\pmb x)$ 表示旋转分量


下面 仍然以4.1.3节中的Maier-Stein系统(35)为例, 取参数$\alpha = 10$, $\mu = 1$, 应用OUM计算系统的拟势以及最优离出路径, 结果如图10 所示. 其中图10(a)给出了拟势的三维图及等高线图, 图10(b)中红线表示从$(-1,0)$点离出(转移)至(1,0)点的最优离出路径.

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10   有序逆风算法求解Maier-Stein系统最优离出路径. (a)拟势三维图及等高线图, (b)最优离出路径(红线)和向量场$\pmb G$ (黑色箭头)


Chen等 (2017)将OUM方法运用到了FitzHugh-Nagumo神经元系统, 并将其结果与gMAM以及action plot方法相比较, 取得了一致的结果. 最近, Cameron团队将OUM算法推广到了各向异性噪声 (Dahiya & Cameron 2018)和高维情形 (Yang et al. 2019), 因而这一方法有望获得较广泛的应用.

4.1.5 最优离出路径数值算法比较

上述4种计算最优离出路径的数值算法各有其优缺点. 其中, 历程概率密度方法有别于其余3种方法, 前者针对有限强度的噪声, 得出具体噪声强度下的最优离出路径, 而后3种方法得出的最优离出路径是针对噪声强度趋于零时的极限情形. 当噪声强度趋于零时, 系统的平均首次离出时间趋于无穷, 采用Monte Carlo模拟将变得不切实际, 造成历程概率密度方法无法实施, 而此时后3种方法可以得出较为精确的结果. gMAM方法由于采用差分方程, 在这4种方法中计算效率最高, 但是采用这一方法需要预先给定初始点和离出点, 这会给初始或者离出边界为其他复杂不变集的情形(如混沌吸引子)带来麻烦; 此外, 若初始曲线选定不当, 也会造成收敛到局部最优解的问题. Action plot方法可以通过枚举所有Hamilton轨迹的方式计算出每一条轨迹的作用量, 从而选出其中作用量值最小的那一条作为最优离出路径, 但它对于初始或者离出边界为复杂不变集的情形同样存在初值选取的问题. OUM方法直接得出系统的全局拟势结构, 因而计算结果更为精确, 得到的最优离出路径一定是全局最优解.

综上所述, 在实际计算最优离出路径过程中, 要配合使用这4种数值计算方法, 以保证结果的可靠性.

4.2 拉格朗日流形拓扑结构的奇异性

在前面的讨论中,没有考虑离出路径或者涨落路径的唯一性. 事实上, 涨落至一点或者离出的路径并不总是唯一的, 这体现在作用量泛函的多值性, 或者从流形的拓扑结构上来讲, 体现在拉格朗日流形(Lagrangian manifold, UM)的奇异性上. 本节 讨论拉格朗日流形拓扑结构的奇异性, 并给出几个实际的例子.

4.2.1 切换线与焦散线

拉格朗日流形是由辅助Hamilton系统(7)的解轨迹张成的不变流形, 或者说是从吸引子处扩展出的不稳定流形. 它反映了受扰系统涨落或者离出轨迹的拓扑结构. 前面2.4节讨论过拟势的可微性, 根据式(2), 拟势为作用量泛函在给定边界条件的所有绝对连续曲线取最小值得到的, 尽管最小值是唯一的, 但是对应的作用量泛函最小值曲面未必是一张光滑曲面, 很多情形下该曲面为多张分段光滑的曲面横截相交形成的, 沿着这些交线, 拟势在其横截方向上一阶偏导数不连续, 也就是说, 在这些点处拟势不可微, 这些点集构成的曲线被称为切换线(switching line, shock line). 切换线具有十分明确的物理意义. 对于切换线两侧距离十分相近的两个终点, 会有两条拓扑性质完全不同的最优路径到达; 当终点落在切换线上时, 会有两条作用量完全相同的不同最优路径抵达. 它的出现反映了系统细致平衡条件的缺失、对称性的破缺, 甚至对于后面讨论的离出位置分布也会有着本质的影响 (Maier & Stein 1993a).

2.3节中将WKB近似代入FPK方程保留至噪声强度的最低阶 得到了拟势所满足的Hamilton-Jacobi方程(6), 而当保留至噪声强度的次低阶时, 可以得到前因子项满足的输运方程(transport equation) (Roy 1994a)

$ \left( {F_i \left(\pmb x \right) + \dfrac{\partial W}{\partial \pmb x_i }} \right)\dfrac{\partial z}{\partial \pmb x_i } + \left( {\dfrac{\partial F_i }{\partial \pmb x_i } + \dfrac{1}{2}\dfrac{\partial ^2W}{\partial \pmb x_i \partial \pmb x_j }} \right)z\left(\pmb x \right) = 0$

根据辅助Hamilton系统(7), 很容易将式(37)改写成

$ \dot {z} = - \left( {\dfrac{\partial F_i }{\partial \pmb x_i } + \dfrac{1}{2}\dfrac{\partial ^2W}{\partial \pmb x_i \partial \pmb x_j }} \right)z$

式中, 拟势的二阶导数$\dfrac{\partial ^2W}{\partial \pmb x_i \partial \pmb x_j }$由如下Riccati方程得到 (Maier & Stein 1993a, 1996)

$ \dfrac{\rm d}{{\rm d}t}\left( {\dfrac{\partial ^2W}{\partial \pmb x_i \partial \pmb x_j }} \right) = - \dfrac{\partial ^2W}{\partial \pmb x_i \partial \pmb x_k }\dfrac{\partial ^2W}{\partial \pmb x_k \partial \pmb x_j } - \dfrac{\partial F_k }{\partial \pmb x_i }\dfrac{\partial ^2W}{\partial \pmb x_k \partial \pmb x_j } - \dfrac{\partial F_k }{\partial \pmb x_j }\dfrac{\partial ^2W}{\partial \pmb x_k \partial \pmb x_i } - \dfrac{\partial ^2F_k }{\partial \pmb x_i \partial \pmb x_j }\dfrac{\partial W}{\partial \pmb x_k }$

通过联立辅助Hamilton系统(7)、输运方程(38)以及Riccati方程(39), 可以计算出前因子项随Hamilton轨迹的变化情况. 实际运算过程中, 拟势的二阶导数项会在相空间的某些位置发散, 指数前因子项也会因此发散 (Kong & Liu 2017, 孔琛 2018). 这一现象与辅助Hamilton系统在其构形空间(configuration space)的拓扑结构有关. 以二维系统为例, 则辅助Hamilton系统为四维, 相应的拉格朗日流形为二维. 根据定义, 拟势的二阶导数即为动量的一阶导数, 因此当拟势的二阶导数发散时, 动量的一阶导数也发散, 这意味着拉格朗日流形的切空间斜率趋于无穷, 切空间指向动量的方向.

对于二维的拉格朗日流形而言, 只存在两种结构稳定的奇异性 (Whitney 1955, Arnold 1984, Smelyanskiy et al. 1997): 尖点(cusp)和折(fold). 具体示意图见图11(a) 所示. 从每一个尖点出发可以生成两个折, 当Hamilton轨迹经过折的时候, 拉格朗日流形切空间的斜率亦即动量的一阶导数趋于无穷. 将拉格朗日流形投影到坐标空间中, 折变成了焦散线, 如图11(b) 所示, Hamilton轨迹变成了相应的极值路径, 此时极值路径在坐标空间上看起来像是被焦散线"反射"回来一般. 两条焦散线之间的区域对应着拉格朗日流形的上、中、下三叶, 也就是说, 其中的每一个点都对应着三条Hamilton轨迹, 具有三个不同的动量. 因此, 该区域中的作用量也是多值的.

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11   拉格朗日流形拓扑结构奇异性示意图(Chen & Liu 2017b). (a)辅助Hamilton系统相空间中拉格朗日流形典型奇异性结构: 尖点(cusp)与折(fold); (b)与拉格朗日流形对应的坐标空间极值路径结构, 辅助Hamilton系统相空间中的折投影到二维平面变成了焦散线(caustic), 虚线表示切换线(switching line)


下面 讨论拉格朗日流形与最优离出路径的关系. 由式(8)可知, 作用量沿着Hamilton轨迹是不减的, 对于拉格朗日流形的上下两叶, 由于它们还没有越过折, 因而它们的作用量一定不超过已经越过折的Hamilton轨迹的作用量. 所以, 拉格朗日流形中叶的作用量最大 (Dykman et al. 1994; Smelyanskiy et al. 1997; Chen & Liu 2016, 2017a). 也正因为如此, 最优路径对应的Hamilton轨迹一定位于拉格朗日流形的上下两叶上, 因而最优路径一定不会触碰焦散线.

4.2.2 奇异性算例

下面结合几个典型的系统展示拉格朗日流形拓扑结构的奇异性. 结果如图12 所示. 这边 只是给出拉格朗日流形拓扑结构奇异性的直观感受, 并没有给出具体的系统方程和参数值, 详细内容请参考相应文献.

图12

新窗口打开| 下载原图ZIP| 生成PPT

图12   几个典型的拉格朗日流形拓扑结构奇异性系统. (a)周期激励单稳非线性振子 (Dykman et al. 1994), (b) IVDP振子 (Beri et al. 2005), (c)受简谐激励的过阻尼杜芬振子 (Dykman et al. 1996), (d) Maier-Stein系统 (Maier & Stein 1993a)


4.3 由拟势的不可微点集诱导的非光滑动力学行为

由4.2节可知拓扑结构中存在奇异性——尖点和折, 因而拉格朗日流形在局部具有三叶甚至多叶结构, 从而导致作用量函数在坐标空间的相应区域内产生多值性; 最优路径总是取得作用量函数的最小值, 然而, 由作用量的最小值所形成的曲面并非处处可微, 一般情况下, 它由若干张连续可微的曲面横截相交形成, 而在交集处, 其光滑性条件退化至连续. 数学上, 此分段可微的作用量最小值曲面, 对应了Hamilton-Jacobi方程的黏性解——拟势; 而在物理上, 拟势的不可微点集则对应了前面介绍的另一种奇异性: 切换线. 如前所述, 切换线两侧的点分别对应了拓扑性质完全不同的最优路径, 这也是其如此命名的原因; 在本节中, 针对一类具体的神经元模型——基于电导的二维Morris-Lecar系统, 讨论其分别在Type-I和Type-II可激性的条件下, 由弱噪声诱导的放电现象; 同时, 详细讨论切换线的存在对最优路径的动力学行为会产生何种影响. 可以发现, 物理学中定义的切换线这一概念, 在力学和动力系统的体系中, 同样有着重要的对应.

Morris-Lecar系统是通过对藤壶肌肉纤维中的电压振荡进行实验研究而建立的具有实际生物学意义神经元模型,形式如下

$ \left. \begin{array}{l} \dot {v} = I_{\rm app} - g_{\rm L} \left( {v - V_{\rm L} } \right) - g_{\rm K} \omega \left( {v - V_{\rm K} } \right) - g_{\rm Ca} m_\infty \left( v \right)\left( {v - V_{\rm Ca} } \right) + \xi _1 (t) \\ \dot {\omega } = \varphi \dfrac{\omega _\infty \left( v \right) - \omega }{\tau _\omega \left( v \right)} + \xi _2 (t) \\ \end{array}\right\}$

其中

$ \left. \begin{array}{l} m_\infty \left( v \right) = \dfrac{1}{2}\left[ {1 + \tanh \left( {\dfrac{v - V_1 }{V_2 }} \right)} \right]\\ \omega _\infty \left( v \right) = \dfrac{1}{2}\left[ {1 + \tanh \left( {\dfrac{v - V_3 }{V_4 }} \right)} \right]\\ \tau _\omega \left( v \right) = \left[ {\cosh \left( {\dfrac{v - V_3 }{2V_4 }} \right)} \right]^{ - 1} \\ \end{array}\right\}$

状态变量$v$表示细胞膜电位, $\omega $表示用以刻画K$^ + $电流活化程度的离子门变量; 方程(40)的第一式右侧依次包括: 外加电流刺激$I_{\rm app}$、欧姆泄漏电流、控制弛豫的电压门K$^ + $电流和控制兴奋的电压门Ca$^{2 + }$电流, 同时, 压敏电流Ca$^{2 + }$具有瞬时活化响应$m_\infty \left( v \right)$; 参数$g_{\rm L} $, $g_{\rm K} $和$g_{\rm Ca} $分别为三种跨膜电流的电导, $V_{\rm L} $, $V_{\rm K} $和$V_{\rm Ca} $分别是它们的反向电位(reversal potential); 参数$\varphi $决定K$^ + $通道状态变化的时间尺度. 噪声$\xi _1 (t)$表示突触噪声——神经元外部的随机输入, 而噪声$\xi _2 (t)$表示离子通道内部的热力学涨落, 满足: $\left\langle {\xi _i (t)} \right\rangle = 0,\;\left\langle {\xi _i (t)\xi _j \left( {t - \tau } \right)} \right\rangle = D\delta _{ij} \delta \left( \tau \right)$, 假设噪声强度$D \ll 1$. 本节重点讨论两种基本的相空间结构, 它们分别对应两种神经元的可激性类型, 即 Type-I可激与Type-II可激. 两种可激性相应的参数值取自计算神经科学领域的参考文献 (Ermentrout 1996, Lücken et al. 2013)

Type-I可激性

$\begin{array}{l} I = 0.05,\quad \phi = 0.2,\quad V_{\rm K} = - 0.7,\quad V_{\rm Ca} = 1,\quad V_{\rm L} = - 0.5,\quad g_{\rm L} = 0.5, \\ g_{\rm K} = 2,\quad g_{\rm Ca} = 1.33,\quad V_1 = - 0.01,\quad V_2 = 0.15,\quad V_3 = 0.1,\quad V_4 = 0.145 \\ \end{array}$

Type-II可激性

$\begin{array}{l} I = 0.05,\quad \phi = 0.2,\quad V_{\rm K} = - 0.7,\quad V_{\rm Ca} = 1,\quad V_{\rm L} = - 0.5,\quad g_{\rm L} = 0.5, \\ g_{\rm K} = 2,\quad g_{\rm Ca} = 1.33,\quad V_1 = - 0.01,\quad V_2 = 0.15,\quad V_3 = 0.017,\quad V_4 = 0.25 \\ \end{array}$

4.3.1 Type-I可激性与擦切环

利用OUM方法(Cameron 2012), 选取$30 000\times 30 000$的网格, 计算出系统(40)在有界区域$\left( {v,\omega } \right) \in \left[ { - 2,2} \right]\times \left[ { - 2,2} \right]$内拟势的分布情况; 沿着鞍点的稳定流形$W_{\rm s}$, 可以找到拟势在边界上的最小值点——鞍点$\pmb x_{\rm s}$, 这意味着, 在零噪声极限的意义下, 最优离出路径的离出位置正是鞍点. 从稳定结点$\pmb x_ {0}$出发的最优离出路径, 对应图13(a)中的绿色粗实线, 注意到, 在穿越边界$W_{\rm s}$之前, 最优离出路径并不沿着不稳定流形的左支运动, 这是因为非平衡态系统(40)不满足细致平衡条件, 因此不具备时间的反向对称性. 此外, 图13(a)中的蓝色箭头表示最优路径满足的向量场.

接下来,计算拟势在$v$零倾线右支上的分布, 如图13(a)的内图所示, 横坐标表示零倾线上某点到不稳定焦点$\pmb x_{\rm u}$的弧长; 同样地, 使用OUM和gMAM两种方法所得的结果非常一致, 保证了计算的可靠性. 由于最优离出路径在穿过边界$W_{\rm s}$后, 将沿着不稳定流形$W_{\rm u}$运动, 那么可想而知, 拟势在$v$零倾线右支上的最小值, 应当出现在$W_{\rm u}$与$v$零倾线的交点位置, 即图13(a)中的紫色正方形; 事实上, 根据大偏差理论, 这一结果也是显然的. 然而, 拟势的值沿着$v$零倾线呈阶梯状分布, 仔细研究这一分布的成因, 可以发现, 在不稳定焦点$\pmb x_{\rm u}$的周围, 出现了一个原系统中并不存在的极限环, 下面详细阐述这一现象.

图13

新窗口打开| 下载原图ZIP| 生成PPT

图13   Type-I可激性 (陈朕 2018). (a)最优路径及其向量场$\pmb G\left( \pmb x \right)$, 在不稳定焦点$\pmb x _{\rm u} $附近出现一个闭轨GC, 内图为 拟势在$v$零倾线右支上的分布; (b) GC附近的局部放大


在计算出拟势后, 可以得到最优路径满足的向量场

$ \dot {\pmb\varphi }_{\rm s} = \dfrac{1}{2}\nabla V\left( {\varphi _{\rm s} } \right) +\pmb l\left( {\varphi _{\rm s} } \right) = \nabla V\left( {\varphi _{\rm s} } \right) +\pmb F\left( {\varphi _{\rm s} } \right) \triangleq \pmb G\left( {\varphi _{\rm s} } \right)$

这里, $\pmb F\left( \pmb x \right)$表示Morris-Lecar系统的确定性向量场; 从任意点出发, 在$t \to + \infty $方向上的最优路径, 可通过向量场$\pmb G\left( \pmb x \right)$的流线获得; 反之, 可通过$ -\pmb G\left( \pmb x \right)$的流线获得. 在$v$零倾线的右支上选择若干个起点沿着向量场$ -\pmb G\left( \pmb x \right)$的流线, 可得分别以这些点为终点的最优路径, 如图13(a)中的黑色实线所示. 在$t \to - \infty $的方向上, 这些最优路径依次沿着稳定流形$W_{\rm s}$向上运动, 绕过不稳定焦点$\pmb x_{\rm u}$, 最后均沿着不稳定流形$W_{\rm u}$回到稳定结点$\pmb x_ {0}$. 但是, 在相反的演化方向$t \to + \infty $上, 从这些点出发的最优路径, 并不会像确定性系统一样, 完成一次放电, 然后弛豫至$\pmb x_ {0}$; 相反地, 如图13(a)中的浅蓝色实线所示, 它们在焦点附近趋于一个闭轨, 姑且称其为"GC", 即图中的红色实线, 具体含义将在后文解释. 通过OUM方法所得到的最优路径, 一旦接触GC, 将永远在GC上运动, 无法脱离, GC附近的局部放大如图13(b) 所示. 为验证GC的可靠性, 将起点取在稳定结点$\pmb x_ {0}$, 终点取在不稳定焦点$\pmb x_{\rm u}$, 通过gMAM方法, 得到了作用量的最小值曲线, 即图13(b)中的黄色实线; 可以明显地看出, 在螺旋地趋于焦点之前, 作用量最小值曲线先趋于GC、并沿GC运动了一段时间.

那么, GC是如何产生的? 要回答这一问题, 不妨先仔细讨论一下拟势沿$v$零倾线的阶梯状分布. 如图13(a)的内图所示, 将此分布划分为四个部分: Part I~Part IV. Part I是一个平台, 其边界恰好对应GC与零倾线的交点, 当噪声强度趋于零时, 所有涨落至Part I的最优路径. 理论上路径包括两个阶段. 第一阶段先趋于GC, 然后螺旋趋于终点. 由于趋于GC的最优路径几乎完全一致, 所以第一阶段将消耗几乎相同的作用量. 第二阶段从GC上剥离, 而后趋于焦点. 这个阶段将消耗额外的能量; 但是, 由于焦点是原系统的不变集, 此过程中的动量趋于零, 因此, 作用量最终趋于一个定值, 从而使得Part I的拟势变化幅度极小, 难以识别. Part II和 Part III对应GC以外的点, 它们的最优路径对应图13(a)的黑色实线. Part III中拟势也表现为一个平台, 主要有以下几方面原因: 第一, 稳定流形$W_{\rm s}$是原系统的不变集, 沿着它运动并不消耗作用量; 第二, 稳定流形两侧的确定性行为是高度发散的, 以至于从$W_{\rm s}$上剥离也几乎不需要消耗作用量; 第三, 最优离出路径从鞍点向$v$零倾线的右支运动时, $\pmb G\left( \pmb x \right) \approx \pmb F\left( \pmb x \right)$. 需要注意的是, 拟势的数值在Part I与Part II之间的变化, 是连续但不光滑的, 发生此突变的位置恰好对应了GC与零倾线的交点. 另一个拟势的突变点出现在Part III与Part IV之间, 经过该点后, 拟势的数值陡然下降. 事实上, 这一不可微点对应了零倾线上的另一个分水岭: 该点以下位置的最优路径, 当$t \to - \infty $时, 直接回溯至鞍点; 而该点以上位置的最优路径, 当$t \to - \infty $时, 先沿着稳定流形向上运动, 绕过焦点, 然后再沿不稳定流形回到鞍点, 如图13(a)最下方的两条黑色曲线所示. Part IV中的拟势, 描述了最优离出路径在穿过鞍点后, 从不稳定流形剥离所需消耗的作用量. 显然, 距离不稳定流形越近, 该值将越小, 故而最小值在不稳定流形的小邻域内取得.

根据以上分析可知, 拟势的不可微点, 对于理解GC的产生机理至关重要. 因此, 需要找到相空间中所有使得$\nabla V\left( \pmb x \right)$不连续的点, 即拟势$V\left( \pmb x \right)$的一阶偏导数不连续点. 事实上, 对于多元函数$V\left( \pmb x \right)$而言, 其偏导数连续是可微的充分条件; 那么$V\left( \pmb x \right)$的不可微点集, 一定包含于$\nabla V\left( \pmb x \right)$的不连续点集中, 如图14中的绿色粗实线所示, 记作$\varSigma _{\rm s} $. 拟势作为Hamilton-Jacobi方程的黏性解, 正是沿着此不可微点集"粘合"而成; 同时, 从物理意义上说, 这些不可微点集, 又对应了作用量最小值曲面中的切换 线——其中的每一点, 都同时对应两条最优路径以其为终点, 并且它们具有完全相同的作用量. 如图14 所示, $\varSigma _{\rm s} $中的每一点都对应两条最优路径(图中黑色曲线), 而拟势的连续性则保证了它们的作用量必然相等. 最后需要指出, 拟势的不可微点集, 与闭轨GC近乎相切于一点, 即图14中的蓝色五角星, 对此, 接下来做详细说明.

图14

新窗口打开| 下载原图ZIP| 生成PPT

图14   Type-I可激性的条件下, 拟势的不可微点集$\varSigma _{\rm s} $ (陈朕 2018). 其中每一点都对应两条最优路径以其为终点(黑色实线); 内图中 不可微点集$\varSigma _{\rm s} $与闭轨GC近乎相切


从图14可以看出, 由式(42)所确定的最优路径在接触切换线后, 并不会穿过切换线, 而是沿着切换线$\varSigma _{\rm s} $运动. 其多值性导致了以下事实: 当一条最优路径穿过切换线时, 它将不再是最优或作用量最小的, 因为在作用量曲面中, 它已穿过了两片最小作用量曲面横截相交位置处的交线, 并位于其中一片最小值曲面的上方. 因此, 从物理意义上讲, 最优路径不论以何种角度运动至切换线, 它都不会穿过, 而将骤然地、不连续地改变方向.

另一方面, 回到方程(42), 从数学上讲, 最优路径的动力学行为在切换线处的非光滑性, 也正是源于拟势的梯度在切换线处的不连续性; 换句话说, 拟势在切换线处的不可微性, 诱导了一个分段光滑的动力系统(42), 该系统含有一个一维的不连续性边界 (Kuznetsov et al. 2003) (discontinuity boundary)

$ \varSigma = \left\{ {\pmb x \in {\rm R}^2\vert h\left( \pmb x \right) = 0} \right\}$

使得在$\varSigma $的两侧, 向量场$\pmb G\left( \pmb x \right)$具有不同的动力学行为: $\pmb G^{(i)}\left( \pmb x \right)$, $i = 1$, 2. 根据Filippov的构造, 将不连续性边界$\varSigma $的子集称为滑移集 (Kuznetsov et al. 2003) (sliding set), 若

$ \varSigma _{\rm s} = \left\{ {\pmb x \in \varSigma \vert \sigma \left( \pmb x \right) \leq 0} \right\}$

其中

$ \sigma \left( \pmb x \right) = \left\langle {\nabla h\left( \pmb x \right),\pmb G^{\left( 1 \right)}\left( \pmb x \right)} \right\rangle \left\langle {\nabla h\left( \pmb x \right),\pmb G^{\left( 2 \right)}\left( \pmb x \right)} \right\rangle$

反之, 若$\sigma \left( \pmb x \right) > 0$, 则称之为穿越集(crossing set). 对于滑移集中的点, 若满足

$ \left\langle {\nabla h\left( \pmb x \right),\pmb G^{\left( 2 \right)}\left( \pmb x \right) - \pmb G^{\left( 1 \right)}\left( \pmb x \right)} \right\rangle = 0$

则称为奇异滑移点 (Kuznetsov et al. 2003) (singular sliding points); 在这些点处, 必有以下事实之一成立: 或者向量场$\pmb G^{\left( 1 \right)}\left( \pmb x \right)$和$\pmb G^{\left( 2 \right)}\left( \pmb x \right)$均切于$\varSigma _{\rm s} $, 或者其中之一为零、而另一个切于$\varSigma _{\rm s} $, 或者两者均为零.

根据Filippov方法, 对于滑移集$\varSigma _{\rm s} $中的点,它的动力学行为满足

$ \dot{\pmb x} = \pmb g \left( \pmb x \right)$

其中, $\left\langle {\nabla h\left( \pmb x \right),\pmb g \left( \pmb x \right)} \right\rangle = 0$, 而$\pmb g \left( \pmb x \right)$是向量场$\pmb G^{(i)}\left( \pmb x \right)$的凸线性组合

$ \pmb g \left( \pmb x \right) = \lambda \pmb G^{\left( 1 \right)}\left( \pmb x \right) + \left( {1 - \lambda } \right)\pmb G^{\left( 2 \right)}\left( \pmb x \right),\qquad \lambda = \dfrac{\left\langle {\nabla h\left( \pmb x \right),\pmb G^{\left( 2 \right)}\left( \pmb x \right)} \right\rangle }{\left\langle {\nabla h\left( \pmb x \right),\pmb G^{\left( 2 \right)}\left( \pmb x \right) - \pmb G^{\left( 1 \right)}\left( \pmb x \right)} \right\rangle }$

这意味着$\pmb g \left( \pmb x \right)$是切于滑移集$\varSigma _{\rm s} $的, 如图15 所示.

图15

新窗口打开| 下载原图ZIP| 生成PPT

图15   Filippov构造下的非光滑系统示意图(陈朕 2018). (a)滑移集, (b)擦切环, (c)滑移环


回到方程(42), 在时间正向演化的方向上, 切换线两侧的最优路径均朝着$\varSigma _{\rm s} $运动, 满足滑移集的定义. 此外, 如图14的内图所示, 滑移集的末端与GC的一小部分重合, 从而迫使附近的最优路径运动至GC上; 同时, 在它们的交集处(图中蓝色五角星), 向量场$\pmb G\left( \pmb x \right)$与滑移集几乎相切; 若成立, 则该点对应式(46)所定义的奇异滑移点, 而GC则对应平面Filippov系统中, 擦切分岔或滑移分岔 (Guardia et al. 2011)导致的擦切环(grazing cycle)或滑移环(sliding cycle), 如图15(b)和图15(c) 所示.

最后, 需要强调的是, 在系统(42)中, 上述的擦切分岔或者滑移分岔是否发生, 以及在系统参数取何值时发生, 都需要对非光滑动力系统(42)进行进一步研究后确定.

4.3.2 Type-II可激性与离出边界的选择

下面考察Type-II可激性条件下的Morris-Lecar系统. 由于相空间中不存在鞍点及其稳定流形作为边界, 因此在这种情况下, 系统不存在明确定义的阈值流形. 一个合理的阈值流形, 应当满足条件: 它是系统(40)的一条特殊的解轨线, 且在它的两侧, 动力学行为是高度发散的. 一般地, 考察如下的奇异摄动系统

$ \dot {v} = f\left( {v,\omega ,\varepsilon } \right),\qquad \dot {\omega } = \varepsilon g\left( {v,\omega ,\varepsilon } \right)$

其中$\left( {v,\omega } \right) \in {\rm {\bf R}}^k\times {\rm R}^m$为系统的状态变量, 假设函数$f$和$g$为$C^\infty $光滑的, 并且$\varepsilon \ll 1$. 当$\varepsilon = 0$时, 上述条件可由系统(49)的临界流形(critical manifold)$S \triangleq \left\{ {\left( {v,\omega } \right)\vert f\left( {v,\omega ,\varepsilon } \right) = 0} \right\}$的法向双曲子集满足; 然而, 法向双曲性(normal hyperbolicity)在临界流形$S$的折点$\pmb{x}^{\ast }$处失效. 所谓折点, 即$f$关于快变量$v$的雅各比矩阵$D_v f\left( {\pmb x^\ast } \right)$至少有一个特征值具有零实部, 这意味着对不同的慢变量$\omega $, 式(49)的子系统产生鞍结分岔. 集合$S$即为$v$零倾线的中间一支, 其左右两侧的动力学行为高度发散, 折点$\pmb{x}^{\ast }$就是位于零倾线右侧的极大值点, 如图16的星号所示.

图16

新窗口打开| 下载原图ZIP| 生成PPT

图16   Type-II可激性条件下, 拟势的不可微点集,最优路径及其向量场 (陈朕 2018)


当$\varepsilon \neq 0$时, 根据Fenichel(1979)理论, 全系统中的法向双曲流形得以保持(要求$\varepsilon $足够小), 并且阈值现象与$S$的折奇异性密切相关; 因此, 分界线可由奇异摄动问题(49)的一组特殊解——Canard轨线给出, 并且, 对一般的Morris-Lecar神经元模型, 大量研究表明, Canard轨线可以为不同的动力学行为提供合适的边界 (Wechselberger et al. 2013). 将分界线(separatrix)标记为图16中黄色粗实线, 可以计算出在Type-II可激性条件下, Morris-Lecar系统的拟势分布; 同样地, 也能找到拟势在分界线上的最小值点, 如图16中的蓝色正方形所示, 它是零噪声极限意义下的穿越分界线的最优离出位置, 以它为终点的最优离出路径, 对应图中的红色曲线; 当终点在分界线上移动时, 从平衡点$\pmb{x}_{0}$出发的最优路径(图中蓝色实线)逐渐靠近最优离出路径. 另一方面, 切换线为图中的绿色粗实线, 内图则给出了极值哈密顿轨线向坐标平面的$\left( {v,\omega } \right)$的投影, 可以清楚地看出, 两条焦散线从尖点(蓝色五角星)出发, 且切换线位于两焦散线之间. 需要强调的是, 内图中的哈密顿轨线只取得作用量泛函的局部极值, 而泛函的全局最小值, 即最优路径, 则对应图中的黑色实线所示, 同样, 切换线上每一点都对应两条最优路径.

下面, 对参考文献 (Newby et al. 2013)中给出的离出边界加以讨论, 并将其与Canard轨线确定的分界线进行比较. 在参考文献 (Newby et al. 2013)中, 离出边界由切换线和拟势通过尖点的等高线给出, 如图中的青色实线所示. 经过计算后可以发现, 在切换线上, 拟势的最小值点确实在尖点处取得, 那么在此意义下, 过尖点的等高线刻画了平衡点附近的一个势阱, 动作电位的产生则相当于从此势阱中离出; 此外, 参考文献中所定义的"bottleneck", 实际就是拟势的等高线与最优离出路径的交点, 对应图中的蓝色三角形. 本质上, 两种取法的出发点不同, 文献 (Newby et al. 2013)从能量的观点出发, 选择拟势的等高线作为离出边界; 而我们 则从动力学的角度出发, 试图找到一条阈值流形, 穿过它后, 系统必定会发生一次典型的动作电位. 事实上, 从拟势的等高线离出后, 许多最优路径只产生了小幅的瞬态响应, 并未完成典型的放电, 因此, 从动力学的观点看, 拟势的等高线并不能作为确定放电与否的边界. 除此之外, 从图中也可以看出, 等高线的一部分与切换线极为接近, 直到切换线上的某点为止, 该点对应图中标记为"Canard"的轨线, 且该轨线经过零倾线的右侧极值点; 由于作用量沿着最优路径是不减的, 因此, 这条"Canard"轨线使得拟势在分界线的左侧区域非常平坦, 以致于上述等高线所确定的势阱边沿, 可被大范围延伸、且数值几乎不变: 上至"Canard"轨线, 右至分界线. 因此, 即便从能量的角度出发, 选择分界线作为边界也是合理的.

4.4 有限噪声强度

大偏差理论描述的最优离出路径要求噪声强度趋于零, 因而通过作用量泛函或者拟势所刻画的最优离出路径也是在零噪声极限意义下的. 比如, IVDP系统中最优离出路径是一条连接平衡点与极限环的异宿轨迹. 但是实际系统中或实验中, 噪声强度总是有限的, 这种情况下, 理论上的最优离出路径和实际结果有可能存在偏差, 产生诸如噪声导致的偏移(noise-induced shift) (Bandrivskyy et al. 2003a)、噪声导致的鞍点避免(noise-induced saddle-point avoidance) (Luchinsky et al. 1999)等现象. 这种情况下, 零噪声极限下的作用量近似式(8)还不够, 需要考虑更高阶的近似. 在WKB近似摄动展开中

$ p_1 \left(\pmb x \right)=\exp \left( { - \dfrac{1}{\varepsilon }\left( {S_0 \left(\pmb x \right) + \varepsilon S_1 \left(\pmb x \right)+o\left( \varepsilon \right)} \right)} \right) \triangleq z\left(\pmb x \right)\exp \left( { - \dfrac{1}{\varepsilon }\left( {S_0 \left(\pmb x \right) + o\left( \varepsilon \right)} \right)} \right)$

因而, 对于作用量泛函, 考虑如下弱噪声强度的一阶近似 (Bandrivskyy et al. 2003a, Beri et al. 2004)

$ S_\varepsilon \left(\pmb x \right) = S_0 \left(\pmb x \right) - \varepsilon \ln \left( {z\left(\pmb x \right)} \right)$

式中, $S_0 \left(\pmb x \right)$为零噪声极限的作用量, $z\left(\pmb x \right)$为WKB近似中的指数前因子项, $\varepsilon $为噪声强度. 由此, 得到了噪声强度一阶近似的修正作用量泛函$S_\varepsilon \left(\pmb x \right)$, 它的全局极小值对应给定弱噪声强度$\varepsilon $下的最优离出路径和离出点.

根据式(51), 要计算修正后的作用量, 需要计算出Hamilton轨迹的指数前因子项, 在前一节的讨论中, 提到当Hamilton轨迹越过折的时候, 拟势的二阶导数即动量的一阶导数会发散, 相应的指数前因子项也会发散, 这给式(51)的计算带来了麻烦. 幸运的是, 由于最优离出路径不会触碰焦散线, 因而当 计算最优离出路径的前因子项时不会出现发散的情形. 所以采用式(51)计算有限噪声强度下的最优离出路径和离出位置是有效的.

接下来 以受简谐激励的过阻尼杜芬振子为例, 研究噪声导致的偏移现象,系统方程如下 (Beri et al. 2005)

$ \dot {x} = x - x^3 + A\sin \left( {\omega t} \right) + \sqrt \varepsilon \xi (t)$

将其化为自治系统

$ \left.\begin{array}{l} \dot {x}_1 = \omega \\ \dot {x}_2 = x_2 - x_{2} ^3 + A\sin \left( {x_1 } \right) + \sqrt \varepsilon \xi (t) \\ \end{array}\right\}$

参数选取$A = 0.1$, $\omega = 1$. 通过action plot方法, 可以得到系统的最优离出路径, 如图17绿色粗实线所示. 根据式(51) 对作用量进行修正, 可以得到有限噪声强度的作用量结果, 针对噪声强度$\varepsilon = 0.01$, $0.05$, 相应的最优离出路径如图17中红色虚线和紫色点划线所示. 从图上可以看出, 有限噪声强度的最优路径不再和不稳定极限环相切, 而是以横截的形式离出.

图17

新窗口打开| 下载原图ZIP| 生成PPT

图17   受简谐激励的过阻尼杜芬振子最优离出路径及有限噪声强度修正. 红色和蓝色实线分别为稳定和不稳定极限环, 黑色实线为涨落路径, 绿色粗实线为零噪声极限最优离出路径, 红色三角和黄色虚线分别为尖点和焦散线, 紫色点划线和红色虚线分别为噪声强度$\varepsilon = 0.05$, 0.01 时修正后的最优离出路径


5 离出位置分布

离出位置分布反映了噪声扰动的离出路径在边界上的离出点的分布情况. 这一问题呈现在许多科技领域内, 包括统计物理、化学物理、流体力学、生态学、通信等 (Luchinsky et al. 1999). 对于具有鞍点的平衡态系统, 随着噪声强度趋于零, 边界上的离出位置分布一般情形下渐近地收敛于以鞍点为中心的高斯分布, 分布的宽度随着噪声强度的减弱而减小. 而对于非平衡态系统(non-equilibrium system), 边界上的离出位置分布一般将不再是高斯的, 而会发生偏斜(skewed), 理论分析、数值模拟和实验都验证了这一点 (Maier & Stein 1993a, 1997; Luchinsky et al. 1999). 此外, 正如4.4节中所看到的, 对于有限的噪声强度扰动, 离出点的位置会发生偏移, 相应的离出位置分布也会产生一定的变化.

离出位置分布一般情况下很难解析求解, 对于复杂的边界更是如此. 国内外的学者研究了各种离出边界情形下的离出概率. 如Cai等 (2017)研究了列维噪声激励下可激系统的离出问题, 数值计算了矩形边界的平均首次离出时间及边界上的首次离出概率. Marshall (2016)研究了带有任意缺口的圆盘边界离出问题, 得出了首达时间和分离概率的解析结果. 求解离出位置分布最直接的方法是采用蒙特卡罗模拟进行蛮力(brute force)求解; 但是由第3节可知, 随着噪声强度的减弱, 离出轨迹的平均首次离出时间将会呈指数级地增大, 模拟过程大部分都浪费在吸引子周围无意义的随机游走上. 为此, 当实际问题不需要考虑最优离出路径和平均首次离出时间时, 需要找出加快模拟的方法.

近年来, 提高模拟效率的数值算法层出不穷, 表1 中列举了部分比较有代表性的方法. 其中前四个算法 (Faradjian & Elber 2004, Allen et al. 2005, Glowacki et al. 2009, Adams et al. 2010)是分子动力学(molecular dynamics, MD)中的一些著名算法, 采用了中间截面, 将离出过程分段, 从而有效缩短了样本离出的时间. DIMS方法 (Beri et al. 2004)将离出样本控制在最优离出路径的稳定流形上. 而有效转移路径采样方法(Crooks & Chandler 2001)则更加直接地模拟噪声过程, 对于非平衡随机动力学模拟更有效. 此外, 广义胞映射方法 (Han et al. 2016)作为一种研究非线性系统全局分 析的有效数值工具, 对于离出过程的数值模拟计算也取得了非常好的效果. 最近, 概率演化算法 (Zhu et al. 2018)综合了FFS方法以及快速蒙特卡罗模拟方法 (Bandrivskyy et al. 2003b)的优点, 对于弱噪声情形下的离出问题非常有效. 这里 对于概率演化算法做简单介绍.

5.1 概率演化算法

概率演化算法(probability evolution method)的示意图如图18 所示. 受到分子动力学模拟方法的启示, 在离出过程的初始位置和最终边界之间采取一系列截面, 然后将模拟过程分段, 在每一对相邻截面之间, 在前一个截面取一定数量的出发点, 然后模拟从这些点出发终止于后一个截面的样本轨迹, 最终可以得到从这些点出发终止于后一个截面的转移概率分布. 为了加快模拟, 不允许样本点越过前一个边界, 如右图中虚线表示的箭头, 每一条跨越前一个截面的样本都会迅速地返回起点, 这一过程称之为回注(reinjection). 通过上述回注的方法, 文献 (Bandrivskyy et al. 2003b)指出从一个点到下一个截面或者边界所需要的时间将会被指数级地减小. 对于任意两个相邻截面, 取初始截面中若干个起始点, 计算每一个点到下一个截面的转移概率分布, 如图18中第1个截面处的第$i$个点$x_i^{\left( 1 \right)} $和第$j$个点$x_j^{\left( 1 \right)} $, 得到他们转移到下一个截面的转移概率分布为图中绿色和红色色块. 最终, 可以得到两个截面之间的转移概率分布矩阵$P\left( {\pmb x^{\left( {i + 1} \right)}\vert \pmb x^{(i)}} \right)$, 其中$\pmb x^{(i)}$表示第$i$个截面上选取的点的集合, 因而矩阵中第$m$行表示$\pmb x^{(i)}$中的第$m$个点到达第$i+1$个截面上的转移概率分布. 假设不同截面之间样本的独立性, 从而通过类似马氏链的形式, 可以得到最终的离出位置分布

$ P\left( {\pmb x^{\rm (end)}\vert \pmb x^{(0)}} \right) = \prod\limits_{i = 0}^n {P\left( {\pmb x^{\left( {i + 1} \right)}\vert \pmb x^{(i)}} \right)}$

其中, $n$表示截面总数, $i=0$和$i=n+1$分别代表初始和最终位置.

表1   随机动力学模拟的部分快速算法

新窗口打开| 下载CSV


图18

新窗口打开| 下载原图ZIP| 生成PPT

图18   概率演化算法示意图(Zhu et al. 2018)


与现有的计算离出位置分布的算法相比, 概率演化算法最大的优点是可以进行大规模的并行计算, 这有赖于不同截面之间计算的独立性, 以及同一截面上不同起始点之间计算的独立性, 因而可以大大削减计算时间, 这与FFS等方法有着本质的区别. 但是, 由于采用了回注的技巧, 当噪声强度较大时, 会产生较大的误差, 因为忽略了穿越前一个截面又重新转移到下一个截面的那些样本轨迹, 这是概率演化算法的一个主要局限性. 但由大偏差理论知道, 当噪声强度充分小时, 那些不沿着最优离出路径的样本轨迹将会被指数级地削减, 而考虑到采用该算法的目的, 当噪声强度足够小时, 误差是可以容忍的. 此外, 由于计算利用到了马氏性, 因而, 对于色噪声(或者说时间相关噪声), 若要采用这里讨论的方法, 需要保证相邻截面之间的平均首次离出时间远大于噪声的相关时间(此时可以采用Markov近似, 不同截面之间近似具有马氏性), 但误差仍然存在, 所以对于色噪声要慎用本方法.

5.2 算例

下面考察两个经典的例子, 验证概率演化算法计算离出位置分布的有效性.

5.2.1 Maier-Stein系统

系统的方程如式(35)所示, 为了与解析解做对比, 仍然采用4.1.3节中的参数, 分别为$\alpha = 6.67$, $\mu = 2$以及$\alpha = 10$, $\mu = 0.67$, 相应的最优离出路径如图8中红色虚线所示. 对于第一种参数情形, 最优离出路径沿着$x$轴, 这种情形类似于对称性未破坏的情况 (Maier & Stein 1993b), 离出位置分布具有以原点为中心的高斯分布 (Maier & Stein 1997)

$ P(y) = \left( {\dfrac{2}{\pi D}} \right)^{1 / 2}\exp \left( { - \dfrac{2y^2}{D}} \right)$

而对于第二种情形, 最优离出路径发生了对称性破缺, 为避免触碰焦散线, 路径偏离了$x$轴, 最终仍然从鞍点离出. 而对于有限噪声情形, 此时离出位置分布形式发生了偏斜, 转变为以原点为中心的对称韦伯分布(Weibull distribution) (Maier & Stein 1993a, 1997)

$ P(y) = \left( {\dfrac{2}{A^{2 / \mu }D\mu }} \right)\left| y \right|^{2 / \mu - 1}\exp \left( { - \left( {\left| y \right| / A} \right)^{2 / \mu } / D} \right)$

式中, $A$为$y / x^\mu $趋于原点时的极限. 两种参数情形下的理论离出位置分布如图19中实线所示. 图中符号为概率演化算法计算的结果. 对于第一种参数情形, 从图19(a)可以看出, 随着噪声强度的减小, 分布的宽度减小而高度增加了. 这与弱噪声极限情形时最优离出路径将会穿过原点的理论结果一致, 离出位置也会相应地集中在鞍点附近. 概率演化算法得到的离出位置分布与式(55)给出的理论结果拟合较好. 对于第二种参数情形, 由图19(b)可以看出, 最优离出路径将会渐进地沿着$y$轴趋于原点. 由WKB近似可知, 由于最优离出路径与$y$轴相切, WKB的管道将会在通过原点前与$y$轴相碰, 从而产生了鞍点避免现象 (Luchinsky et al. 1999). 这边需要注意的是, 式(56)为噪声强度趋于零时候的结果, 对于有限噪声强度会有偏差, 因而图19(b)中的理论与数值解的偏差可以理解, 当噪声强度减小时, 理论与数值结果会更加吻合.

图19

新窗口打开| 下载原图ZIP| 生成PPT

图19   Maier-Stein系统两种参数条件下离出位置分布 (Zhu et al. 2018). 图中实线为理论结果, 符号为概率演化算法的数值计算结果. (a) $\alpha = 6.67$, $\mu = 2$. 噪声强度分别为: $D=0.01$ (红色), 0.03 (绿色), 0.05 (蓝色); (b) $\alpha = 10$, $\mu = 0.67$. 噪声强度分别为: $D =0.015$ (红色), 0.04 (绿色)


5.2.2 IVDP系统

在4.1.2节中, 采用action plot方法计算了IVDP系统的最优离出路径, 它是一条连接平衡点和极限环的异宿轨迹, 对于有限噪声强度的情形, 最优离出路径将会发生偏移 (Beri et al. 2004). 本节采用4.4节中的作用量修正方法作为理论值与概率演化算法的数值结果做对比, 研究IVDP系统中噪声导致的离出路径偏移现象, 结果如图20示. 当噪声强度较小的时候$(D=0.05)$, 理论与数值计算结果吻合地较好, 而当噪声强度较大时$(D=0.3$和$D=0.5)$, 它们之间的差异会增大. 造成这种误差的来源在5.1节中已经做过讨论, 这是回注导致的直接结果. 幸运的是, 这种误差会随着噪声强度的减小而减小, 这与 使用概率演化算法的初心相一致.

图20

新窗口打开| 下载原图ZIP| 生成PPT

图20   IVDP系统噪声导致的最优离出路径偏移 (Zhu et al. 2018). 图中黑色粗实线为极限环, 灰色细实线为概率演化算法中采用的截面, 绿色实线、蓝色虚线和紫色点划线分别对应噪声强度$D=0.05$, 0.3, 0.5的理论修正后的最优离出路径, 相应颜色的符号为概率演化算法得到的离出位置分布峰值


6 展望

大偏差理论是求解离出问题的强大工具, 随着人类对于自然的描述更加精细, 实际问题对于大偏差理论的应用范围提出了很多新的挑战. 从以下几个方面提出一些潜在的开放性问题:

(1) 由于低维系统的常见性以及易于可视化的特点, 大偏差理论被广泛地应用于低维系统, 尤其是二维系统. 而随机微分方程理论表明 (Bernt Øksendal 2010), 二维布朗运动具有常返性(recurrence), 而三维及以上布朗运动不具有常返性, 因而可以预见高维动力系统离出问题除了计算上的复杂度, 还有现象上的丰富性.

(2) 边界的复杂性也会给离出问题带来丰富的行为. 比如窄口逃逸问题(narrow escape problem) (Holcman & Schuss 2004), 窄口逃逸问题在数学上描述为: 一个布朗粒子(比如离子、分子、蛋白质)被限制在一个具有反射边界的有界区域内(比如细胞内), 只有一个小窗口可供逃逸, 这一问题在细胞生物学中具有巨大的应用前景 (Schuss et al. 2007, Bressloff & Newby 2013, Agranov & Meerson 2018). 此外, 分形边界也给离出问题带来不小的挑战 (Chen et al. 2016, Chen & Liu 2017b).

(3) 从3.2节的讨论中可以发现时间尺度对于噪声导致的离出具有重要影响, 噪声在不同时间尺度上的作用并不如通常想象的那样平凡 (Lee DeVille et al. 2005), 由于如神经元系统在内的诸多实际系统中普遍存在多时间尺度, 因而研究噪声在快慢动力系统中的作用具有非常重要的现实意义 (朱金杰 2018).

(4) 尽管在数学处理上, 白噪声常常被用来作为随机扰动的典型, 但是实际系统中的噪声往往更加复杂, 比如列维噪声(Lévy noise) (Cai et al. 2017)、泊松噪声(Poisson noise) (Dykman 2010)、色噪声 (Schuecker et al. 2015)、分数阶噪声 (Wang et al. 2017)等. 今后的研究中, 要建立更加符合实际的噪声模型, 更好地应用大偏差理论来分析实际问题.

(责任编委: 丁千)

致谢

国家自然科学基金资助项目 (11772149, 11472126).

参考文献

陈朕 . 2018.

基于大偏差理论的几类非线性随机系统动力学行为研究. [博士论文]

南京: 南京航空航天大学

URL     [本文引用: 4]

( Chen Z . 2018.

Dynamical behaviors of several nonlinear stochastic systems based on the large deviation theory

[PhD Thesis]. Nanjing: Nanjing University of Aeronautics and Astronautics).

URL     [本文引用: 4]

孔琛 . 2018.

噪声扰动下非线性动力系统离出行为研究. [博士论文]

南京: 南京航空航天大学

[本文引用: 1]

( Kong C . 2018.

On the exit problems in nonlinear dynamical systems driven by random perturbations

[PhD Thesis]. Nanjing: Nanjing University of Aeronautics and Astronautics).

[本文引用: 1]

孔琛, 刘先斌 . 2014.

受周期和白噪声激励的分段线性系统的吸引域与离出问题研究

力学学报, 46:447-456

DOI      URL     [本文引用: 2]

离出行为是随机非线性系统的重要现象之一，而离出问题是除随机动力系统理论以外考察随机非线性系统随机稳定性的另一种重要的方法.分段线性系统是一个经典的非线性动力学模型，受随机激励后成为随机系统，但并不是严格的随机动力系统，因而此时随机动力系统理论也不适用.为了研究同时受周期和白噪声激励的分段线性系统，首先使用Poincar&#233;截面模拟其在无噪声时确定性的动力学行为，然后使用Monte Carlo模拟对其在白噪声激励下的离出行为进行了数值仿真分析.其次，为了考察离出问题中的重要参数，系统的平均首次通过时间（mean first-passage time，MFPT），使用van der Pol变换，随机平均法，奇异摄动法和射线方法进行了量化计算.通过对理论结果与模拟结果的对比分析，得到结论：当系统吸引子对应的吸引域边界出现碎片化时，理论结果与模拟结果的误差极大；而当吸引域边界足够光滑的以后，理论结果与模拟结果才会相当吻合.

( Kong C, Liu X B . 2014.

Research for attracting region and exit problem of a piecewise linear system under periodic and white noise excitations

Chinese Journal of Theoretical and Applied Mechanics, 46: 447-456).

DOI      URL     [本文引用: 2]

离出行为是随机非线性系统的重要现象之一，而离出问题是除随机动力系统理论以外考察随机非线性系统随机稳定性的另一种重要的方法.分段线性系统是一个经典的非线性动力学模型，受随机激励后成为随机系统，但并不是严格的随机动力系统，因而此时随机动力系统理论也不适用.为了研究同时受周期和白噪声激励的分段线性系统，首先使用Poincar&#233;截面模拟其在无噪声时确定性的动力学行为，然后使用Monte Carlo模拟对其在白噪声激励下的离出行为进行了数值仿真分析.其次，为了考察离出问题中的重要参数，系统的平均首次通过时间（mean first-passage time，MFPT），使用van der Pol变换，随机平均法，奇异摄动法和射线方法进行了量化计算.通过对理论结果与模拟结果的对比分析，得到结论：当系统吸引子对应的吸引域边界出现碎片化时，理论结果与模拟结果的误差极大；而当吸引域边界足够光滑的以后，理论结果与模拟结果才会相当吻合.

刘先斌, 陈虬, 陈大鹏 . 1996.

非线性随机动力系统的稳定性和分岔研究

力学进展, 26:437-452

URL     [本文引用: 1]

在随机动力系统中的分岔──噪声导致的跃迁行为，是一种有别于确定性系统分岔与混沌的独特的非线性复杂现象．本文全面评述非线性随机系统的稳定性问题、离出问题、随机动力系统理论和随机分岔等各项研究的发展历史、基本的思想方法以及主要的研究成果．

( Liu X B, Chen Q, Chen D P . 1996.

The researches on the stability and bifurcation of nonlinear stochastic dynamical systems

Advances in Mechanics, 26: 437-452).

URL     [本文引用: 1]

在随机动力系统中的分岔──噪声导致的跃迁行为，是一种有别于确定性系统分岔与混沌的独特的非线性复杂现象．本文全面评述非线性随机系统的稳定性问题、离出问题、随机动力系统理论和随机分岔等各项研究的发展历史、基本的思想方法以及主要的研究成果．

孙建桥, 熊夫睿 . 2017.

非线性动力学系统全局分析之外的胞映射方法新发展

力学进展, 47:150-177

URL     [本文引用: 1]

( Sun J Q, Xiong F R . 2017.

Cell mapping methods-beyond global analysis of nonlinear dynamic systems

Advances in Mechanics, 47: 150-177).

URL     [本文引用: 1]

徐伟, 孙春艳, 孙建桥, 贺群 . 2013.

胞映射方法的研究和进展

力学进展, 43:91-100

DOI      URL     [本文引用: 1]

介绍了胞映射方法的研究和进展. 归纳了目前胞映射方法的几种主要研究方法, 主要包括简单胞映射、广义胞映射、图胞映射、图胞映射的符号分析方法、图胞映射的面向集合方法、邻接胞映射、庞加莱型的简单胞映射、插值胞映射以及胞参照点映射方法, 分析了各类方法的基本特点和特色, 简述了这几种胞映射方法的最新国内外进展, 综述了胞映射方法在控制及相关领域的应用研究及进展, 给出了胞映射方法研究的一些展望, 提出了胞映射方法研究可能率先突破的几个研究方向.

( Xu W, Sun C Y, Sun J Q, He Q . 2013.

Development and study on cell mapping methods

Advances in Mechanics, 43: 91-100).

DOI      URL     [本文引用: 1]

介绍了胞映射方法的研究和进展. 归纳了目前胞映射方法的几种主要研究方法, 主要包括简单胞映射、广义胞映射、图胞映射、图胞映射的符号分析方法、图胞映射的面向集合方法、邻接胞映射、庞加莱型的简单胞映射、插值胞映射以及胞参照点映射方法, 分析了各类方法的基本特点和特色, 简述了这几种胞映射方法的最新国内外进展, 综述了胞映射方法在控制及相关领域的应用研究及进展, 给出了胞映射方法研究的一些展望, 提出了胞映射方法研究可能率先突破的几个研究方向.

徐伟, 岳晓乐, 韩群 . 2017.

胞映射方法及其在非线性随机动力学中的应用

动力学与控制学报, 15:200-208

URL     [本文引用: 1]

( Xu W, Yue X L, Han Q . 2017.

Cell mapping method and its applications in nonlinear stochastic dynamical systems

Journal of Dynamics and Control, 15: 200-208).

URL     [本文引用: 1]

许勇, 裴斌, 徐伟 . 2017.

随机平均原理研究若干进展

动力学与控制学报, 15:193-199

URL     [本文引用: 1]

( Xu Y, Pei B, Xu W . 2017.

Some recent developments of stochastic averaging principle

Journal of Dynamics and Control, 15: 193-199).

URL     [本文引用: 1]

朱位秋 . 1987.

随机平均法及其应用

力学进展, 17:342-352

URL     [本文引用: 1]

( Zhu W Q . 1987.

Stochastic averaging methods and their applications

Advances in Mechanics, 17: 342-352).

URL     [本文引用: 1]

朱位秋 . 1992. 随机振动. 北京: 科学出版社

[本文引用: 1]

( Zhu W Q. 1992. Stochastic Vibration. Beijing: Science Press).

[本文引用: 1]

朱位秋 . 2003. 非线性随机动力学与控制——Hamilton理论体系框架. 北京: 科学出版社

[本文引用: 2]

( Zhu W Q. 2003. Nonlinear Stochastic Dynamics and Control—the Framework of Hamilton Theory. Beijing: Science Press).

[本文引用: 2]

朱位秋, 蔡国强 . 2017. 随机动力学引论. 北京: 科学出版社

[本文引用: 1]

( Zhu W Q, Cai G Q. 2017. Introduction to Stochastic Dynamics. Beijing: Science Press).

[本文引用: 1]

朱金杰 . 2018.

神经元同步、共振及离出问题研究. [博士论文]

南京: 南京航空航天大学

[本文引用: 6]

( Zhu J J . 2018.

Synchronization, resonance and exit problem for neuronal dynamical systems

[PhD Thesis]. Nanjing: Nanjing University of Aeronautics and Astronautics).

[本文引用: 6]

Adams D A, Sander L M, Ziff R M. 2010.

The barrier method: A technique for calculating very long transition times

Journal of Chemical Physics, 133:124103.

DOI      URL     PMID      [本文引用: 1]

In many dynamical systems, there is a large separation of time scales between typical events and

Agazzi A, Dembo A, Eckmann J P. 2017.

Large deviations theory for Markov jump models of chemical reaction networks

Annals of Applied Probability, 28:1821-1855.

DOI      URL     [本文引用: 1]

Agranov T, Meerson B. 2018.

Narrow Escape of Interacting Diffusing Particles

Physical Review Letters, 120:120601.

DOI      URL     PMID      [本文引用: 1]

The narrow escape problem deals with the calculation of the mean escape time (MET) of a Brownian particle from a bounded domain through a small hole on the domain's boundary. Here we develop a formalism which allows us to evaluate the nonescape probability of a gas of diffusing particles that may interact with each other. In some cases the nonescape probability allows us to evaluate the MET of the first particle. The formalism is based on the fluctuating hydrodynamics and the recently developed macroscopic fluctuation theory. We also uncover an unexpected connection between the narrow escape of interacting particles and thermal runaway in chemical reactors.

Allen R J, Warren P B, Ten Wolde P R. 2005.

Sampling rare switching events in biochemical networks

Physical Review Letters, 94:018104.

DOI      URL     PMID      [本文引用: 1]

Bistable biochemical switches are widely found in gene regulatory networks and signal transduction pathways. Their switching dynamics are difficult to study, however, because switching events are rare, and the systems are out of equilibrium. We present a simulation method for predicting the rate and mechanism of the flipping of these switches. We apply it to a genetic switch and find that it is highly efficient. The path ensembles for the forward and reverse processes do not coincide. The method is widely applicable to rare events and nonequilibrium processes.

Arnold V I. 1984.

Catastrophe Theory

Berlin: Heidelberg: Springer-Verlag.

[本文引用: 1]

Bandrivskyy A, Beri S, Luchinsky D G. 2003 a.

Noise-induced shift of singularities in the pattern of optimal paths

Physics Letters A, 314:386-391.

DOI      URL     [本文引用: 2]

Bandrivskyy A, Beri S, Luchinsky D G, Mannella R, McClintock P V E. 2003 b.

Fast Monte Carlo simulations and singularities in the probability distributions of nonequilibrium systems

Physical Review Letters, 90:210201.

DOI      URL     PMID      [本文引用: 2]

A numerical technique is introduced that reduces exponentially the time required for Monte Carlo simulations of nonequilibrium systems. Results for the quasistationary probability distribution in two model systems are compared with the asymptotically exact theory in the limit of extremely small noise intensity. Singularities of the nonequilibrium distributions are revealed by the simulations.

Ben-Jacob E, Bergman D J, Matkowsky B J, Schuss Z. 1982.

Lifetime of oscillatory steady-states

Physical Review A, 26:2805-2816.

DOI      URL     [本文引用: 1]

Benzi R, Parisi G, Sutera A, Vulpiani A. 1982.

Stochastic resonance in climatic change

Tellus, 34:10-16.

DOI      URL     [本文引用: 1]

Benzi R, Sutera A, Vulpiani A. 1981.

The mechanism of stochastic resonance

Journal of Physics A: Mathematical and General, 14:L453.

DOI      URL     [本文引用: 2]

Beri S, Mannella R, Luchinsky D G, Silchenko A N, McClintock P V E. 2005.

Solution of the boundary value problem for optimal escape in continuous stochastic systems and maps

Physical Review E, 72:036131.

DOI      URL     [本文引用: 9]

Beri S, Mannella R, McClintock P V E. 2004.

Dynamic importance sampling for the escape problem in nonequilibrium systems: Observation of shifts in optimal paths

Physical Review Letters, 92:020601.

DOI      URL     PMID      [本文引用: 3]

The activation problem is investigated in two-dimensional nonequilibrium systems. A numerical approach based on dynamic importance sampling (DIMS) is introduced. DIMS accelerates the simulations and allows the investigation to access noise intensities that were previously forbidden. The escape path is observed to be shifted compared to a heteroclinic trajectory calculated in the limit of zero-noise intensity. A theory to account for such shifts is presented and shown to agree with the simulations for a wide range of noise intensities.

Bernt Øksendal. 2010.

Stochastic Differential Equations.

Berlin, Heidelberg: Springer-Verlag.

[本文引用: 2]

Bobrovsky B Z, Schluss Z. 1982.

A singular perturbation method for computation of the mean first passage time in a nonlinear filter

SIAM J. Appl. Math., 42:174-187.

DOI      URL     [本文引用: 1]

Bouchet F, Laurie J, Zaboronski O. 2014.

Langevin dynamics, large deviations and instantons for the quasi-geostrophic model and two-dimensional Euler equations

Journal of Statistical Physics, 156:1066-1092.

DOI      URL     [本文引用: 1]

We investigate a class of simple models for Langevin dynamics of turbulent flows, including the one-layer quasi-geostrophic equation and the two-dimensional Euler equations. Starting from a path integral representation of the transition probability, we compute the most probable fluctuation paths from one attractor to any state within its basin of attraction. We prove that such fluctuation paths are the time reversed trajectories of the relaxation paths for a corresponding dual dynamics, which are also within the framework of quasi-geostrophic Langevin dynamics. Cases with or without detailed balance are studied. We discuss a specific example for which the stationary measure displays either a second order (continuous) or a first order (discontinuous) phase transition and a tricritical point. In situations where a first order phase transition is observed, the dynamics are bistable. Then, the transition paths between two coexisting attractors are instantons (fluctuation paths from an attractor to a saddle), which are related to the relaxation paths of the corresponding dual dynamics. For this example, we show how one can analytically determine the instantons and compute the transition probabilities for rare transitions between two attractors.

Bray A J, McKane A J. 1989.

Instanton calculation of the escape rate for activation over a potential barrier driven by colored noise

Physical Review Letters, 62:493-496.

DOI      URL     PMID      [本文引用: 2]

Bressloff P, Newby J. 2013.

Stochastic models of intracellular transport

Reviews of Modern Physics, 85:135-196.

DOI      URL     [本文引用: 1]

The interior of a living cell is a crowded, heterogenuous, fluctuating environment. Hence, a major challenge in modeling intracellular transport is to analyze stochastic processes within complex environments. Broadly speaking, there are two basic mechanisms for intracellular transport: passive diffusion and motor-driven active transport. Diffusive transport can be formulated in terms of the motion of an overdamped Brownian particle. On the other hand, active transport requires chemical energy, usually in the form of adenosine triphosphate hydrolysis, and can be direction specific, allowing biomolecules to be transported long distances; this is particularly important in neurons due to their complex geometry. In this review a wide range of analytical methods and models of intracellular transport is presented. In the case of diffusive transport, narrow escape problems, diffusion to a small target, confined and single-file diffusion, homogenization theory, and fractional diffusion are considered. In the case of active transport, Brownian ratchets, random walk models, exclusion processes, random intermittent search processes, quasi-steady-state reduction methods, and mean-field approximations are considered. Applications include receptor trafficking, axonal transport, membrane diffusion, nuclear transport, protein-DNA interactions, virus trafficking, and the self-organization of subcellular structures. DOI: 10.1103/RevModPhys.85.135

Cai R, Chen X, Duan J, Kurths J, Li X. 2017.

Lévy noise-induced escape in an excitable system

Journal of Statistical Mechanics: Theory and Experiment, 2017: 063503.

[本文引用: 2]

Cameron M K. 2012.

Finding the quasipotential for nongradient SDEs

Physica D: Nonlinear Phenomena, 241:1532-1550.

DOI      URL     [本文引用: 7]

Chan H B, Dykman M I, Stambaugh C. 2008.

Paths of fluctuation induced switching

Physical Review Letters, 100:130602.

DOI      URL     PMID      [本文引用: 1]

We demonstrate that the paths followed by a system in fluctuation-activated switching form a narrow tube in phase space. A theory of the path distribution is developed and its direct measurement is performed in a micromechanical oscillator. The experimental and theoretical results are in excellent agreement, with no adjustable parameters. We also demonstrate the lack of time-reversal symmetry in switching of systems far from thermal equilibrium.

Chatterjee M, Robert M E. 2001.

Noise enhances modulation sensitivity in cochlear implant listeners: Stochastic resonance in a prosthetic sensory system

JARO - Journal of the Association for Research in Otolaryngology, 2:159-171.

DOI      URL     PMID      [本文引用: 1]

Cochlear implants restore auditory sensitivity to the profoundly hearing-impaired by means of electrical stimulation of residual auditory nerve fibers. Sensorineural hearing loss results in a loss of spontaneous activity among the remaining auditory neurons and is accompanied by a reduction in the normal stochastic nature of neural firing in response to electric stimulation. It has been hypothesized that the natural stochasticity of the neural response is important for auditory signal processing and that introducing some optimal amount of noise into the stimulus may improve auditory perception through the implant. In this article we show that, for soft but audible stimuli, an optimal amount of

Chen L C, Deng M L, Zhu W Q. 2009.

First passage failure of quasi integrable-Hamiltonian systems under combined harmonic and white noise excitations

Acta Mechanica, 206:133-148.

DOI      URL     [本文引用: 1]

The first passage failure of multi-degree-of-freedom (MDOF) quasi integrable-Hamiltonian systems under combined harmonic and white noise excitations in the case of external resonance is studied. First, a stochastic averaging method for quasi integrable-Hamiltonian systems under combined harmonic and white noise excitations using generalized harmonic functions is reviewed briefly. Then, a backward Kolmogorov equation governing the conditional reliability function and a Pontryagin equation governing the conditional mean of the first passage time are established from the averaged It? equations, respectively. The conditional reliability function, and the conditional probability density and conditional mean of the first passage time are obtained from solving these equations together with suitable initial condition and boundary conditions. The comparison between the analytical results and those from Monte Carlo simulation for an example shows that the proposed method works very well. It is also shown by using Monte Carlo simulation that the reliability of the system in the case of external resonance is much lower than that in the non-resonant case.

Chen Z, Li Y, Liu X. 2016.

Noise induced escape from a nonhyperbolic chaotic attractor of a periodically driven nonlinear oscillator

Chaos, 26:935-992.

[本文引用: 1]

Chen L, Li Z, Zhuang Q, Zhu W. 2013.

First-passage failure of single-degree-of-freedom nonlinear oscillators with fractional derivative

Journal of Vibration and Control, 19:2154-2163.

DOI      URL     [本文引用: 1]

Chen Z, Liu X. 2016.

Patterns and singular features of extreme fluctuational paths of a periodically driven system

Physics Letters A, 380:1953-1958.

DOI      URL     [本文引用: 2]

Chen Z, Liu X. 2017 a.

Noise induced transitions and topological study of a periodically driven system

Communications in Nonlinear Science and Numerical Simulation, 48:454-461.

DOI      URL     [本文引用: 1]

Chen Z, Liu X. 2017 b.

Subtle escaping modes and subset of patterns from a nonhyperbolic chaotic attractor

Physical Review E, 95:012208.

DOI      URL     PMID      [本文引用: 4]

Noise-induced escape from the domain of attraction of a nonhyperbolic chaotic attractor in a periodically excited nonlinear oscillator is further investigated. Deviations are found to be amplified at the primary homoclinic tangency from which the optimal force begins to fluctuate dramatically. Escaping trajectories turn out to possess several modes to pass through the saddle cycle on the basin boundary, and each mode corresponds to a certain type of value of the action plot, respectively. A subset of the pattern of fluctuational paths from the chaotic attractor is obtained, showing the existence of complicated singularities.

Chen L, Zhu W Q. 2010 a.

First passage failure of quasi non-integrable generalized Hamiltonian systems

Archive of Applied Mechanics, 80:883-893.

DOI      URL     [本文引用: 1]

The first passage failure of quasi non-integrable generalized Hamiltonian systems is studied. First, the generalized Hamiltonian systems are reviewed briefly. Then, the stochastic averaging method for quasi non-integrable generalized Hamiltonian systems is applied to obtain averaged It? stochastic differential equations, from which the backward Kolmogorov equation governing the conditional reliability function and the Pontryagin equation governing the conditional mean of the first passage time are established. The conditional reliability function and the conditional mean of first passage time are obtained by solving these equations together with suitable initial condition and boundary conditions. Finally, an example of power system under Gaussian white noise excitation is worked out in detail and the analytical results are confirmed by using Monte Carlo simulation of original system.

Chen L C, Zhu W Q. 2010 b.

Reliability of quasi integrable generalized Hamiltonian systems

Probabilistic Engineering Mechanics, 25:61-66.

DOI      URL     [本文引用: 1]

Abstract

The reliability of quasi-integrable generalized Hamiltonian systems is studied. An m-dimensional integrable generalized Hamiltonian system has M Casimir functions C1,…,CM and n (n=(m−M)/2) independent first integrals HM+1,…,HM+n in involution. When an integrable generalized Hamiltonian system is subjected to light dampings and weakly stochastic excitations, it becomes a quasi-integrable generalized Hamiltonian system. The averaged Itô equations for slowly processes C1,…,CM,HM+1,…,HM+n can be obtained by using stochastic averaging method, from which a backward Kolmogorov equation governing the conditional reliability function and a Pontryagin equation governing the conditional mean of the first passage time are established. The conditional reliability function and the conditional mean of first passage time are obtained by solving these equations together with suitable initial condition and boundary conditions. Finally, an example of a 5-dimensional quasi-integrable generalized Hamiltonian system is worked out in detail and the solutions are confirmed by using a Monte Carlo simulation of the original system.

Chen L, Zhu W. 2010 c.

First passage failure of dynamical power systems under random perturbations

Science China Technological Sciences, 53:2495-2500.

DOI      URL     [本文引用: 1]

Chen Z, Zhu J, Liu X. 2017.

Crossing the quasi-threshold manifold of a noise-driven excitable system

Proceedings of the Royal Society A, 473:20170058.

[本文引用: 3]

Cohen J, Lewis R. 1967.

A ray method for the asymptotic solution of the diffusion equation

IMA J Appl Math, 3:266-290.

DOI      URL     [本文引用: 1]

Crandall M G, Evans L C, Lions P L. 1984.

Some properties of viscosity solutions of Hamilton-Jacobi equations

Transactions of the American Mathematical Society, 282:487-502.

DOI      URL     [本文引用: 1]

Crandall M G, Lions P L. 1983.

Viscosity solutions of Hamilton-Jacobi equations

Transactions of the American Mathematical Society, 277:1-42.

DOI      URL     [本文引用: 1]

Crooks G E, Chandler D. 2001.

Efficient transition path sampling for nonequilibrium stochastic dynamics

Physical Review E, 64:026109.

DOI      URL     [本文引用: 1]

Dahiya D, Cameron M K. 2018.

An ordered line integral method for computing the quasi-potential in the case of variable anisotropic diffusion

Physica D: Nonlinear Phenomena, 382-383:33-45.

DOI      URL     [本文引用: 1]

Deng M, Zhu W. 2009.

Some applications of stochastic averaging method for quasi Hamiltonian systems in physics

Science in China, Series G: Physics, Mechanics and Astronomy, 52:1213-1222.

DOI      URL     [本文引用: 1]

Dykman M I. 2010.

Poisson-noise-induced escape from a metastable state

Physical Review E, 81:051124.

DOI      URL     [本文引用: 1]

Dykman M I, Luchinsky D G, McClintock P V E, Smelyanskiy V N. 1996.

Corrals and critical behavior of the distribution of fluctuational paths

Physical Review Letters, 77:5229-5232.

DOI      URL     PMID      [本文引用: 2]

Dykman M I, McClintock P V E, Smelyanski V N, Stein N D, Stocks N G. 1992.

Optimal paths and the prehistory problem for large fluctuations in noise-driven systems

Physical Review Letters, 68:2718-2721.

DOI      URL     PMID      [本文引用: 6]

Dykman M I, Millonas M M, Smelyanskiy V N. 1994.

Observable and hidden singular features of large fluctuations in nonequilibrium systems

Physics Letters A, 195:53-58.

DOI      URL     [本文引用: 2]

Einstein A. 1905.

On the motion of small particles suspended in liquids at rest required by the molecular-kinetic theory of heat

Annalen der Physik, 322:549-560.

DOI      URL     [本文引用: 1]

Ermentrout B. 1996.

Type I membranes, phase resetting curves, and synchrony

Neural Computation, 8:979-1001.

DOI      URL     PMID      [本文引用: 1]

Type I membrane oscillators such as the Connor model (Connor et al. 1977) and the Morris-Lecar model (Morris and Lecar 1981) admit very low frequency oscillations near the critical applied current. Hansel et al. (1995) have numerically shown that synchrony is difficult to achieve with these models and that the phase resetting curve is strictly positive. We use singular perturbation methods and averaging to show that this is a general property of Type I membrane models. We show in a limited sense that so called Type II resetting occurs with models that obtain rhythmicity via a Hopf bifurcation. We also show the differences between synapses that act rapidly and those that act slowly and derive a canonical form for the phase interactions.

Faradjian A K, Elber R. 2004.

Computing time scales from reaction coordinates by milestoning

Journal of Chemical Physics, 120:10880-10889.

DOI      URL     PMID      [本文引用: 1]

An algorithm is presented to compute time scales of complex processes following predetermined milestones along a reaction coordinate. A non-Markovian hopping mechanism is assumed and constructed from underlying microscopic dynamics. General analytical analysis, a pedagogical example, and numerical solutions of the non-Markovian model are presented. No assumption is made in the theoretical derivation on the type of microscopic dynamics along the reaction coordinate. However, the detailed calculations are for Brownian dynamics in which the velocities are uncorrelated in time (but spatial memory remains).

Fenichel N. 1979.

Geometric singular perturbation theory for ordinary differential equations

Journal of Differential Equations, 31:53-98.

DOI      URL     [本文引用: 1]

Freidlin M I, Wentzell A D. 2012.

Random Perturbations of Dynamical Systems

Berlin Heidelberg: Springer.

[本文引用: 11]

Gammaitoni L, Hänggi P, Jung P, Marchesoni F. 1998.

Stochastic resonance

Reviews of Modern Physics, 70:223-287.

DOI      URL     [本文引用: 3]

Glowacki D R, Paci E, Shalashilin D V. 2009.

Boxed molecular dynamics: A simple and general technique for accelerating rare event kinetics and mapping free energy in large molecular systems

Journal of Physical Chemistry B, 113:16603-16611.

DOI      URL     [本文引用: 1]

Grafke T, Grauer R, Schäfer T. 2015.

The instanton method and its numerical implementation in fluid mechanics

Journal of Physics A: Mathematical and Theoretical, 48:333001.

DOI      URL     [本文引用: 2]

Gu X, Zhu W. 2014.

A stochastic averaging method for analyzing vibro-impact systems under Gaussian white noise excitations

Journal of Sound and Vibration, 333:2632-2642.

DOI      URL     [本文引用: 1]

A new stochastic averaging method for predicting the response of vibro-impact (VI) systems to random perturbations is proposed. First, the free VI system (without damping and random perturbation) is analyzed. The impact condition for the displacement is transformed to that for the system energy. Thus, the motion of the free VI systems is divided into periodic motion without impact and quasi-periodic motion with impact according to the level of system energy. The energy loss during each impact is found to be related to the restitution factor and the energy level before impact. Under the assumption of lightly damping and weakly random perturbation, the system energy is a slowly varying process and an averaged Ito stochastic differential equation for system energy can be derived. The drift and diffusion coefficients of the averaged Ito equation for system energy without impact are the functions of the damping and the random excitations, and those for system energy with impact are the functions of the damping, the random excitations and the impact energy loss. Finally, the averaged Fokker-Plank-Kolmogorov (FPK) equation associated with the averaged Ito equation is derived and solved to yield the stationary probability density of system energy. Numerical results for a nonlinear VI oscillator are obtained to illustrate the proposed stochastic averaging method. MonteCarlo simulation (MCS) is also conducted to show that the proposed stochastic averaging method is quite effective. (C) Elsevier Ltd.

Guardia M, Seara T M, Teixeira M A. 2011.

Generic bifurcations of low codimension of planar Filippov systems

Journal of Differential Equations, 250:1967-2023.

DOI      URL     [本文引用: 2]

In this article some qualitative and geometric aspects of non-smooth dynamical systems theory are discussed. The main aim of this article is to develop a systematic method for studying local (and global) bifurcations in non-smooth dynamical systems. Our results deal with the classification and characterization of generic codimension-2 singularities of planar Filippov Systems as well as the presentation of the bifurcation diagrams and some dynamical consequences. (c) 2010 Elsevier Inc.

Gutkin B S, Jost J, Tuckwell H C. 2009.

Inhibition of rhythmic neural spiking by noise: The occurrence of a minimum in activity with increasing noise

Naturwissenschaften, 96:1091-1097.

DOI      URL     PMID      [本文引用: 1]

The effects of noise on neuronal dynamical systems are of much current interest. Here, we investigate noise-induced changes in the rhythmic firing activity of single Hodgkin-Huxley neurons. With additive input current, there is, in the absence of noise, a critical mean value mu=mu(c) above which sustained periodic firing occurs. With initial conditions as resting values, for a range of values of the mean micro near the critical value, we have found that the firing rate is greatly reduced by noise, even of quite small amplitudes. Furthermore, the firing rate may undergo a pronounced minimum as the noise increases. This behavior has the opposite character to stochastic resonance and coherence resonance. We found that these phenomena occurred even when the initial conditions were chosen randomly or when the noise was switched on at a random time, indicating the robustness of the results. We also examined the effects of conductance-based noise on Hodgkin-Huxley neurons and obtained similar results, leading to the conclusion that the phenomena occur across a wide range of neuronal dynamical systems. Further, these phenomena will occur in diverse applications where a stable limit cycle coexists with a stable focus.

Han Q, Xu W, Yue X. 2016.

Exit location distribution in the stochastic exit problem by the generalized cell mapping method.

Chaos, Solitons and Fractals, 87:302-306.

DOI      URL     [本文引用: 1]

Han Q, Xu W, Yue X, Zhang Y. 2015.

First-passage time statistics in a bistable system subject to Poisson white noise by the generalized cell mapping method

Communications in Nonlinear Science and Numerical Simulation, 23:220-228.

DOI      URL    

Heymann M, Vanden-Eijnden E. 2008.

The geometric minimum action method: A least action principle on the space of curves

Communications on Pure and Applied Mathematics, 61:1052-1117.

DOI      URL     [本文引用: 6]

Higham D J. 2001.

An algorithmic introduction to numerical simulation of stochastic differential equations

SIAM Review, 43:525-546.

DOI      URL     [本文引用: 1]

Holcman D, Schuss Z. 2004.

Escape through a small opening: Receptor trafficking in a synaptic membrane

Journal of Statistical Physics, 117:975-1014.

DOI      URL     [本文引用: 1]

We model the motion of a receptor on the membrane surface of a synapse as free Brownian motion in a planar domain with intermittent trappings in and escapes out of corrals with narrow openings. We compute the mean confinement time of the Brownian particle in the asymptotic limit of a narrow opening and calculate the probability to exit through a given small opening, when the boundary contains more than one. Using this approach, it is possible to describe the Brownian motion of a random particle in an environment containing domains with small openings by a coarse grained diffusion process. We use the results to estimate the confinement time as a function of the parameters and also the time it takes for a diffusing receptor to be anchored at its final destination on the postsynaptic membrane, after it is inserted in the membrane. This approach provides a framework for the theoretical study of receptor trafficking on membranes. This process underlies synaptic plasticity, which relates to learning and memory. In particular, it is believed that the memory state in the brain is stored primarily in the pattern of synaptic weight values, which are controlled by neuronal activity. At a molecular level, the synaptic weight is determined by the number and properties of protein channels (receptors) on the synapse. The synaptic receptors are trafficked in and out of synapses by a diffusion process. Following their synthesis in the endoplasmic reticulum, receptors are trafficked to their postsynaptic sites on dendrites and axons. In this model the receptors are first inserted into the extrasynaptic plasma membrane and then random walk in and out of corrals through narrow openings on their way to their final destination.

Horsthemke W, Lefever R. 1984.

Noise-induced Transitions: Theory and Applications in Physics, Chemistry, and Biology

Berlin: Springer Verlag.

[本文引用: 1]

Huang Y, Liu X. 2012.

Stochastic stability of viscoelastic system under non-Gaussian colored noise excitation.

Science China: Physics, Mechanics and Astronomy, 55:483-492.

DOI      URL     [本文引用: 1]

Huang D, Yang J, Zhang J, Liu H. 2018.

An improved adaptive stochastic resonance method for improving the efficiency of bearing faults diagnosis

Journal of Mechanical Engineering Science, 232:2352-2368.

[本文引用: 1]

Ji P, Lu W, Kurths J. 2018.

Stochastic basin stability in complex networks

EPL, 122:1-6.

[本文引用: 1]

Khovanov I A, Khovanova N A, McClintock P V E. 2003.

Noise-induced failures of chaos stabilization: Large fluctuations and their control

Physical Review E, 67:051102.

DOI      URL     [本文引用: 1]

Klosek-Dygas M M, Matkowsky B J, Schuss Z. 2006.

Stochastic stability of nonlinear oscillators

SIAM Journal on Applied Mathematics, 48:1115-1127.

DOI      URL     [本文引用: 2]

Knessl C, Matkowsky B, Schuss Z, Tier C. 1985.

An asymptotic theory of large deviations for Markov jump processes

SIAM Journal on Applied Mathematics, 45:1006-1028.

DOI      URL    

Kong C, Gao X, Liu X. 2016.

On the global analysis of a piecewise linear system that is excited by a Gaussian white noise

Journal of Computational and Nonlinear Dynamics, 11:051029.

[本文引用: 2]

Kong C, Liu X. 2017.

Noise-induced chaos in a piecewise linear system

International Journal of Bifurcation and Chaos, 27:1750137.

DOI      URL     [本文引用: 4]

Kougioumtzoglou I A, Spanos P D. 2013.

Response and first-passage statistics of nonlinear oscillators via a numerical path integral approach

Journal of Engineering Mechanics, 139:1207-1217.

DOI      URL     [本文引用: 1]

A numerical path integral solution approach is developed for determining the response and first-passage probability density functions (PDFs) of nonlinear oscillators subject to evolutionary broad-band stochastic excitations. Specifically, based on the concepts of statistical linearization and of stochastic averaging, the system response amplitude is modeled as a one-dimensional Markov diffusion process. Further, using a discrete version of the Chapman-Kolmogorov equation and the associated first-order stochastic differential equation, the response amplitude and first-passage PDFs are derived. The main concept of the approach relates to the evolution of the response PDF in short time steps, assuming a Gaussian form for the conditional response PDF. A number of nonlinear oscillators are considered in the numerical examples section including the versatile Preisach hysteretic oscillator. For this oscillator, first-passage PDFs are derived for the first time to the authors' best knowledge. Comparisons with pertinent Monte Carlo data demonstrate the reliability of the approach.

Kramers H A. 1940.

Brownian motion in a field of force and the diffusion model of chemical reactions

Physica, 7:284-304.

DOI      URL     [本文引用: 2]

Kraut S, Feudel U. 2003 a.

Enhancement of noise-induced escape through the existence of a chaotic saddle

Physical Review E, 67:015204.

DOI      URL     [本文引用: 1]

Kraut S, Feudel U. 2003 b.

Noise-induced escape through a chaotic saddle: lowering of the activation energy

Physica D: Nonlinear Phenomena, 181:222-234.

DOI      URL     [本文引用: 1]

Kubo R. 1966.

The fluctuation-dissipation theorem

Reports on Progress in Physics, 29:255-284.

DOI      URL     [本文引用: 1]

Kuznetsov Y A, Rinaldi S, Gragnani A. 2003.

One-parameter bifurcations in planar Filippov systems

International Journal of Bifurcation and Chaos, 13:2157-2188.

DOI      URL     [本文引用: 3]

Lee DeVille R E, Vanden-Eijnden E, Muratov C B. 2005.

Two distinct mechanisms of coherence in randomly perturbed dynamical systems

Physical Review E, 72:31105.

DOI      URL     [本文引用: 2]

Li W, Chen L, Trisovic N, Cvekovic A, Zhao J. 2015.

First passage of stochastic fractional derivative systems with power-form restoring force.

International Journal of Non-Linear Mechanics, 71:83-88.

DOI      URL     [本文引用: 1]

Lin Y K, Cai G Q. 1995.

Probabilistic Structural Dynamics: Advanced Theory and Applications.

Boston: McGraw-Hill.

[本文引用: 1]

Longtin A. 1997.

Autonomous stochastic resonance in bursting neurons

Physical Review E, 55:868-876.

DOI      URL     [本文引用: 1]

Lu S, He Q, Zhang H, Kong F. 2017.

Rotating machine fault diagnosis through enhanced stochastic resonance by full-wave signal construction

Mechanical Systems and Signal Processing, 85:82-97.

DOI      URL     [本文引用: 1]

Luchinsky D G, Beri S, Mannella R, McClintock P V E, Khovanov I A. 2002.

Optimal fluctuations and the control of chaos

International Journal of Bifurcation and Chaos, 12:583-604.

DOI      URL     [本文引用: 1]

Luchinsky D G, Maier R S, Mannella R, McClintock P V E, Stein D L. 1997.

Experiments on critical phenomena in a noisy exit problem

Physical Review Letters, 79:3109-3112.

DOI      URL     [本文引用: 3]

Luchinsky D G, Maier R S, Mannella R, McClintock P V E, Stein D L. 1999.

Observation of saddle-point avoidance in noise-induced escape

Physical Review Letters, 82:1806.

DOI      URL     [本文引用: 4]

Luchinsky D G, McClintock P V E. 1997.

Irreversibility of classical fluctuations studied in analogue electrical circuits

Nature, 389:463-466.

DOI      URL    

Luchinsky D G, McClintock P V E, Dykman M I. 1998.

Analogue studies of nonlinear systems

Reports on Progress in Physics, 61:889-997.

DOI      URL     [本文引用: 1]

Lücken L, Yanchuk S, Popovych O V, Tass P A. 2013.

Desynchronization boost by non-uniform coordinated reset stimulation in ensembles of pulse-coupled neurons

Frontiers in Computational Neuroscience, 7:63.

DOI      URL     PMID      [本文引用: 1]

Several brain diseases are characterized by abnormal neuronal synchronization. Desynchronization of abnormal neural synchrony is theoretically compelling because of the complex dynamical mechanisms involved. We here present a novel type of coordinated reset (CR) stimulation. CR means to deliver phase resetting stimuli at different neuronal sub-populations sequentially, i.e., at times equidistantly distributed in a stimulation cycle. This uniform timing pattern seems to be intuitive and actually applies to the neural network models used for the study of CR so far. CR resets the population to an unstable cluster state from where it passes through a desynchronized transient, eventually resynchronizing if left unperturbed. In contrast, we show that the optimal stimulation times are non-uniform. Using the model of weakly pulse-coupled neurons with phase response curves, we provide an approach that enables to determine optimal stimulation timing patterns that substantially maximize the desynchronized transient time following the application of CR stimulation. This approach includes an optimization search for clusters in a low-dimensional pulse coupled map. As a consequence, model-specific non-uniformly spaced cluster states cause considerably longer desynchronization transients. Intriguingly, such a desynchronization boost with non-uniform CR stimulation can already be achieved by only slight modifications of the uniform CR timing pattern. Our results suggest that the non-uniformness of the stimulation times can be a medically valuable parameter in the calibration procedure for CR stimulation, where the latter has successfully been used in clinical and pre-clinical studies for the treatment of Parkinson's disease and tinnitus.

Ludwig D. 1975.

Persistence of dynamical systems under random perturbations

SIAM Review, 17:605-640.

DOI      URL     [本文引用: 5]

Maier R S, Stein D L. 1993 a.

Effect of focusing and caustics on exit phenomena in systems lacking detailed balance

Physical Review Letters, 71:1783-1786.

DOI      URL     PMID      [本文引用: 7]

Maier R S, Stein D L. 1993 b.

Escape problem for irreversible systems

Physical Review E, 48:931-938.

DOI      URL     [本文引用: 2]

Maier R S, Stein D L. 1996.

A scaling theory of bifurcations in the symmetric weak-noise escape problem

Journal of Statistical Physics, 83:291-357.

DOI      URL     [本文引用: 4]

Maier R S, Stein D L. 1997.

Limiting exit location distributions in the stochastic exit problem

SIAM Journal on Applied Mathematics, 57:752-790.

DOI      URL     [本文引用: 3]

Marshall J S. 2016.

Analytical solutions for an escape problem in a disc with an arbitrary distribution of exit holes along its boundary

Journal of Statistical Physics, 165:920-952.

DOI      URL     [本文引用: 1]

Matkowsky B J, Schuss Z. 1977.

The exit problem for randomly perturbed dynamical systems

SIAM Journal on Applied Mathematics, 33:365-382.

DOI      URL     [本文引用: 1]

Matkowsky B, Schuss Z. 1982.

Diffusion across characteristic boundaries

SIAM Journal on Applied Mathematics, 42:822-834.

DOI      URL     [本文引用: 2]

Matkowsky B J, Schuss Z, Ben-Jacob E. 1982.

A singular perturbation approach to Kramers' duffusion problem

SIAM Journal on Applied Mathematics, 42:835-849.

DOI      URL     [本文引用: 1]

Matkowsky B J, Schuss Z, Tier C. 1984.

Uniform expansion of the transition rate in Kramers' problem

Journal of Statistical Physics, 35:443-456.

DOI      URL     [本文引用: 1]

Menck P J, Heitzig J, Kurths J, Schellnhuber H J. 2014.

How dead ends undermine power grid stability

Nature Communications, 5:3969.

DOI      URL     PMID      [本文引用: 1]

The cheapest and thus widespread way to add new generators to a high-voltage power grid is by a simple tree-like connection scheme. However, it is not entirely clear how such locally cost-minimizing connection schemes affect overall system performance, in particular the stability against blackouts. Here we investigate how local patterns in the network topology influence a power grid's ability to withstand blackout-prone large perturbations. Employing basin stability, a nonlinear concept, we find in numerical simulations of artificially generated power grids that tree-like connection schemes--so-called dead ends and dead trees--strongly diminish stability. A case study of the Northern European power system confirms this result and demonstrates that the inverse is also true: repairing dead ends by addition of a few transmission lines substantially enhances stability. This may indicate a topological design principle for future power grids: avoid dead ends.

Menck P J, Heitzig J, Marwan N, Kurths J. 2013.

How basin stability complements the linear-stability paradigm

Nature Physics, 9:89-92.

URL     [本文引用: 1]

The human brain(1,2), power grids(3), arrays of coupled lasers(4) and the Amazon rainforest(5,6) are all characterized by multistability(7). The likelihood that these systems will remain in the most desirable of their many stable states depends on their stability against significant perturbations, particularly in a state space populated by undesirable states. Here we claim that the traditional linearization-based approach to stability is too local to adequately assess how stable a state is. Instead, we quantify it in terms of basin stability, a new measure related to the volume of the basin of attraction. Basin stability is non-local, nonlinear and easily applicable, even to high-dimensional systems. It provides a long-sought-after explanation for the surprisingly regular topologies(8-10) of neural networks and power grids, which have eluded theoretical description based solely on linear stability(11-13). We anticipate that basin stability will provide a powerful tool for complex systems studies, including the assessment of multistable climatic tipping elements(14).

Muratov C B, Vanden-Eijnden E, E W. 2005.

Self-induced stochastic resonance in excitable systems

Physica D: Nonlinear Phenomena, 210:227-240.

DOI      URL     [本文引用: 5]

Naeh T, Klosek M M, Matkowsky B J, Schuss Z. 1990.

A direct approach to the exit problem

SIAM Journal on Applied Mathematics, 50:595-627.

DOI      URL     [本文引用: 3]

Newby J M, Bressloff P C, Keener J P. 2013.

Breakdown of fast-slow analysis in an excitable system with channel noise

Physical Review Letters, 111:128101.

DOI      URL     PMID      [本文引用: 3]

We consider a stochastic version of an excitable system based on the Morris-Lecar model of a neuron, in which the noise originates from stochastic sodium and potassium ion channels opening and closing. One can analyze neural excitability in the deterministic model by using a separation of time scales involving a fast voltage variable and a slow recovery variable, which represents the fraction of open potassium channels. In the stochastic setting, spontaneous excitation is initiated by ion channel noise. If the recovery variable is constant during initiation, the spontaneous activity rate can be calculated using Kramer's rate theory. The validity of this assumption in the stochastic model is examined using a systematic perturbation analysis. We find that, in most physically relevant cases, this assumption breaks down, requiring an alternative to Kramer's theory for excitable systems with one deterministic fixed point. We also show that an exit time problem can be formulated in an excitable system by considering maximum likelihood trajectories of the stochastic process.

Nolting B C, Abbott K C. 2016.

Balls, cups, and quasi-potentials: Quantifying stability in stochastic systems

Ecology, 97:850-864.

DOI      URL     PMID      [本文引用: 1]

When a system has more than one stable state, how can the stability of these states be compared? This deceptively simple question has important consequences for ecosystems, because systems with alternative stable states can undergo dramatic regime shifts. The probability, frequency, duration, and dynamics of these shifts will all depend on the relative stability of the stable states. Unfortunately, the concept of

Pei B, Xu Y, Yin G. 2017.

Stochastic averaging for a class of two-time-scale systems of stochastic partial differential equations.

Nonlinear Analysis, Theory, Methods and Applications, 160:159-176.

DOI      URL     [本文引用: 1]

Pikovsky A S, Kurths J. 1997.

Coherence resonance in a noise-driven excitable system

Physical Review Letters, 78:775-778.

DOI      URL     [本文引用: 2]

Qiao Z, Lei Y, Lin J, Jia F. 2017.

An adaptive unsaturated bistable stochastic resonance method and its application in mechanical fault diagnosis

Mechanical Systems and Signal Processing, 84:731-746.

DOI      URL     [本文引用: 1]

Reimann P, Schmid G J, H鋘ggi P. 1999.

Universal equivalence of mean first-passage time and Kramers rate

Physical Review E, 60:R1-R4.

DOI      URL     [本文引用: 1]

Rodrigo G, Stocks N G. 2018.

Suprathreshold stochastic resonance behind cancer

Trends in Biochemical Sciences, 43:483-485.

DOI      URL     PMID      [本文引用: 1]

Noise in gene expression is pervasive and, in some cases, even fulfills a functional role. Cancer cell populations exploit noise to increase heterogeneity as a defense against therapies. What lies behind this picture is a phenomenon of stochastic resonance led by the collective, rather than by individual cells.

Roy R V. 1993.

Noise perturbations of nonlinear dynamical systems

Computational Stochastic Mechanics, 79:125-148.

[本文引用: 1]

Roy R V. 1994 a.

Asymptotic analysis of first-passage problems.

International Journal of Non-Linear Mechanics, 32:173-186.

DOI      URL     [本文引用: 2]

Roy R V. 1994 b.

Noise perturbations of a non-linear system with multiple steady states.

International Journal of Non-Linear Mechanics, 29:755-773.

DOI      URL     [本文引用: 1]

Roy R V. 1995.

Noise-induced transitions in weakly nonlinear oscillators near resonance

Journal of Applied Mechanics, 62:496-504.

DOI      URL     [本文引用: 1]

Schuecker J, Diesmann M, Helias M. 2015.

Modulated escape from a metastable state driven by colored noise

Physical Review E, 92:052119.

DOI      URL     [本文引用: 1]

Schultz P, Menck P J, Heitzig J, Kurths J. 2017.

Potentials and limits to basin stability estimation

New Journal of Physics, 19:023005.

DOI      URL     [本文引用: 1]

Schuss Z, Matkowsky B. 1979.

The exit problem: A new approach to diffusion across potential barriers

SIAM Journal on Applied Mathematics, 36:604-623.

DOI      URL     [本文引用: 1]

Schuss Z, Singer A, Holcman D. 2007.

The narrow escape problem for diffusion in cellular microdomains

Proceedings of the National Academy of Sciences, 104:16098-16103.

DOI      URL     [本文引用: 1]

Schuss Z, Spivak A. 1998.

Where is the exit point?

Chemical Physics, 235:227-242.

DOI      URL     [本文引用: 2]

Sethian J A, Vladimirsky A. 2001.

Ordered upwind methods for static Hamilton-Jacobi equations

Proceedings of the National Academy of Sciences, 98:11069-11074.

DOI      URL     [本文引用: 1]

Sethian J A, Vladimirsky A. 2003.

Ordered upwind methods for static Hamilton-Jacobi equations: Theory and algorithms

SIAM Journal on Numerical Analysis, 41:325-363.

DOI      URL     [本文引用: 1]

Sidney Redner. 2001. A Guide to First-Passage Processes. Cambridge: Cambridge University Press.

[本文引用: 1]

Smelyanskiy V N, Dykman M I. 1997.

Optimal control of large fluctuations

Physical Review E, 55:2516-2521.

DOI      URL     [本文引用: 3]

Smelyanskiy V N, Dykman M I, Maier R S. 1997.

Topological features of large fluctuations to the interior of a limit cycle

Physical Review E, 55:2369-2391.

DOI      URL     [本文引用: 1]

Spivak A, Schuss Z. 2002 a.

Analytical and numerical study of Kramers' exit problem I.

Applied Mathematics E-Notes, 2:132-140.

Spivak A, Schuss Z. 2002 b.

The exit distribution on the stochastic separatrix in Kramers' exit problem

SIAM Journal on Applied Mathematics, 62:1698-1711.

DOI      URL     [本文引用: 1]

Spivak A, Schuss Z. 2003.

Analytical and numerical study of Kramers' exit problem II.

Applied Mathematics E-Notes, 3:147-155.

[本文引用: 1]

Stocks N G, Allingham D, Morse R P. 2002.

The application of suprathreshold stochastic resonance to cochlear implant coding

Fluctuation and Noise Letters, 02:L169-L181.

DOI      URL     [本文引用: 1]

Sun J Q, Hsu C S. 1988.

First-passage time probability of non-linear stochastic systems by generalized cell mapping method

Journal of Sound and Vibration, 124:233-248.

DOI      URL     [本文引用: 1]

Tél T, Lai Y C. 2010.

Quasipotential approach to critical scaling in noise-induced chaos

Physical Review E, 81:56208.

DOI      URL     [本文引用: 1]

Tél T, Lai Y C, Gruiz M. 2008.

Noise-induced chaos: A consequence of long deterministic transients

International Journal of Bifurcation and Chaos, 18:509-520.

DOI      URL     [本文引用: 1]

Touchette H. 2009.

The large deviation approach to statistical mechanics

Physics Reports, 478:1-69.

DOI      URL     [本文引用: 1]

Tuckwell H C, Jost J. 2012.

Analysis of inverse stochastic resonance and the long-term firing of Hodgkin-Huxley neurons with Gaussian white noise

Physica A: Statistical Mechanics and Its Applications, 391:5311-5325.

DOI      URL     [本文引用: 1]

Wang J. 2015.

Landscape and flux theory of non-equilibrium dynamical systems with application to biology

Advances in Physics, 64:1-137.

DOI      URL     [本文引用: 1]

Wang W, Yan Z, Liu X. 2017.

The escape problem and stochastic resonance in a bistable system driven by fractional Gaussian noise

Physics Letters A, 381:2324-2336.

DOI      URL     [本文引用: 1]

Weber J. 1956.

Fluctuation dissipation theorem

Physical Review, 101:1620-1626.

DOI      URL     [本文引用: 1]

Wechselberger M, Mitry J, Rinzel J.

2013. Canard theory and excitability

Lecture Notes in Mathematics, 2102: 89-132.

[本文引用: 1]

Wentzell A D, Freidlin M I. 1970.

On small random perturbations of dynamical systems

Russian Mathematical Surveys, 25:1-55.

[本文引用: 2]

Whitney H. 1955.

On singularities of mappings of Euclidean spaces. I. Mappings of the plane into the plane

Annals of Mathematics, 62:374-410.

DOI      URL     [本文引用: 1]

Wuehr M, Boerner J C, Pradhan C , et al. 2017.

Stochastic resonance in the human vestibular system - Noise-induced facilitation of vestibulospinal reflexes

Brain Stimulation, 11:261-263.

DOI      URL     PMID      [本文引用: 1]

BACKGROUND: There is strong evidence that the presence of noise can enhance information processing in sensory systems via stochastic resonance (SR). OBJECTIVES: To examine the presence of SR in human vestibulospinal reflex function. METHODS: Healthy subjects were stimulated with 1 Hz sinusoidal GVS of varying amplitudes (0-1.9 mA). Coherence between GVS input and stimulation-induced motion responses was determined and psychometric function fits were subsequently used to determine individual vestibulospinal reflex thresholds. This procedure was repeated with additional application of imperceptible white noise GVS (nGVS). RESULTS: nGVS significantly facilitated the detectability of weak subthreshold vestibular inputs (p < 0.001) and thereby effectively lowered the vestibulospinal threshold in 90% of participants (p < 0.001, mean reduction: 17.5 +/- 14.6%). CONCLUSION: This finding provides evidence for the presence of SR-dynamics in the human vestibular system and gives a functional explanation for previously observed ameliorating effects of low-intensity vestibular noise stimulation on balance control in healthy subjects and patients with vestibular hypofunction.

Xu M. 2018.

First-passage failure of linear oscillator with non-classical inelastic impact

Applied Mathematical Modelling, 54:284-297.

DOI      URL     [本文引用: 1]

Xu Y, Guo R, Liu D, Zhang H, Duan J. 2013.

Stochastic averaging principle for dynamical systems with fractional brownian motion

Discrete and Continuous Dynamical Systems - Series B (DCDS-B), 19:1197-1212.

[本文引用: 1]

Yang S, Potter S F, Cameron M K. 2019.

Computing the quasipotential for nongradient SDEs in 3D

Journal of Computational Physics, 379:325-350.

DOI      URL     [本文引用: 1]

Young H P. 2015.

The evolution of social norms

Annual Review of Economics, 7:359-387.

DOI      URL     [本文引用: 1]

Zhu J, Chen Z, Liu X. 2018.

Probability evolution method for exit location distribution

Physics Letters A, 382:771-775.

Zhu W Q, Huang Z L, Deng M L. 2002.

Feedback minimization of first-passage failure of quasi non-integrable Hamiltonian systems.

International Journal of Non-Linear Mechanics, 37:1057-1071.

DOI      URL     [本文引用: 1]

Abstract

The non-linear stochastic optimal control of quasi non-integrable Hamiltonian systems for minimizing their first-passage failure is investigated. A controlled quasi non-integrable Hamiltonian system is reduced to an one-dimensional controlled diffusion process of averaged Hamiltonian by using the stochastic averaging method for quasi non-integrable Hamiltonian systems. The dynamical programming equations and their associated boundary and final time conditions for the problems of maximization of reliability and of maximization of mean first-passage time are formulated. The optimal control law is derived from the dynamical programming equations and the control constraints. The dynamical programming equations for maximum reliability problem and for maximum mean first-passage time problem are finalized and their relationships to the backward Kolmogorov equation for the reliability function and the Pontryagin equation for mean first-passage time, respectively, are pointed out. The boundary condition at zero Hamiltonian is discussed. Two examples are worked out to illustrate the application and effectiveness of the proposed procedure.

Zhu W Q, Wu Y J. 2003.

First-passage time of Duffing oscillator under combined harmonic and white-noise excitations

Nonlinear Dynamics, 32:291-305.

DOI      URL     [本文引用: 1]

The first-passage time of Duffing oscillator under combined harmonic andwhite-noise excitations is studied. The equation of motion of the system is firstreduced to a set of averaged Itô stochastic differential equations by using thestochastic averaging method. Then, a backward Kolmogorov equation governing theconditional reliability function and a set of generalized Pontryagin equationsgoverning the conditional moments of first-passage time are established. Finally, theconditional reliability function, and the conditional probability density and momentsof first-passage time are obtained by solving the backward Kolmogorov equation andgeneralized Pontryagin equations with suitable initial and boundary conditions.Numerical results for two resonant cases with several sets of parameter values areobtained and the analytical results are verified by using those from digital simulation.

/

〈 〉

Baidu
map