科学杂志文章-图灵斑图动力学(欧阳颀)

斑图动力学
斑图动力学是研究复杂系统中各种不同类型的斑点或模式如何形成、演变和消失的一门学科。

它主要用于研究多相流、流体力学、化学反应和生物学等领域。

斑图动力学是一个非常复杂的领域，研究的对象包括物理、化学和生物等多个领域。

斑图动力学的研究方法主要包括数学建模、计算机模拟和实验研究。

斑图动力学可以帮助我们理解复杂系统中各种不同类型的斑点或模式的形成机制、演化规律和消失机制等。

在物理学领域，斑图动力学主要用于研究多相流中的模式形成。

例如，在流体力学中，斑点可能是由不同密度、温度或速度的流体相互作用而形成的。

在化学领域中，斑点可能是由不同化学物质在反应过程中形成的。

在生物学领域，斑点可能是由生物体在生长过程中形成的。

斑图动力学的研究还可以帮助我们解释复杂系统中的许多现象，例如生物多样性、环境污染、气候变化等。

斑图动力学是一个跨学科领域，它涉及到数学、物理、化学、生物学等多个学科。

该论文发表在《中国物理快报》（Chinese Physics letters）18卷第6期834-836页。

作者简介：毛晓明，北京大学物理学院物理系九八级本科生。

从2000年4月开始在北京大学物理系介观与人工微结构国家重点实验室下属的非线性科学及生物技术实验室工作。

本文发表于Chinese Physics letters，为第二作者。

有两篇论文正在准备中，分别是美国UIC大学廖秀北教授指导的DNA与蛋白质相互作用几何特性与识别模式研究的论文，和与本实验室孙锴同学合作的一篇关于社会演化的论文。

感悟与寄语：从儿时开始的梦想，终于在这一刻逐渐变成现实。

在心中神圣殿堂面前，我们的每一天，都感到幸福和激动，握有无穷的力量。

我们不畏惧前路的荆棘，因为希望就在前方。

吾生也有涯，而知也无涯，我们正在准备着把自己的人生投入到人类探索自然的洪流中去。

指导教师简介：欧阳颀教授是北京大学物理系教授、长江学者。

1989年于法国波尔多第一大学获博士学位,后一直从事非线性科学的基础理论与实验研究。

主要研究方向是非线性动力学中的斑图自组织行为。

近十年来在该领域取得了一系列重大成果，被国际同行公认为斑图动力学领域的实验科学带头人之一。

迄今为止在各类科学杂志上共发表论文近四十篇，其中包括《自然》杂志三篇，《科学》杂志两篇，《物理通讯快报》三篇；八次应邀在不同国际会议上作专题报告，其中包括美国物理学会，加拿大化学学会，美国工业与应用数学学会；十余次在各个大学和研究单位作专题讲座。

1996年受聘于日本电器公司（NEC）在美国的研究中心，从事生物计算机的研究开发和其他一些生物基因工程问题。

1997年在《科学》上发表了DNA计算机的论文。

文章引起了国际新闻界的广泛注意。

包括英国《新科学家》，美国《新闻周刊》，日本教育电视台在内的数家新闻媒介作了报道。

1998年6月到北京大学物理系从事非线性科学与生物芯片技术开发工作。

Turing 不稳定性及斑图形成摘要:在这篇文中，我们借助于浮游植物-浮游动物的数学模型来研究Turing 不稳定是如何产生的.首先介绍了Turing 不稳定产生的内在机理，给出了详细的过程，并且最终得出了产生Turing 不稳定的参数空间.然后在结合含有扩散项的浮游植物、浮游动物的捕食模型来研究该模型是否能够产生Turing 不稳定现象. 关键词：Turing 不稳定，捕食模型1.Turing 不稳定性1952年Turing 在文中《The chemical basis of morphogenesis 》一文中提出：如果参加相互反应的化学物质自身不存在扩散作用，经过一段时间反应后，它们会达到一定的平衡状态，即这些化学物质的浓度将会变得均匀. 但如果这些化学物质具有扩散作用的话，那么在某种条件下，这种均匀的平衡态将会被打破，变成不均匀的平衡态，这边是Turing 不稳定现象. 换句话说在同一个正常数平衡解处的常微风模型是稳定的，但对于加入扩散作用的偏微分方程模型却是不稳定的.本文借助于数学模型来说明发生Turing 不稳定性的条件. 海洋中存在着多种浮游植物和浮游动物，它们的关系非常的复杂，这里我们仅分别考虑一种浮游植物、一种浮游动物，并且这种浮游动物主要以这种浮游植物为食. 浮游植物会产生毒素，可以杀死一定量的浮游动物，进而来保护自己免受捕食.并且还考虑两种浮游生物在二维平面上的空间分布，从而引入其含有Laplacian 算子的扩散项。

