Bar模型双稳态螺旋波研究.docx

Br模型双稳态螺旋波研究 高 见, 王 群，吕华平 （江苏师范大学 物理与电子工程学院，江苏 徐州 221100） 摘要：螺旋波是自然界中最为常见的时空斑图结构，在诸多生化反应以及数值模拟中都可以观察到，但是对于双稳态介质中螺旋波的时空动力学行为的研究却鲜有提及。本文通过Br模型模拟研究了螺旋波的动力学行为，并解释了螺旋波的形成机理；建立了双稳态介质存在螺旋波所要满足的方程；同时预测了在特定参数下会产生成稳定螺旋波、线状和带状漫游螺旋波；得出了振荡和可激发介质都是双稳态介质的极限状态，并且通过数值模拟进行了验证。关键词：螺旋波; Br模型; 双稳态介质 The research of the br model in bistable state Gao Jian（School of Physics & Electrical Engineering, Jiangsu Normal University, Xuzhou 221100, Jiangsu, China）Abstract：Spiral wave, which is the mostcommon, is the time and space pattern structure. Could be observed in a variety of biochemical reaction and numerical modeling. Yet, the research of spatio-temporal dynamics of the br model in bistable state is scarcely mentioned. This paper has researched the dynamics of spiral wave, got the formation mechanism, built the equation that describe formation condition of spiral wave in this bistable state，predicted the existence of stable spiral waves, threadiness and zonary meandering spiral waves, and got the conclusion that oscillation media and excitable media are the limiting condition of bistable state.Key words： Spiral wave；br model；bistable state螺旋波在自然界和实验中广泛存在，是远离平衡态系统中最为常见的耗散结构之一，各类螺旋波结构一直都是斑图动力学研究的热点1-3。在很多生化反应扩散系统的实验过程中，包括CO在Pt表面的催化氧化过程4-7，BZ反应8, 9过程和糖酵解10-12等过程中都发现了螺旋波。产生螺旋波的介质可分为可激发介质，振荡介质和双稳态介质三种。目前的研究都集中于可激发介质，简单振荡介质以及个别倍周期振荡周期2系统中的螺旋波。已经在实验和数值模拟中发现了诸如波尖漫游而导致的超螺旋波13-15，碎片螺旋波16-19，以及倍周期振荡系统所对应的线缺陷螺旋波20-22等，但是对于双稳态介质中的螺旋波的报道鲜见报道。双稳系统的主要特征是有两个稳定的定态并且都有各自的吸引域。在相空间中，当相点的轨道运动到某个稳态的吸引域范围内时，就会被吸引到该稳态上，如果没有受到足够大的干扰，相点会永远停留在该稳态上；当扰动足够大时，相点会被激发出该稳态，进入另一个稳态的吸引域，最终停留在另一个稳态上。如果系统中有足够强的连续扰动，相点会在这两个稳态上振荡23。木文通过对一个二变量模型描述的具有双稳态特性的反应散系统中所呈现的螺旋波形式的研究，观测了螺旋波的动力学规律，解释了螺旋波的形成机理，建立了双稳态介质存在螺旋波所要满足的方程，同时预测了在特定参数下会产生成稳定螺旋波、线状漫游和面状漫游螺旋波，得出了双稳态介质在极限状态下会过渡到振荡介质与可激发介质。1 模型介绍Bar 模型由Bar 和Eiswirthr于1993年提出，该模型的方程为： ut=1u1-uu-v+ba+2u (1a) vt=fu-v (1b) 其中 fu=0 0u1 (1c) 方程含有3个参数：，a和b。参数是一个比较小的数,揭示了u和v之间的时间尺度关系。由于通常情况下比较小,因此,u是快变量,v是慢变量。a和b的取值决定系统的类型，当ab+1而b0时,系统是振荡介质,而当a0或a b +1且b b +1而b 0时,则是双稳介质。 当系统是可激发介质时，刻画了系统可激发性的强弱。若较小，则系统的可激发性较强，随着的增大，系统的可激发性逐渐下降以致波不能在该系统中传播。不同螺旋波的状态将不一样。下图为系统在不同的参数（，b）下的相图。图1.系统在切空间中的状态6Fig.1 The status of the system in tangent space6图1画出了在激发介质区域系统在不同参数大小下的相图，其中 N（no wave）是无波区间，F（flat wave）是平面波区域，S（spiral with rigid rotation）是周期螺旋波区域，M（meandering spiral）是漫游螺旋波区域，T1 是由多普勒不稳定性导致的混沌区域，T2 是由回火（backfiring）所导致的混沌区域，两条虚线分别是sl（saddle loop bifurcation）鞍节点分岔和hopf（Hopf bifurcation）霍普夫分岔。鞍节点分岔产生了一个环绕不稳定固定点的极限环，而霍普夫分岔使得这个不稳定的固定点变为稳定的固定点。2系统存在螺旋波解的理论及数值模拟本文下面的数值模拟中，系统的参数a和b固定在ab+1，b0的区域内，即系统控制在双稳态领域内，算法采用欧拉法，时间步长取为t=0.02，空间差分采用五点式，步长取h=1，采用零流边界条件，系统尺寸为LxLy=256256。2.1介质中行波的形成条件与机理去掉方程（1）中的扩散项，方程变为： ut=1u1-uu-v+ba (2a) vt=fu-v (2b) 其中 fu=0 0u0.068时，其反应项为正，且在一定范围内随着u值的增加而增加。当u1/3时，v开始增加，最终达到稳态u=1,v=1（高稳态）。于此同时，该点相邻的处于高稳态的点，由于扩散作用而逐渐减小，当低于闸值0.97时，便会迅速的降低到低稳态。为了简化问题，令方程（1）中的扩散项2u=2u2x，方程（1a）变为： ut=1u1-uu-v+ba+2u2x (3) 我们利用含有100个点的一位数列进行数值模拟，前50个点的初值取u=1,v=1,后50个点的初值取u=0,v=0, 取a=1.1,b=0.07和=0.05，空间步长取h=1，时间步长取为t=0.02，采用无流边界条件。数值模拟的结果如图3（a）所示，在两种稳态相邻处，高稳态的点，被激发到低稳态；低稳态的点被激发到高稳态。最终在低稳态区域激发出一个波峰，并且向两边扩展。u值能否被激发，由它对应的v值决定，因为v值决定了u的闸值。如图2（b）所示，随着v值的增加，高稳态的u值被激发的闸值就越小；低稳态的u值被激发的闸值就越大。图4描述了相邻的处于两种稳态的点，通过相互激发，在相空间中的运动情况。为了验证该猜想，我们在方程7的扩散项前加上扩上系数Du，通过不断减小扩散系数而减弱系统中的扩散作用。 ut=1u1-uu-v+ba+Du2u2x (4) 图2（b）描述了u的反应项大小在相空间中的分布；（c）描述了当v=0和v=1时u的反应项大小与u的关系。我们定义u的反应项大小为： fm=utm=1u1-uu-v+ba (5) dfmdu=Au2+Bu+C （6） A=-3 B=21+v+ba （7） C=-v+ba 令dfmdu=0，得 fm极值点横坐标： u=-BB2-4AC2A （8）低稳态u值的闸值为： ug=-B+B2-4AC2A （9）u的扩散项大小为： fp= ujtp=Du2uj2x=Duuj+1-2uj+uj-1h2 (10) 首先求低稳态的点被激发到高稳态的条件。取参数a=1.1,b=0.07,=0.05，u的闸值为u=0.031296，反应项为 fm=-0.019609。当扩散项 fp与反应项相同时，对应的扩散系数为Du=0.020918。而高稳态的点被激发到低稳态所对应的闸值、反应项和扩散系数分别为u=0.986458， fm=-0.00668，Du=0.003771。按照上文所述的原理，理论上随着扩散系数的减小，系统的可激发性会减弱，最终系统将停止反应而稳定于这两种稳态。为了更清晰模拟点被激发的临界条件，我们用含6个点的数组进行模拟。首先令前三个点固定于高稳态，后三个点处于低稳态。如图3（b）所示，取Du=0.0209时，第4个点的u值增长到u=0.031296,然后停止增长，此时扩散项和反应项达到平衡；当把扩散系数稍微增加，取Du=0.0210时，便能够把第四个点激发起来。然后令前三个点处于高稳态，后三个点处于低稳态。如图11（c）所示，取Du=0.003771时，第3个点的u值降低到u=0.0.9865, 然后停止降低，此时扩散项和反应项达到平衡；当把扩散系数稍微增加，取Du=0.003772时，该点依然没有被激发，当Du=0.00388时，同样没有完成激发。原因在于，与低稳态附近v值的行为相同，如图10（b）所示，在高稳态附近v值的反应项为负，导致在该点u值的反应项和扩散项竞争时，而v值在不停的减小，直到减小到v对应的零线上。由方程9可知，随着v值的减小，u值的反应项会增加，故扩散作用的阻力变大，最终达到平衡。如果令第3个点的v值固定于v=1时，当Du=0.003772时，则可以完成激发。在临界情况下，u值的反应项和扩散项作用的合力效果与v值的反应项大小相仿。考虑到临界情况下高稳态被激发的特殊性：相点不是沿着直线v=1缓慢变化，而是沿着零线变化。把u和v的关系代入方程9得： fm=utm=1u1-uu-1-6.75uu-12+ba (11) dfmdu=1 -56.75au4+126.75au3-31+36.75au2+21+6.75+1+bau-1a+ba 12令(df_m)/du=0，得高稳态的点的v值不受控制时能够被激发到低稳态所对应的闸值、反应项和扩散系数分别为u=0.984226， fm=0.004037，Du=0.004169。利用上述数值模拟方法，取Du=0.004169时，第3个点的u值降低到u=0.984226, 然后停止降低，此时扩散项和反应项达到平衡，如图3（d）所示；当把扩散系数稍微增加，取Du=0.004170时，该点被激发到低稳态。图3. 参数与图2中相同；（a）：当Du=1时，两种稳态交界处的动力学过程；(b)：处于低稳态的u值被激发的临界状态；（c）：处于高稳态（v固定于1）的u值被激发的临界状态；（d）：处于高稳态的u值被激发的临界状态；Fig.3 Parameters are the same in Fig.2;(a):The dynamics of the junction when Du=1;(b):The critical state of u could be inspired staying at low stable state;(c)&(d): The critical state of u could be inspired staying at upper stable state上述中，在模拟低稳态的点被激发的临界条件时，数组中的前3个点是固定在高稳态的，原因是高稳态的点被激发的临界条件远远小于低稳态的临界条件。所以第4个点被激发到高稳态附近后又被激发到低稳态，导致该点做上下振荡运动。如图4（a）所示，该点的u值从0缓慢的增加；当达到临界点时，其增加的速度趋于0；当突破临界点后，u值迅速增加；达到高稳态附近后，又被后快的激发到低稳态。图4（b）更形象的描述了该点在相空间中的运动。图4.（a）：初态处于低稳态的点被反复激发，在原点振荡，图为对应u值的振荡状态；（b）：该点在相空间中的轨线Fig.4 (a):The point being inspired repeatedly is oscillating;(b): The path in phase space 数值模拟与理论相吻合，证实了猜想的正确性。明确了激发的机理，便可以建立系统能够支撑螺旋波的条件。理论和实验证明，在双稳态系统中振荡的周期可分为3部分：静息态过程用时、激发过程用时和动力学过程用时。当激发比较容易时，激发过程用时可以忽略，周期主要由动力学用时组成；当激发难度较大时，激发用时便不能忽略。图5中，红线为u值从低值反应到高值的过程，其用时记为Tlow；绿线为u值从高值反应到高低值的过程，用时为Thigh。图5.稳态被激发后，在相空间中的轨线；把轨线分为4段，并取近似函数Fig.5 The path in phase space of the point staying at steady state be inspired is divided into four segments and find the approximate function结合动力学方程求其值： Tlow=0.550-vdv+011u(1-u)u-f2u+ba+ Du2 du (13) Thigh=0.051dv1-v+101u(1-u)u-f1u+ba+ Du2du (14)支持螺旋波的条件为： 0Tlsfpu+fmudt=umin (15a) fpumin+fmumin=0 (15b) 0Thsfpu+fmudt=umax (15c) fpumax+fmumax=0 (15d) fpumin=-2umin (15e) fpumax=2-2umax (15f) TlsThigh (15g) ThsTlow (15h)方程15a表示低稳态被自由激发的条件，经过时间Tls，在扩散项和反应项的共同作用下，当u值增加到umin，且TlsThigh时，处于低稳态的点才能被激发到高稳态。方程15c表示高稳态被自由激发的条件，经过时间Ths，在扩散项和反应项的共同作用下，当u值下降到umax，且ThsTlow时，处于高稳态的点才能被激发到低稳态。由方程5、15b和15e可以求出umin；同理，由方程5、15c和15f可以求出umax。 2.2 螺旋波波头的形成条件和机理图6.(a)：螺旋波形成的过程，参数与图2中相同，第一列分别是u和v的初始条件，在该初始条件下会形成两个方向相反的行波 ，两者交织在一起，最终形成螺旋波；(b)：黑色区域是螺旋波波头在相空间中的值域，红线和绿线是方程的零线，黑色区域里的交点是稳定焦点Fig.6 (a): The forming process of spiral wave,parameters are the same as in Fig.2;(b):The black area is the range of the spiral wave tip 如图6（b）所示，稳定焦点周围的行为比较复杂，该点把红色区域和蓝色区域的交界线分为两段，下半段交界线附近红色区域的点，由于其对应的v值比较小，故相应的u的闸值比较大，不容易被激发到低稳态附近；而下半段交界线附近蓝色区域的点，同样由于其对应的v值比较小，相应的u的闸值比较小，所以容易被激发到高稳态附近。同理，上半段交界线附近红色区域的点，对应的v值比较大，u容易被激发到低稳态附近；而蓝色区域的点则不容易被激发到高稳态附近。所以该焦点便成为划分两段行波的点。以该点为坐标原点作极坐标，由于行波上个点的行进速度相同，故距离原点越远的点的角速度越小，最终形成螺旋波。但是，螺旋波的形成与该焦点无关，螺旋波波头与之也没必然联系，如图7(a)(b)所示。图7.(a)：改变参数b，使稳定焦点变为切点，以及使之相离，其它参数与图2中相同；（b）：当b取不同值时，所形成的螺旋波；(c)：稳定的螺旋波；（d）：图（c）中波头的放大图；（e）：图（d）中波头处对应的慢变量分布；（f）：黑色区域脱离同步反应 Fig.7 (a)&(b): Zero curves and spiral waves in different parameter b;(c):The stable spiral wave;(d):Spiral wave tip of the spiral wave in (c);(e):Recovery variable near spiral tip;(f):Special area of the spiral wave现在研究螺旋波波头附近的行为。在稳定的螺旋波中，如图7（c）所示，该螺旋波波头稳定于一定点。（d）是波头的放大图，图中黑点为波头所在点。可见，波头两侧可看做两行方向相反的行波，在波头附近，行波的前进会受到慢变量的抑制。如图（e），是图 （d）中慢变量的大小分布，距离波头越近，慢变量的延迟效应越明显，以至于阻止了波头邻域中的点参与反应。如图（f），黑色区域脱离了同步反应。但是该区域会受到外部扩散作用的影响，除中心点以外，其它点的快变量依然做周期变化，由于动力学项，也伴随着慢变量的周期运动，在相空间中的轨线集中在图6（b）的灰色区域。2.3 螺旋波波头漫游的过程本文已经描述了波头周围的动力学行为，如图9所示。但是一般情况下高稳态和低稳态的激发难易程度会有一定的差异，就会导致波头漫游。下面我们讨论波头漫游的机理。 如图15所示，描述了波头产生及漫游的过程，与图7（b）和图9中描述的螺旋波不同，该系统中，高稳态容易被激发到低稳态，而低稳态较难激发，所以波头两侧行波的波速不同，如图（a）和（b）所示，波头右侧的行波会传播到左侧行波的后面，由于慢变量的延迟性而造成最终波速相等。如图（c）所示，图中圈出的波已经不能再继续向右传播，达到相互激发的临界点，由于该区域下方附近的波在传播，激发性较强，造成该区域下方的波开始反向传播，全部反向 ，在此过程中波头不断向上移动，如图（d）所示。图8.螺旋波形成及漫游的过程，a=1.32,b=0.30,=0.05,空间步长dx=0.05,时间步长dt=0.0005；（c）：椭圆圈出的行波突然改变原来的传播方向，反向传播；（d）：漫游后的波头；漫游呈跳跃式的迁移方式Fig.8 The formation and meandering of spiral wave, a=1.32,b=0.30,=0.05,space grid dx=0.05,time step dt=0.0005;(c):The marked traveling wave with oval reverses direction suddenly;(d):The spiral wave after meandering 3 螺旋波结构预测及数值模拟方程是否有螺旋波解，是由扩散项和动力学项的竞争结果而决定的。那么，当系统满足螺旋波解的条件，且低稳态和高稳态的激发过程用时都远小于动力学过程用时的条件下，系统支撑稳定的螺旋波；当高稳态或低稳态的激发过程用时与动力学过程用时相近时，由于两种稳态的激发难易程度不同，螺旋波波头会漫游。当低稳态较难激发时，那么高稳态被激发到低稳态的速度大于低稳态被激发的速度，导致波峰很窄，形成线状螺旋波，波峰宽度由慢变量v决定；当高稳态较难激发时，相反，波谷则很窄，形成带状螺旋波，宽度同样由慢变量v决定。而当高稳态或低稳态的激发过程用时其中一个远小于动力学过程用时，即Ths0或Tls0时，系统演变为可激发介质；当高稳态和低稳态的激发过程用时都远小于动力学过程用时，即Tls0和Ths0时，系统演变为振荡介质。振荡介质和可激发介质都是双稳态介质的极限情况，所以方程15依然是振荡介质和可激发介质支撑螺旋波的条件。可见，振荡介质中必然有螺旋波解。数值模拟证实了所预测的现象。如图9所示，描述了螺旋波的3种状态：稳定螺旋波、线状漫游螺旋波和面状漫游螺旋波。图9（a）是当低稳态和高稳态的激发过程用时都远小于动力学过程用时所产生的稳定螺旋波；（b）是当低稳态的激发过程用时与动力学过程用时相近时产生的线状漫游螺旋波；（c）是当高稳态的激发过程用时与动力学过程用时相近时产生的带状漫游螺旋波。图9.当=0.05时，a, b 分别取不同值时，系统所支撑的螺旋波；（a）：当a=1.064, b=0.053时，产生稳定的螺旋波；（b）：当a=1.316, b=0.303时，产生线状漫游螺旋波；（c）：当a=1.149, b=0.057时，产生带状漫游螺旋波Fig.9 =0.05，(a):We got this stable spiral wave when a=1.064, b=0.053;(b): got this wirelike meandering spiral wave when a=1.316, b=0.303;(c): got this zonary meandering spiral wave when a=1.149, b=0.057图10.当=0.05时，（a）：a=0.769, b=0.154，振荡介质中产生稳定的螺旋波；（b）：当a=1.205, b=0.265时,可激发介质中产生线状漫游螺旋波；（c）：当a=1.005, b=-0.051时，可激发介质中产生带状漫游螺旋波Fig.10 =0.05，(a): a=0.769, b=0.154 ,we got this stable spiral wave in oscillation media;(b): a=1.205, b=0.265 ,got this wirelike meandering spiral wave;(c): a=1.005,b=-0.051 , got this zonary meandering spiral wave图10（a）是振荡介质产生的稳定螺旋波；图10（b）和（c）是两种激发介质中产生的螺旋波，它们分别与图9中（b）和（c）十分相似。充分证明了上述推理。4 结论动力学方程1中的反应项和扩散项的竞争结果决定了该系统能否支持螺旋波。通过改变高稳态与低稳态的闸值，改变激发难易程度，可以得到不同类型的螺旋波。本文建立了双稳态介质存在螺旋波所要满足的方程，解释了波头漫游、线状漫游螺旋波和面状漫游螺旋波的产生机理。推导及证明了振荡介质和可激发介质都是双稳态介质的极限情况。目前还没有得出详细描述漫游螺旋波波头运动的方程，这正是下一步研究工作的内容。参考文献1I.R.Epstein. An introduction to nonlinear chemical dynamics: oscillations, waves, patterns, and chaosM. New York: Oxford University Press 1998.2I.R.Epstein, K.Showalter. Nonlinear chemical dynamics: oscillations, patterns, and chaosJ. Phys Chem 1996.3Kapral R, K.Showalter. Chemical waves and patternsM. Kluwer Academic Publishers,Boston 1995.4G.Ertl, Mikhailov AS. Nonequilibrium structures in condensed systemsJ. Science 1996:1596-1597.5G.Ertl, Rotermund H-H. Spatiotemporal pattern formation in reaaction at surfacesJ. Curr Opinion in Solid State and Mat Sci 1996:5.6Br M, Eiswirth M. Turbulence due to spiral breakup in a continuous excitable mediumJ. Phys Rev E 1993:3.7杨灵法, 中怀 侯, 文 辛厚. 一氧化碳在Pt_110_面催化氧化过程中化学振荡和化学混沌J. 化学物理学报 1999;12:4.8Belousov BP. A periodic reaction and its mechannism, in collection of short papers on radiation medicne for 1958J. Med Publ 1959.9Belousov BP, Field RJ, Buger M. Oscillations and Traveling Waves in Chemical SystemsM. New York: Wiley 1985.10 P.Ball. The self-made tapestry:pattern formation in natureM. New York: Oxford University Press 1999.11Westermark PO, Lansner A. A Model of Phosphofructokinase and Glycolytic Oscillations in the Pancreatic b-cellJ. Biophysical Journal 2003;85:126139.12Wierschem K, Bertram R. Complex bursting in pancreatic islets: a potential glycolytic mechanismJ. Journal of theoretical biology 2004;228:513-521.13 Barkley D. Enclideam sy