Spatiotemporal dynamics toxic-phytoplankton-zooplankton model ：1P P aPZ rP t K P m Z bPZ cPZ dZ t P m P m∂⎛⎫=-- ⎪∂+⎝⎭∂=--∂++（1） 这里的参数均为正常数，其中()()=,,,,P P x y t Q x y t =，分别是能够产生毒素的浮游植物、浮游动物在t 时刻(),x y 处的密度，并且浮游植物产生的毒素可以杀死浮游动物且满足第二类功能性反应函数. 浮游植物服从Logistic 的增长方式，r为其内禀增长率，K 为其环境容纳量. 浮游动物捕食浮游植物满足第二类功能性反应函数，a 为捕食率，m 为半饱和常数. b 为浮游动物捕食浮游植物转化为自身增长的效率，d 为浮游动物的死亡率，c 为浮游植物产生毒素杀死浮游动物的概率，显然要满足b c >.对于模型（1）的各个平衡点处的稳定性在文献[1]中已经研究，这里不再详细介绍，仅仅在下面简单分析其正平衡态存在、稳定的条件. 下面我们在模型（1）的基础上，考虑其扩散项，从而得到如下的模型.Spatiotemporal dynamics in a reaction-diffusion toxic-phytoplankton-zooplankton model ：()()11221,,,P P aPZ rP D P f P Z D P t K P m Z bPZ cPZ dZ D Z g P Z D Z t P m P m∂⎛⎫=--+∆+∆ ⎪∂+⎝⎭∂=--+∆+∆∂++@@（2） 且满足非零的初始条件()()(),,00,,,00,[0,][0,]P x y Q x y x y Lx Lx >>∈Ω=⨯ 以及零边界条件()0,P Q x y n n ∂∂==∈∂Ω∂∂其中,Lx Ly 分别是模型（1）在,x y 方向上的一段，向量n 是边界∂Ω上的单位外法向量，零边界条件也就说明了这个系统没有外部的输入，此时可以认为模型是独立的. 12,D D 分别表示浮游植物和浮游动物的扩散系数. ∆为二维空间上拉布拉斯算子.2222=x y∂∂∆+∂∂ 本文研究的是Turing 不稳定性，所以只需关心正平衡态，从模型（1）可以计算出本系统存在唯一的一个正平衡态为()***,E P Z =其中：()()()()**2,rm b c K b c d md md P Z b c d aK b c d ----==----并且满足：()/0K md b c d >-->.模型（1）在正平衡点*E 处的线性化模型为：P P J t Z Z ⎛⎫⎛⎫∂= ⎪ ⎪ ⎪ ⎪∂⎝⎭⎝⎭其中**,P P P Q Q Q =-=--，矩阵J 为()()()()()()()111221220rd K b c d m b c d ad J J K b c b c d b c J J J r K b c d md aK ⎛⎫----+- ⎪----⎛⎫ ⎪= ⎪ ⎪⎝⎭--- ⎪ ⎪⎝⎭@ 则由二维系统的Routh-Hurwitz 判据[1]可得正平衡点稳定的冲要条件为()()()()11221221det 0rd K b c d md J J J J J K b c ---=-=>-（3） ()()()()()()11220rd K b c d m b c d tr J J J K b c b c d ----+=+=<---（4）联合（3）、（4）式可解出参数范围为：()m b c d md K b c d b c d-+<<----（5） 接下來研究Turing 不稳定性，即是由于扩散系统引起的不稳定性. 因此，我们总假设条件（3）、（4）成立，也即式（5）式是恒成立的. 下面考虑含有扩散的模型（2），做与上述相同的平移变换,并把新的变量,P Q 仍记为,P Q ，这里的,P Q 表示模型（2）在平衡点*E 附近的扰动. 可得：1112121222P J P J Q D P t Q J P J Q D Q t ∂⎧=++∆⎪⎪∂⎨∂⎪=++∆⎪∂⎩（6） 又因为模型（6）的任意解都可以展开成下述的Fourier 级数：()()()()()(),0,0,0,0,,sin ,,cos ij ij i j i j ij ij i j i j P r t u r t t krQ r t v r t t kr αα∞∞==∞∞==⎧==⎪⎪⎨⎪==⎪⎩∑∑∑∑（7）这里向量(),r x y =，且0,0x L xy L y <<<<. 向量(),i j k k k =，且/,/i j k i Lx k j Ly ππ==，,i j k k 称为波数.把（7）式带入（6）式可得：()()211112221222ij ij ij ij ij ij J D k J t J J D k tααββαβ∂⎧=-+⎪⎪∂⎨∂⎪=+-⎪∂⎩（8） 其中222i j k k k =+，这里是因为()()()222sin sin sin sin sin i j i i j j i j kr k x k y k k x k y k k x k y k kr∆=∆+=-+-+=-模型（8）是一个常系数微分方程组，其解的形式为1212t t c e c e λλ+，其中12,c c 为常数，是由初始条件所确定. 12,λλ是其系数矩阵1J 的特征值2111121221222J D k J J J J D k ⎡⎤-=⎢⎥-⎣⎦（9） 求得此系数矩阵的行列式、迹分别为：()()()4211211222111221221det J D D k J D J D k J J J J =-++-（10）()()()21121122tr J D D k J J =-+++（11）为了研究是由于扩散发生的不稳定，系数矩阵的特征值12,λλ至少有一个是具有正实部，也就说条件()()11det 0,0J tr J ><至少有一个不成立. 有假设条件（3）、（4）恒成立，可知11220J J +<恒成立，所以得到()10tr J <恒成立，所以要使Turing 不稳定发生，存在一个参数空间使得()1det 0J <成立.令()()()2421211222111221221G k D D k J D J D k J J J J =-++-（12）这是一个关于未知数2k 的一元二次函数，由条件（3）知112212210J J J J ->成立.显然，()20G k <在()20,k ∈+∞成立的必要条件是1122210J D J D +>（13）在条件（13）成立的前提下，要使()20G k <成立，即要求()20G k =有两个实根，则必须满足系数判别式：()()211222112112212214J D J D D D J J J J +>-（14） 在满足条件（13），（14），函数（12）将会存在两个正的实根22,k k ，当满足 222k k k <<（15）时，有()20G k <，即模型（8）的系数矩阵的特征值12,λλ至少有一个是具有正实部，则模型（2）的平衡点*E 是不稳定的，此时平衡点*E 的不稳定性是由于扩散项∆算子的特征值也成波数的k 所引起的，所以称（15）式为Turing 不稳定空间. 得到Turing 不稳定的参数空间后，可以选取输入参数空间的各个参数，使得模型在这些参数下发生Turing 不稳定，进而会形成各种斑图，对于具体形成斑图，这里不做介绍.综上，可得发生Turing 不稳定性的充分必要条件是：式子（3）、（4）、（13）、（14），也即：()()()()11221221112211222121122211211221221det 0004J J J J J tr J J J J D J D J D J D D D J J J J ⎧=->⎪=+<⎪⎨+>⎪⎪+>-⎩ 对于模型数学模型（2），根据上面的Truing 不稳定的充要条件来求其Truing 不稳定的参数空间. 前面已经得到求解其雅克比矩阵J ，其中220J =，由（4）式可知()110tr J J =<，再由（13）式可得1120J D >，而这里的20D >，从而可得对于模型（2）来说，不满足上述的条件，所以并不会发生Turing 不稳定现象.通过（4）、（13）式可知，1122,J J 必须是异号的，并且负值的绝对值要大于正值的绝对值，在模型（2）中，220J =，所以其不会发生Turing 不稳定现象.2. 总结这是最近看到的一篇关于反应扩散微分方程的文章，原文中也是简单介绍个各个理论，我有利用生物数学课堂上学过的知识，进行整理。

不同图灵模作用的几种斑图白占国;董丽芳【摘要】Mechanisms of pattern formation and pattern selection with different Turing modes interaction are investigated by using a two-layer coupled CIMA model. It is shown that hexagonal superlattice and simple hexagon arise respectively in subsysteml and subsystem 2 under the condition that two subsystems locate at supercritical or subcritical bifurcation point. Both of them in two subsystems cannot interact when the two Turing modes are supercritical and one simple stripe pattern in each of sub-systems emerges spontaneously. The identical 'bean' patterns is selected in the two subsystems when two Turing modes are subcriticl. In addition, the bifurcation types of the Turing modes also affect the spatial symmetry of the e-merging patterns in system.%采用双层耦合的CIMA模型,研究了不同图灵模相互作用时斑图的选择、形成机制.结果表明:当2个子系统分别处在超临界和次临界分岔点附近时,超临界图灵模和次临界图灵模相互作用产生耦合,得到六边形和超六边斑图；当2个子系统激发的图灵模均为超临界模时,二者之间不发生耦合,每个子系统各自形成简单的条纹斑图；当2个子系统激发的图灵模均为次临界模时,2个模产生相互作用,系统最终选择完全相同的“豆角”斑图.此外,图灵模的分岔类型还改变斑图的空间对称性.【期刊名称】《河北大学学报（自然科学版）》【年(卷),期】2012(032)002【总页数】4页(P140-143)【关键词】图灵模;超点阵;超临界和次临界【作者】白占国;董丽芳【作者单位】河北大学物理科学与技术学院,河北保定071002;河北大学物理科学与技术学院,河北保定071002【正文语种】中文【中图分类】O461.2斑图（pattern）是一种典型的非线性自组织现象［1－2］，广泛地存在于自然界，也可以在不同的实验室系统中进行研究.其研究内容涉及物理学、数学、化学、生态学等各个学科，而且在心脏病的防治、材料处理和局域生长以及等离子体光子晶体等方面具有广阔的应用前景，近年来引起人们极大的兴趣.国内外学者在实验［3－7］、尤其是理论上［8－14］做了大量的研究，得到了种类丰富的斑图.例如，杨灵法等人［8－10］研究了系统处在超临界和次临界分岔时的超点阵斑图和叠加斑图的形成机理和空间共振条件.发现三波共振和相同的对称性对超点阵斑图的形成起着重要作用.Fineberg等人［11－12］研究了次临界图灵模与一个流体表面的超临界法拉第波之间的相互作用，结果表明该系统出现的自组织四边形和多种超点阵斑图都是由多个非线性的波矢构成，而且不同空间模峰值同时发生是三波共振的必要条件.Bachir小组［13］和Epstein小组［14］分别得到零模与不同空间模的相互作用，及超临界模与MASK模相互作用时的图灵斑图.从以往的研究看，前人工作大多集中研究超临界模和次临界模的相互作用，对2个图灵模均为超临界或次临界研究较少.为了进一步理解斑图形成和选择的物理机制，推进非线性科学的发展、加快其实际应用的进程，研究系统处于不同分岔点时图灵模之间的相互作用尤为重要.本工作针对此现状，采用双层耦合的CIMA模型，细致研究2个子系统激发3种分岔类型的图灵模相互作用时斑图的形成机理.采用双层耦合的CIMA模型［14］，在无量纲条件下方程可写成如下表达式：其中，i）为系统的局部动力学.式中u和v分别表示变量活化子与禁阻子，Du和D v为二变量的扩散系数，a和b是系统的控制参数，此方程组存在均匀定态解通过作线性稳定性分析得到：当控制参数b＞b H＝3a25－a125时，系统出现霍普分岔，当的条件下，系统处于图灵空间，该模型包含2个子系统：系统1（u1，v1）和系统2（u2，v2），α和β为2个子系统之间活化子和禁阻子的耦合强度，本文固定控制参量a＝15，b＝9，并选取合适的参数使2个子系统均在图灵空间，研究2个图灵模的相互作用.其他条件选择格点数为128×128，时间步长和空间步长分别为0.01和1单位进行数值模拟.图1是超临界与次临界分岔点耦合系统的色散关系及出现的斑图.从系统的色散关系（如图1a所示）可以看出，子系统1处于超临界分岔点，激发超临界图灵模k1，振幅较大为基模又叫主动模，占主导地位；子系统2则处于次临界分岔点，产生1个次临界图灵模k2，振幅较小，是次谐振模，又称从动模，二者具有不同空间尺度.超临界图灵模k1是不稳定的，在次临界图灵模k2作用下，激发出1个新的图灵模k3，三者满足三波共振关系k1＋k2＝k3，使2子系统之间发生非线性共振.长波模调制短波模，在子系统1出现超六边斑图（如图1b所示），子系统2仍然呈现简单的大点六边形斑图，如图1c所示.观察超六边的傅里叶谱发现，超点阵能量的空间分布较大点六边形更为复杂，包含3个不同空间尺度的波矢.通过调节控制参量，使得k2模由次临界的稳定模穿过虚轴变为超临界的不稳定模，子系统2经历一个非平衡相变，原来稳定的不动点变为不稳定的焦点，这种不稳定的焦点叫动力学系统的“排斥子”［1］.2个“排斥子”互相排斥，导致2个子系统不发生耦合，如图2所示.2个超临界模k1和k2之间没有相互作用，每个子系统各自出现简单的条纹斑图，而非六边形斑图.继续调节系统的控制参量，使2个不稳定的超临界模变为稳定的次临界模，考察其相互作用时对系统斑图的影响，如图3所示.从图3a可以看出，此时，2个子系统的不动点均是稳定的，都是“吸引子”.2个图灵模地位相当，二者互相竞争，相互影响，2个子系统之间出现较强的耦合现象，并且相互调制，最终出现完全相同的“豆角”斑图，即当系统处于次临界／次临界分岔点时图灵模的作用与前面2种情况均不相同.比较图1、图2和图3发现，2个子系统是否能够耦合及耦合强度的大小，敏感依赖其图灵模分岔类型：当系统处于超临界与次临界分岔点时，由于耦合作用，子系统1出现超六边斑图，子系统2则不受耦合因素影响，仍然呈现简单斑图；当2子系统激发的图灵模均为超临界模时，层与层之间不发生耦合；如果2子系统激发的图灵模均为次临界模时，2个模地位相当，二者互相竞争，相互影响，使得2个子系统强烈耦合.此外，斑图的空间对称性也因图灵模分岔类型的不同而改变，其中当系统处于超临界分岔点时系统选择的斑图空间对称性最低，是条纹对称性；次临界分岔点次之，为类四边形对称，当2个子系统分别处于超临界与次临界分岔点时，系统形成的斑图空间对称性最高，具有六边形对称性.通过对双层耦合的反应扩散方程进行线性稳定性分析，得到系统的分岔条件，选择不同的控制参数，使2个子系统分别处于超临界和次临界、超临界和超临界、次临界和次临界3种不同分岔点，在随机的初始条件下模拟了斑图选择和时空演化.模拟结果表明，不同图灵模的相互作用在斑图的选择和形成过程中起着重要作用.当2个子系统激发的图灵模分别为超临界模和次临界模时，长波模调制短波模，二者相互作用使2个子系统发生耦合，满足三波共振条件，子系统1出现超六边斑图，同时子系统2出现简单六边.2个子系统是否能够耦合，敏感依赖其图灵模分岔类型：如果2子系统激发的图灵模均为超临界模时，层与层之间不发生耦合，每层内部满足空间共振条件，每个子系统各自形成不同尺度空间的简单条纹斑图；值得注意的是，当2个图灵模均为次临界模时，2个模地位相当，其相互作用不同于前面2种情况，二者互相竞争，相互影响，使得2个子系统最终选择完全相同的“豆角”斑图模式.此外，还发现，图灵模的分岔类型还改变斑图的空间对称性.本结果对深刻理解斑图形成和选择的物理机制，推动斑图动力学的发展具有一定意义.【相关文献】［1］欧阳颀.非线性科学与斑图动力学导论［M］.北京：北京大学出版社，2010：130－131. OUYANG Q.Nonlinear science and instroduction of pattern dynamics［M］.Beijing：Peking University Press，2010：130－131.［2］CROSS M C，HOHENBERG P C.Pattern formation outside of equilibrium［J］.Rev Mod Phys，1993，65（3）：851－1086.［3］董丽芳，谢伟霞，赵海涛，等.氩气／空气介质阻挡放电中超六边斑图［J］.物理学报，2009，58：4806－4811.DONG Lifang，XIE Weixia，ZHAO Haitao，et al.Experimental study on self－organized hexagonal superlattice pattern in dielectric barrier discharge in argon／air［J］.Acta Phys Sin，2009，58：4806－4811.［4］董丽芳，赵海涛，谢伟霞，等.介质阻挡放电中四边形超晶格斑图的实验研究［J］.物理学报，2008，57：5768－5772.DONG Lifang，ZHAO Haitao，XIE Weixia，et al.Experimental investigation of square superlattice pattern formation in a dielectric barrier discharge［J］.Acta Phys Sin，2008，57：5768－5772.［5］贺亚峰，董丽芳，尹增谦，等.介质阻挡放电中斑图的傅里叶分析［J］.河北大学学报：自然科学版，2003，23（2）：137－140.HE Yafeng，DONG Lifang，YIN Zengqian，et al.Fourier analysis of patterns in dielectric barrier dischage［J］.Journal of Hebei University：Natural Science Edition，2003，23（2）：137－140.［6］宋倩，董丽芳，李媛媛，等.超六边形斑图的4种形成途径［J］.河北大学学报：自然科学版，2010，30（4）：371－374.SONG Qian，DONG Lifang，LI Yuanyuan，et al.Four pathway to formed hexagonal supperlattice［J］.Journal of Hebei University：Natural Science Edition，2010，30（4）：371－374.［7］李媛媛，董丽芳，宋倩，等.超点阵斑图形成前放电丝时空特征［J］.河北大学学报：自然科学版，2010，30（6）：643－646.LI Yuanyuan，DONG Lifang，SONG Qian，et al.Spatial and temporal characteristic of filaments before formed patterns［J］.Journal of Hebei University：Natural Science Edition，2010，30（6）：643－646.［8］YANG L F，DOLNIK M，ZHABOTINSKY A M，et al.Turing patterns beyond hexagons and stripes［J］.Chaos，2006，16：037114.［9］YANG L F，DOLNIK M，ZH ABOTINSKY A M，et al.Spatial resonances and superposition patterns in a reaction－diffusion model with interaction Turing modes ［J］.Phys Rev Lett，2002，88：208303.［10］BERENSTEIN I，YANG L F，DOLNIK M，et al.Dynamic mechanism of photochemical induction of Turing superlattices in the chlorine dioxide－iodine－malonic acid reaction－diffusion system［J］.J Phys Chem A，2005，109：5382－5387.［11］EPSTEIN T，FINEBERG J.Necessary conditions for mode interactions in parametrically excited waves［J］.Phys Rev Lett，2008，100：134101.［12］ARBELL H，FINEBERG J.Pattern formation in two－frequency forced parametric waves［J］.Phys Rev E，2002，65：036224.［13］BACHIR M，METENS S，BORCKMANS P，et al.Formation of rhombic and superlattice patterns in bistable systems［J］.Europhys Lett，2001，54：612－618. ［14］LENGYEL I，EPSTEIN I R.Modeling of Turing structures in the chlorite－iodide－malonic－acid－starch reaction system［J］.Science，1991，251：650－652.。

反应扩散系统中的图灵斑图动力学介绍：———非线性科学
专题之十
欧阳颀
【期刊名称】《物理通报》
【年(卷),期】1999(000)005
【摘要】1 图灵斑图 1952年,被后人称为计算机科学之父的著名英国数学家阿兰·图灵(Alan Turing)把他的目光转向生物学领域。

在他的一篇著名论文《形态形成的化学基础》中,图灵用一个反应扩散模型成功地说明了某些生物体表面所显示的图纹,如斑马身上的斑图是怎样产生的。

想象在生物胚胎发育的某个阶段,生物体内某些被称为'形态子'的生物大分子与其他反应物发生生物化学反应,同时在机体内随机扩散。

图灵表明在适当的条件下,这些原来浓度均匀分布的'形态子'会在空间自发的组织成一些周期性的结构。

【总页数】4页(P4-7)
【作者】欧阳颀
【作者单位】北京大学非线性科学中心
【正文语种】中文
【中图分类】O414
【相关文献】
1.图灵斑图动力学的数学机制 [J], 刘迎东
2.具有非线性发生率的传染病模型图灵斑图研究 [J], 王涛;靳祯;孙桂全
3.斑图动力学：非线性科学专题之九 [J], 夏蒙棼
4.螺旋波的斑图动力学：非线性科学专题之十一 [J], 欧阳颀
5.双层耦合非对称反应扩散系统中的振荡图灵斑图 [J], 刘雅慧;董梦菲;刘富成;田淼;王硕;范伟丽
因版权原因，仅展示原文概要，查看原文内容请购买。

物理系学术报告Physics Department Colloquium 12月18日，周五，16:00-17:00，教12-201Reverse engineering of biologicalnetworksProfessor Qi Ouyang (欧阳颀)北京大学理论生物学中心北京大学－加州大学定量生物医学联合研究中心教育部长江学者特聘教授ABSTRACTThe interdependence of topology, dynamics, and function of a biological control network is one of central questions in systems biology. Biological control networks have to perform their functions robustly. A robust function may have preferences in the topological structures of the underlying network. In this talk, we will discuss three questions on the reverse engineering of biological networks: (1) what kind of network can perform a giving biological function? (2) How many of them can perform it robustly?(3) How Nature selects those networks? We will use the segmentation polarity network in fruit fly and the cell cycle control network in budding yeast as examples to discuss those questions.Our results show that the functional and dynamical constrains largely limit the selection of biological networks.欧阳颀教授简介：北京大学物理学院长江学者特聘教授，北京大学理论生物学中心副主任，北京大学－加州大学定量生物医学联合研究中心副主任。

反应扩散模型在图灵斑图中的应用及数值模拟张荣培;王震;王语;韩子健【摘要】反应扩散方程模型常被用于描述生物学中斑图的形成.从反应扩散模型出发,理论推导得到Gierer-Meinhardt模型的斑图形成机理,解释了非线性常微分方程系统的稳定常数平衡态在加入扩散项后会发生失稳并产生图灵斑图的过程.通过计算该模型,得到图灵斑图产生的参数条件.数值方法中采用一类有效的高精度数值格式,即在空间离散条件下采用Chebyshev谱配置方法,在时间离散条件下采用紧致隐积分因子方法.该方法结合了谱方法和紧致隐积分因子方法的优点,具有精度高、稳定性好、存储量小等优点.数值模拟表明,在其他条件一定的情况下,系统控制参数κ 取不同值对于斑图的产生具有重要的影响,数值结果验证了理论结果.%Turing proposed a model for the development of patterns found in nature in 1952. Turing instability is known as diffusion-driven instability, which states that a stable spatially homogeneous equilibrium may lose its stability dueto the unequal spatial diffusion coefficients. The Gierer–Mainhardt modelis an activator and inhibitor system to model the generating mechanism of biological patterns. The reaction-diffusion system is often used to describe the pattern formation model arising in biology. In this paper, the mechanism of the pattern formation of the Gierer-Meinhardt model is deduced from the reactive diffusion model. It is explained that the steady equilibrium state of the nonlinear ordinary differential equation system will be unstable after adding of the diffusion term and produce the Turing pattern. The parameters of the Turing pattern are obtained by calculating the model. There are a variety of numerical methods including finitedifference method and finite element method. Compared with the finite difference method and finite element method, which have low order precision, the spectral method can achieve the convergence of the exponential order with only a small number of nodes and the discretization of the suitable orthogonal polynomials. In the present work, an efficient high-precision numerical scheme is used in the numerical simulation of the reaction-diffusion equations. In spatial discretization, we construct Chebyshev differentiation matrices based on the Chebyshev points and use these matrices to differentiate the second derivative in the reaction-diffusion equation. After the spatial discretization, we obtain the nonlinear ordinary differential equations. Since the spectral differential matrix obtained by the spectral collocation method is full and cannot use the fast solution of algebraic linear equations, we choose the compact implicit integration factor method to solve the nonlinear ordinary differential equations. By introducing a compact representation for the spectral differential matrix, the compact implicit integration factor method uses matrix exponential operations sequentially in every spatial direction. As a result, exponential matrices which are calculated and stored have small sizes, as those in the one-dimensional problem. This method decouples the exact evaluation of the linear part from the implicit treatment of the nonlinear reaction terms. We only solve a local nonlinear system at each spatial grid point. This method combines with the advantages of the spectral method and the compact implicit integration factor method, i.e., high precision, good stability, and small storage and soon. Numerical simulations show that it can have a great influence on the generation of patterns that the system control parameters take different values under otherwise identical conditions. The numerical results verify the theoretical results.【期刊名称】《物理学报》【年(卷),期】2018(067)005【总页数】10页(P50-59)【关键词】反应扩散方程;Gierer-Meinhardt模型;图灵斑图;Chebyshev谱方法【作者】张荣培;王震;王语;韩子健【作者单位】沈阳师范大学数学与系统科学学院,沈阳 110034;山东科技大学数学与系统科学学院,青岛 266590;沈阳师范大学数学与系统科学学院,沈阳 110034;沈阳师范大学数学与系统科学学院,沈阳 110034【正文语种】中文1 引言斑图是在空间或时间上具有某种规律性的非均匀宏观结构,普遍存在于自然界.1952年,著名的英国数学家图灵把他的目光转向生物学领域,用一个反应扩散系统成功地说明了某些生物体表面图纹产生的原理[1].图灵从数学角度表明,在反应扩散系统中,稳定状态会在某些条件下失稳,并自发产生空间定态图纹,此斑图通常称为图灵斑图. 经过多年的研究,各界学者利用反应扩散系统预测得到了更多的图灵斑图,在理论和实验方面取得了许多重要成果.他们证实了化学系统中图灵斑图的形成[2],讨论自催化反应中的动力学行为,探讨此类耦合反应扩散体系中影响图灵斑图的因素[3].给出Gray-Scott模型、Brusselator模型等系统扩散引起不稳定的数学机理[4],并描述了Gierer-Meinhardt,Lengyel-Epstein等模型的某些动力学行为(性质)[5,6].最近几年,图灵斑图在实验方面取得一系列最新的进展,Copie等[7]运用实验在一个双稳态被动非线性共振器中探讨了图灵调制和法拉第参数不稳定性的相互作用;Tompkins等[8]利用微流体化学室证实图灵理论体系,并观测到第七种时空模式;Lacitignola[9]研究了图灵不稳定现象的发生条件,论述了具体形态的电化学反应扩散模型在一个球面上的图案形成的特性;Gaskins等[10]在二氧化氯碘丙二酸反应实验中,通过添加卤化钠盐溶液得到新的图灵斑图.在这些系统中存在两种化学反应物质,它们不仅能相互作用,而且还能进行独自扩散.事实上,图灵斑图的产生对应的是一个非线性反应动力学过程与一种特殊扩散过程的耦合.这个特殊的扩散过程由于两种因子的扩散速度不同会发生失稳,这就是图灵斑图产生的机理.在数学上,图灵斑图可以用无量纲化的反应扩散方程组描述[11],即式中u和v是系统变量,分别代表参与化学反应的两种物质的浓度;c和d是扩散系数,t是时间变量,f(u,v)和g(u,v)表示反应项.设Ω为RN中带有光滑边界的有界区域,Ω=[0,a]×[0,b],边界为∂Ω,边界条件为齐次Neumann边界条件,即其中n表示边界上单位外法向.由于(1)式为耦合的非线性反应扩散方程,很难得到其精确解.近年来,许多学者用有限差分方法、有限元方法、谱方法等[12−14]多种数值方法求解(1)式,这些方法各有特点.相比于有限元方法和有限差分方法的低阶精度,谱方法[14]仅用少量的节点,采用Legendre,Chebyshev等适合的正交多项式离散即可达到指数阶收敛的谱精度.图灵斑图在空间上的结构具有一定的规律,且解比较光滑,因此采用谱方法离散是可行的.常用的谱配置方法主要有Fourier配置法[15],Chebyshev配置法[16],Hermite配置法等[17].由于本文考虑的(1)式边界条件为齐次Neumann边界条件,因此采用Chebyshev配置方法求解(1)式.对(1)式进行空间离散后,得到的是刚性的非线性常微分方程组(ODEs).显式时间离散方法虽可以用迭代的方法求解,但其对时间步长有严格的约束;隐式方法虽然可以允许大的时间步长,但是对于阶数非常大的非线性方程组的求解问题十分复杂,这对于全隐式方法来说是一个巨大的挑战.由于谱配置法所得到的谱微分矩阵是满的,显然利用追赶法等代数线性方程组的快速解法是不合适的,因此交替方向隐式方法在这里并不适用.本文采用紧致隐积分因子(compact implicit integration factor,cIIF)方法求解ODEs.2006年Nie等[18]以隐积分因子(IIF)方法为基础发展了cIIF方法.传统的隐积分因子方法在求解高维问题时,离散矩阵的指数运算的存储量和运算量非常大,导致运算速度缓慢.紧致隐积分因子方法[19]通过引入离散矩阵的紧致表达式并在各个方向进行矩阵的指数运算,使得中央处理器(CPU)的存储大大降低,计算速度也得到了显著提高.本文内容安排如下:第2节对反应扩散方程组进行线性分析,通过特征值解释图灵斑图的数学机理,然后以Gierer-Meinhardt模型为例分析系统处于稳定状态和不稳定状态时各参数需要满足的条件,进而探索斑图形成需要满足的条件;第3节研究数值方法,在空间离散条件下采用Chebyshev谱方法,时间离散条件下采用紧致隐积分因子方法,用MATLAB进行编程求解;第4节给出大量数值实验并对理论分析结果进行验证.2 图灵斑图的形成2.1 斑图形成的数学机理首先考虑(1)式没有扩散项,假设存在惟一的均匀定态解(u0,v0),即常数u0,v0满足令U=u−u0,V=v−v0,并在(u0,v0)处线性化后得到如下系统:式中c11=fu(u0,v0), c12=fv(u0,v0),c21=gu(u0,v0),c22=gv(u0,v0).均匀定态解(u0,v0)在没有扩散时是稳定的,这等价于相应的特征值问题的矩阵的特征值实部是负数.考虑加入扩散项后的反应扩散方程组((1)式).如果此时产生斑图,即(u0,v0)是不稳定的,要求特征值有正实部.所谓不稳定,体现为两种反应物的扩散速度不同,从而引起失稳.对(1)式作线性化处理,研究特征值正实部引起的线性不稳定性,进而推导出原方程的不稳定性.对均匀定态解(u0,v0)作一个微扰,可得线性微扰方程为求解如下方程可得相应的特征值:式中λ为特征值.只要(5)式中的特征值有正实部,则(u0,v0)对于(1)式是不稳定的.考虑到齐次Neumann边界条件,得到(5)式所对应的特征值为具体推导过程见附录A.2.2 Gierer-Meinhardt模型生物的发育过程是复杂的,其中重要的是形态形成阶段,与之对应的是生物体内器官的形成.由于该阶段的重要性,渐渐形成一个新的领域——形态学,主要研究导致细胞分化和定位因素的浓度对组织器官的影响.Gierer-Meinhardt模型是由Gierer和Meinhardt在研究激活物和抑制剂两种不同物质的产生和扩散时建立的[20],之后Gierer和Meinhardt利用数值方法导出一维和二维空间区域中上述系统产生多样斑图的条件.Gierer-Meinhardt模型被广泛应用于形态形成过程中一些基本现象的研究,最近的一些工作可以参见文献[21—23].以Gierer-Meinhardt模型为例,结合上述理论分析,计算产生斑图时需要满足的条件.取(1)式中其中系数κ,η,ε为系统的控制参数,固定η=0.1,ε=0.04.由此得到线性化系统(3)式中的系数为易得该系统的特征值为λ1= −1.2984,λ2=−7.7016,此时系统是稳定的.加入扩散项后,原方程组对应的特征问题为相应的特征方程为为使(8)式含有正实部的特征值,需要考虑两种情况.第一种情况是两个特征值异号,则应满足图1 特征值的实部Re(λ)随参数的变化(a)κ=0.0128;(b)κ=0.0152;(c)κ=0.008Fig.1.Real part Re(λ)of eigenvalues varying with parameters:(a)κ=0.0128;(b)κ=0.0152;(c)κ=0.008.经过化简可以得到此时应满足得0.0093248.第二种情况是两个特征值都是正的,应满足此时κ无解.由于反应扩散方程组联系于解析半群,所以线性化后的正实部特征值引起的不稳定性可以推导出原方程组的不稳定性.故当κ＞κ0=0.0093248时,系统处于不稳定状态,因而系统能够产生斑图.特征值的实部Re(λ)在参数κ取不同值时的变化如图1所示.从图1可以看出,当κ = 0.0128＞κ0和κ=0.0152＞κ0时,特征值的实部会出现正值,此时系统不稳定;当κ=0.008＜κ0时,特征值的实部始终为负,系统最后会达到稳定状态.第3节将用数值算例验证该结论.3 数值方法3.1 Chebyshev谱配置法将求解区域[−1,1]2离散为Gauss-Lobatto网格,即其中Nx和Ny是正整数.对于一般的求解区域Ω=[a,b]×[c,d],可以采用公式将区域转化为[−1,1]2.在网格Th中将u(x,y)数值解定义为矩阵形式,U∈R(Nx−1)×(Ny−1),式中ui,j表示u在网格点(xi,xj)的数值解.引入Chebyshev一阶微分矩阵和二阶微分矩阵(具体推导过程见附录B).则u(x,y)关于x的二阶偏导数在配置点的值,可以用矩阵乘积的形式近似为矩阵Ax是在Chebyshev二阶微分矩阵基础上考虑Neumann边界条件得到的,其中同样地,对于y的二阶偏导数,有UAy,其中矩阵Ay定义同Ax.借助谱微分矩阵,可将方程中的Laplace算子离散成矩阵乘积的形式,即将Chebyshev谱配置方法应用于反应扩散方程,得到其半离散形式为式中3.2 紧致隐积分因子法将对空间离散后得到的非线性常微分方程组((12)式)采用紧致隐积分因子方法进行时间离散.定义时间步长为τ=Δt,第n层时间步为tn=nτ,n=0,1,2,···. 在(12)式两端同时左乘指数矩阵e−Axt,右乘指数矩阵e−Ayt.为描述方便,取(10)式c=1,d=1,可将(12)式中第一个等式写为将时间离散为0=t0＜t1＜···,将(13)式在一个时间步长内关于时间积分,并用梯形公式近似可得二阶紧致隐积分因子格式为进一步化简得在非线性方程组(13)式中,右端第一项可以通过矩阵乘积得到,右端第二项采用Picard迭代方法求解:同理处理(12)式中第二个等式可得该方法中矩阵eAxΔt和eAyΔt的阶数分别为Nx×Nx和Ny×Ny.在空间网格剖分量很大时,该方法可以降低存储量和运算量,使计算速度更快.4 数值算例对于前述Gierer-Meinhardt模型,取Ω =(−1,1)× (−1,1),η=0.1,c=0.04,κ 是不固定的参数.设其中图2 取κ=0.0128时Gierer-Meinhardt模型形成的斑图(a)t=20;(b)t=80;(c)t=170;(d)t=270;(e)t=320;(f)t=340;(g)t=500;(h)t=600;(i)t= 900Fig.2.Turing patterns in Gierer-Meinhardt model whenκ=0.0128:(a)t=20;(b)t=80;(c)t=170;(d)t=270;(e)t=320;(f)t=340;(g)t=500;(h)t =600;(i)t=900.4.1 数值算例I取κ=0.0128,N=100,h=2/100=0.02,τ=0.1h,t取图2所示各值时,得到对应的图像.由图2可知,随着时间的推移,初始扰动不断增强扩大,最终形成清晰的斑图.4.2 数值算例II取κ=0.0152,t取图3所示各值,其他参数与算例I相同,可得到t取不同值时对应的图像.由图3可知,随着时间的推移,初始扰动不断增强扩大,最终形成清晰的斑图. 图3 取κ=0.0152时Gierer-Meinhardt模型形成的斑图(a)t=30;(b)t=80;(c)t=90;(d)t=140;(e)t=160;(f)t=290;(g)t=520;(h)t=620;(i)t=9 90Fig.3.Turing patterns in Gierer-Meinhardt model whenκ=0.0152:(a)t=30;(b)t=80;(c)t=90;(d)t=140;(e)t=160;(f)t=290;(g)t=520;(h)t= 620;(i)t=990.4.3 数值算例III取κ=0.008,其他取值与算例II相同,t取不同值时对应的图像如图4所示.由图4可知,随着时间的推移,系统达到稳定状态,反应扩散模型不能形成斑图.由数值模拟结果来看,其他条件一定的情况下,κ取不同值对于产生斑图有重要的影响.数值模拟结果与理论结果一致.此外,我们也对周期性边界条件的Gierer-Meinhardt模型采用Fourier谱方法进行数值求解,结果显示周期边界条件对斑图的形状几乎没有影响.5 结论介绍了图灵斑图形成的数学机理,并结合Gierer-Meinhardt模型,分析系统不稳定状态的各系数需要满足的条件,即产生斑图的条件.运用紧致隐积分因子方法大大减少了存储和CPU运算时间,该方法对于大时间数值模拟是一个高效、高精度的数值方法.数值算例模拟了斑图形成的过程,验证了理论分析结果.这些结论还可应用于求解带有分数阶的反应扩散方程组.图4 取κ=0.008时Gierer-Meinhardt模型形成的斑图(a)t=30;(b)t=80;(c)t=90;(d)t=140;(e)t=160;(f)t=220;(g)t=290;(h)t=270;(i)t=9 90Fig.4.Turing patterns in Gierer-Meinhardt model whenκ=0.008:(a)t=30;(b)t=80;(c)t=90;(d)t=140;(e)t=160;(f)t=220;(g)t=290;(h)t=2 70;(i)t=990.附录A 图灵斑图的形成机理首先在区域Ω⊂RN(N=1,2)内考虑带有齐次Neumann边界条件的Laplace算子的特征值问题.一维情况下,特征值问题为式中a∈R+.特征值问题可表示为µ2−λ=0,解得只有λ＜ 0时可解得特征值λk= −(kπ/a)2,且特征值所对应的特征函数为在二维情况下,特征值问题为式中a,b∈R+,应采用分离变量法求解特征值. 设u=X(x)Y(y),代入方程得设解得故特征值为λk,l=且特征值所对应的特征函数为考虑方程组的特征值问题,令代入原方程组可得设则方程组可化为当方程组(A3)有非零解,满足此时方程组所对应的特征值为附录B 谱微分矩阵定义在[−1,1]上的标准k阶Chebyshev多项式Tk(x)为 Tk(x)=cos(k arccosx),k=0,1,2,···. 令x=cosz,则有Tk=coskz,满足如下递推关系:Tk(x)在[−1,1]上的N+1个Gauss-Lobatto点值为零:设N阶多项式uN(x)∈PN在上述配置点xj满足uN(xj)=u(xj),则有式中hj(x)为N阶Lagrange基函数.用配置法求解未知量在网格点处的值,需要表示配置点处的导数值.对(B3)式求p阶导数,得式中系数从而可得一阶谱微分矩阵其中这里二阶谱微分矩阵可以由一阶谱微分矩阵平方得到,即参考文献[1]Turing A M 1952 Philos.Trans.R.Soc.Lond.B 2 37[2]Li X Z,Bai Z G,Li Y,Zhao K,He Y F 2013 Acta Phys.Sin.62 220503(inChinese)[李新政,白占国,李燕,赵昆,贺亚峰2013物理学报62 220503][3]Zhang L,Liu S Y 2007 Appl.Math.Mec.28 1102(in Chinese)[张丽,刘三阳2007应用数学和力学28 1102][4]Li B,Wang M X 2008 Appl.Math.Mec.29 749(in Chinese)[李波,王明新2008应用数学和力学29 749][5]Hu W Y,Shao Y Z 2014 Acta Phys.Sin.63 238202(in Chinese)[胡文勇,邵元智2014物理学报 63 238202][6]Peng R Wang M 2007 Sci.China A 50 377[7]Copie F,Conforti M,Kudlinski A,Mussot A,Trillo S 2016 Phys.Rev.Lett.116 143901[8]Tompkins N,Li N,Girabawe C,Heymann M,Ermentrout G B,Epstein IR,Fraden S 2014 A 111 4397[9]Lacitignola D,Bozzini B,Frittelli M,Sgura I 2017 Commun.Nonlinear Sci.Numer.Simul.48 484[10]Gaskins D K,Pruc E E,Epstein I R,Dolnik M 2016 Phys.Rev.Lett.117 056001[11]Zhang R P,Yu X J,Zhu J,Loula A 2014 Appl.Math.Model.38 1612[12]Zhang R P,Zhu J,Loula A,Yu X J 2016 put.Appl.Math.302 312[13]Bai Z G,Dong L F,Li Y H,Fan W L 2011 Acta Phys.Sin.60 118201(in Chinese)[白占国,董丽芳,李永辉,范伟丽2011物理学报60 118201][14]Zhang R,Zhu J,Yu X,Li M,Loula A F D 2017 put.310 194[15]Lv Z Q,Zhang L M,Wang Y S 2014 Chin.Phys.B 23 120203[16]Wang H 2010 mun.181 325[17]Hoz F D L,Vadillo F 2013 put.Phys.14 1001[18]Nie Q,Zhang Y T,Zhao R 2006 put.Phys.214 521[19]Nie Q,Wan F Y M,Zhang Y T,Liu X F 2008 put.Phys.227 5238[20]Gierer A,Meinhardt H 1972 Kybernetik 12 30[21]Ward M J,Wei J 2003 J.Nonlinear Sci.13 209[22]Wei J,Winter M 2004 J.Math.Pures Appl.83 433[23]Li H X 2015 J.Northeast Normal University 3 26(in Chinese)[李海侠2015东北师大学报3 26]。