[数据结构与算法] 图灵与图灵斑图

他是战争中的英雄。谈及图灵时不得不提他在二战时为盟军所作的杰出贡献。当时他在布莱切利公园（Bletchley Park）担任解码专家，于 1940 年创造出可以破译德军密报的机器 Bombe，为盟军的胜利立下了汗马功劳。二战时期，图灵带领团队通过破译德军密码（恩尼格玛密码机），让二战时间至少缩短了两年，拯救了国家，拯救了超过 1400 万人的生命…

计算机科学之父

1936 年，图灵在伦敦权威的数学杂志上发表了一篇划时代的重要论文《可计算数字及其在判断性问题中的应用》。文章里，图灵超出了一般数学家的思维范畴，完全抛开数学上定义新概念的传统方式，独辟蹊径，构造出一台完全属于想象中的 “计算机”，数学家们把它称为“图灵机”。

图灵提出用计算机去处理二进制的概念，并定义了最基本的计算机数学模型，冯诺依曼依据图灵的理论来设计出了物理模型，并实现之。现在的计算机，就是由当年的“图灵机”演变而来。1946年，图灵发表论文阐述存储程序计算机的设计。

人工智能之父

他在 19 世纪 50 年代发表论文《计算机器和智能》，提出用于判定机器是否具有智能的“图灵测试”，以及那句著名的提问： “机器会思考吗？ ”（Can Machines Think?）。他提出著名的“图灵测试”，指出如果第三者无法辨别人类与人工智能机器反应的差别，则可以论断该机器具备人工智能。

图灵斑图

他在这方面的贡献，鲜为人知。从 1952 年直到去世，图灵最后一两年一直在数理生物学方面做研究。他在 1952 年发表了一篇论文《形态发生的化学基础》(The Chemical Basis of Morphogenesis)，揭示了图灵斑图形成的机理。这也是今天我们要讲的重头戏。
在这里插入图片描述

同性恋、迫害与平反

图灵是个同性恋。早在在中学时期，他与另一个天资聪慧的男同学克里斯朵夫·莫尔卡姆（Christopher Morcom）成为至交，甚至有人认为，克里斯多夫是图灵的“初恋”。然而这位好朋友 1930 年不幸因肺结核去世。痛失挚友的图灵从此更加专注于科学研究，希望实现好朋友的遗志。

1951 年 12 月，接连发生的一系列事件对日后产生了很大的影响。图灵在曼彻斯特的街上遇见了一个年轻人阿诺德?穆雷。工人出身的穆雷正处于偷窃罪缓刑期，也没有工作。图灵和穆雷共进了午餐，一起回到了图灵的家里。在接下来的一个月，他们还相会了几次。

1952 年 1 月底，图灵发现住所遭窃。他报了警，警察检查了现场的指纹。图灵指控阿诺德?穆雷行窃，而穆雷声称自己是无辜的，并指认真正的罪犯是自己的一个旧相识哈里。警方也在图灵的住所找到了哈里的指纹。哈里彼时因为其他一些事情正在坐牢，在被问到图灵一案时，哈里向警方揭发了图灵和其朋友间一些很私密的情况。拔出萝卜带出泥，图灵把自己也带进去了。

1952 年 2 月 7 日，就在乔治六世驾崩，他的长女伊丽莎白继位的隔日，警方传讯了阿兰?图灵。在几轮审讯后，图灵承认了与穆雷之间的关系。这个供认让图灵遭受了牢狱之灾，因为根据 1885年的刑法修正案第11节：“任何男性，公开或私下，组织或参与组织，引诱或试图引诱其他男性进行严重猥亵的行为，都应该视为不法行为，并理应被法庭判处不超过2年、可带劳役或不带劳役的监禁。”

在公审中，法官给了图灵两个选择：要么坐牢，要么进行化学阉割来进行“治疗”——也就是雌激素注射。图灵选择了被阉割。“治疗” 的效果甚至影响到了图灵的行为和思考能力，他的手不断发抖，报纸上的字谜游戏，他再也没法解出…

可想而知，图灵生活在沉重的心理压力之下。我们不知道 1954 年 6 月 7 日的晚上发生了什么不一样的故事。我们也不知道是什么驱使图灵在睡觉前，将每晚都要吃的苹果蘸上了氰化物。第二天早晨，图灵被发现已经过世，年仅 41 岁。苹果创始人乔布斯终其一生对这位伟人极为崇拜，而至于苹果公司的 LOGO 是否是要致敬图灵，却众说纷纭——苹果公司从未公开承认和否认。

2009 年，英国计算机科学家康明（John Graham-Cumming）发起了为图灵平反的在线请愿，短短几日，请愿签名人数就超过了 3 万。在舆论的压力下，2009 年 9 月 10 日，时任英国首相戈登·布朗发表了给图灵的正式道歉声明。

2012 年 12 月，霍金、纳斯（Paul Nurse，诺贝尔医学奖得主）、里斯（Martin Rees，英国皇家学会会长）等11位重要人士向时任英国首相卡梅伦写去信件，要求正式为图灵平反。

2013 年 12 月 24 日，在英国司法大臣克里斯?格雷灵（Chris Grayling）的要求下，英国女王终于向图灵颁发了“皇家赦免”，取消了当年的指控，图灵也成为了二战以来第四名被皇家赦免的人！

2017年，艾伦·图灵法案生效，所有和图灵同时代因「同性恋」收到法律制裁的 4.9 万人被赦罪。

怪人和运动天赋

在布莱切利庄园这样由数学家和古典文学学者构成的不寻常人群中，图灵依然因性格怪异而博得了一定的名声。每年六月的第一周，图灵都会得一场严重的枯草热病，他会戴着军用防毒面具来遮蔽花粉，然后骑自行车去办公室。他的自行车有个毛病，车链每隔一定时间就会脱落。图灵不去修理它，而是数脚踏板转的圈数，然后赶在车链掉下之前下车，用手调整车链。后来他还在脚蹬子旁边安装了一个小巧的计数器，代替心中计数。

图灵的头衔有很多，传奇性的数学家、逻辑学家、现代计算机之父。但很少人知道他也是一位出色的跑者，甚至差点站上了奥运会的舞台。图灵自幼喜爱运动，战后更加热爱越野长跑。他经常参加业余高手们的越野长跑训练和竞赛。即使去参加学术会议，同事们都是搭乘公共交通工具，而图灵则舍代步而跑步，而且还是第一个到达会场。要不是因为受伤，他会代表英国去参加 1948 年的奥林匹克运动会的。图灵的最好成绩在当年（1947 年）足以排名世界前三。
在这里插入图片描述

图灵原本很有希望跑上奥运赛道，但因持续训练和略显僵硬的跑步姿势，他遭遇了腿部和髋部的伤病。1948 年 6 月 19 日，业余田协举办伦敦奥运马拉松选拔赛，图灵没有参加，只是在看台上观看了比赛。

英镑和电影

图灵登上了新版 50 英镑，在今年，也就是 2021 年他的 109 年诞辰，6 月 23 日开始发行。当初和他竞争的有，物理学家斯蒂芬·霍金、古生物学家玛丽·安宁、量子力学的奠基者之一保罗·狄拉克、物理化学家罗莎琳·富兰克林、拉马努金等等。

除了图灵的肖像之外，还印有图灵的标志性名言：This is only a foretaste of what is to come, and only the shadow of what is going to be。（这不过是将来之事的前奏，也是将来之事的影子。）这是 1949 年 6月11日图灵在接受《泰晤士报》采访时曾说过的话。而名言上面的签名是 1947 年，图灵在布莱切利公园的游客手册上留下的。

影片曾经想提到苹果的 Logo。多年来，人们一直猜测，Mac 上的苹果 Logo 与图灵相关。这是一个有趣的轶事，电影制作者也讨论过是否提这件事情。“我们讨论过，是否在影片最后显示一行相关的文字。一直探讨了好几个月，” Moore 说，“这是一个精彩的故事，乔布斯曾经公开否认过数次，但是其他人认为这是真事…… 我们希望这是真事，但是，我们没有足够的证据，可以安心地把它放在影片后面。从个人来说，我觉得这可能是真的，但是，我们无法证明它，所以需要慎重对待。”

图灵斑图（Turing Pattern）

自然和物理现象

**虎纹豹斑。**图灵看到了老虎和斑马的条纹、猎豹的斑点、鳄鱼牙齿排列的间距，具有非凡观察力的图灵就思考如何用数学解释这些现象。在 1952 年发表的最后一篇论文中，他创造性地使用反应 - 扩散系统的数学模型来描述自然界中的斑图，比如虎纹、豹斑。这类可以用反应 - 扩散方程描述的斑图，被后人称为图灵斑图（Turing Pattern）。
在这里插入图片描述

它的中文——图灵斑图，是 2019 年，经全国科学技术名词审定委员会审定发布的。

**沙漠波纹。**沙漠的波纹和图灵斑图非常类似。沙堆由风吹的沙粒沉积而成，随着一道沙堆变大，会从空气中吸收更多沙子，从而进一步促进自身生长。不过在长大过程中，沙堆从风中留下了沙子，也抑制了其他沙堆在旁边形成。沙漠中沙子涟漪纹路形成过程也类似于“激活剂一抑制剂”模型。在沙漠中，沙粒靠风移动，较多沙粒堆积的区域形成了脊。更高的脊容易从风中捕获沙粒，但同时也抑制了附近其他区域形成脊。这两个过程之间的平衡让沙子形成了均匀分布的涟漪图案。
在这里插入图片描述

**B-Z 振荡反应。**这样的现象在生物世界大都发展很慢，不易观察，但在化学上有一些直观的演示反应。比如 B-Z 振荡反应，常见以铈离子催化，丙二酸在稀硫酸水溶液中被溴酸盐氧化。在大容积的器皿内，我们会观察到溶液从黄色忽然变为透明，从透明又忽然变黄，如此周期振荡上千次才将底物消耗完，这展示了两组中间产物的对抗关系。

**胚胎。**激活剂能促进激活剂和抑制剂的合成过程，而抑制剂则抑制合成过程。图灵不清楚胚胎中的这两种物质到底是什么，可能是激素，也可能是基因，但他相信打破胚胎发育对称性，让胚胎自发形成特殊图案的就是这两种因子，并将其统称为“形态发生素”。

**铋晶体条纹。**在铋晶体生长的过程中，条纹的形成是由铋原子和下方的金属之间的作用力驱动的。铋原子想要填入金属分子晶格上的特定位置，但这些位置的间距对铋原子来说过于接近，这就像把一张照片塞进不够大的相框里一样，铋原子薄层会发生弯曲。弯曲带来的张力遵循一种波状起伏的模式，使一些原子凸起，形成条纹。在图灵方程中，铋原子在垂直方向上的移动（即偏离晶体平面）充当了激活剂，而铋原子在平面上的移动充当了抑制剂。这个过程中的形态发生素是原子的位移，而不是化学分子。这个最近的发现被发表在了自然物理杂志上面。
在这里插入图片描述

**种群。**人类的社交行为和群体的移动模式，也会无意之中展现出图灵斑图样。犯罪活动的分布就是一个典型例子。每座城市几乎都有一些地区的犯罪发生频率远超其他地区，这些区域被称为犯罪热点。美国加利福尼亚大学的数学家马丁·硕特在洛杉矶警局和长滩警局的犯罪活动地图数据的基础上，建立了一个犯罪扩散模型。这个模型中有两类犯罪热点：“超临界犯罪热点”和“亚临界犯罪热点”。硕特解释说，同一个地区的罪犯彼此之间会竞争，一种犯罪也会催生其他犯罪，这相当于“激活剂”。但长期看来，警察打击犯罪行为会抑制犯罪数量，所以警察活动就像“抑制剂”。在某个地区加强警力部署可以有效地消除亚临界犯罪热点，但只能将超临界犯罪热点转移到周围区域。因此，这个模型可以告诉警方，在某一区域严格巡查能不能起到打击罪犯的作用。模型还能在一定程度上预测犯罪热点将在何处形成，并给出较为合理的巡逻路线方案。

形成机理形而上的理解

底物消耗系统，其中一种化学物质（催化剂、活化剂、激活剂）消耗另一种抑制性化学物质（抑制剂、底物），这就像一片森林由于干旱有许多起火点，如果迅速扑灭火苗，就只会让一些区域沦为焦土，火势不会蔓延到整片森林。这种催化剂，能自动激活、并随后产生抑制剂，使得它比催化剂的传播速度更快，导致催化剂中途停止，最终生成相应的图案。

再举一个形象的例子，如果把激活剂比作兔子，那么抑制剂就好像狐狸。兔子多了，能生更多的兔子，狐狸的食物也多，狐狸的数量会增加。但如果狐狸太多，把兔子吃光了，狐狸的数量最终依然会下降。而且狐狸往往比兔子跑得更快。这两者的数量之比是动态变化的。兔子和狐狸数量增加，它们所占的栖息地面积也会增加。为了获得足够的食物，狐狸会彼此分散开来，避免与同类竞争。狐狸向外迁徙的速度要大于兔子迁徙的速度。用不了多久，狐狸会优先占据外围区域，共同將兔子包围起来。最终，就形成了一个个圆形的激活剂（兔子）区域，和围绕这些圆形图案的抑制剂（狐狸）区域。如果“兔子”的数量超过了周围的“狐狸”，在动物皮肤上就会出现一些诸如色素沉积的变化。图灵认为，通过不断重复这个过程，就能形成重复的同样斑纹。

案例来说，如果系统是稳定的，形成的图像一定具有某种对称性。然而，可以观察到，虎斑豹纹总是有参差不齐的美，其实在这个过程中一点不平衡就能产生强烈的不平衡，所以才形成了这种参差的美。

反应扩散系统

形成条件

只要改变模型中的一个变量，比如激活剂和抑制剂的扩散速度，或者反应区域的总面积，就能得到包括圆点、条纹和六边形在内的几乎所有自然界生物所拥有的图案。这正是图灵模型的奇妙之处。要形成图案必须具备一个条件：抑制剂在溶液中的扩散速度必须要比激活剂更快。图灵结构产生的必要条件，就是两个反应物的扩散系数之差要达到一个数量级以上。

想要实现斑图就需要考虑其时空性质，即时间和空间的转变导致斑图的形成。时空离散的捕食者系统斑图形成的条件具体来说为以下三点：

图灵斑图的数值方法之参数化有限元方法

简介

基本想法是利用表面有限元数值方法求解演化曲面上的反应扩散方程。不断增长的生物表面上的模式形成，是通过求解曲面上的反应扩散方程得到的，可以参考图灵的经典文章 ¹。这种方程表现出空间均匀结构的扩散驱动不稳定性，导致空间不均匀的图案。这项工作的进一步的发展，有了更多的数学模型，其中的一些应用包括九头蛇的图案形成、脊椎动物的连续图案形成、动物的皮肤表面的颜色标记、蝴蝶翅膀上的色素沉着图案等等。可以参考综述文章 ²。

图灵的原始工作没有考虑几何的增长和变化。然而，这些因素在自然界中所观察到的模式的发展中很重要。增长的区域上的模式形成已被广泛研究。例如 Crampin 等人 ³ 考虑一维区域的区域增长，并表明它可能是一种机制，用于增加模式形成时，相对于初始条件的稳健性。在了解在一种名叫 Pomacanthus 的鱼中观察到的条纹进化的背景下（[^近藤和浅井 1995])，二维增长平面区域被考虑⁴。也可以看看⁵，其中图灵扩散驱动的不稳定性条件被推广到具有缓慢、各向同性增长的反应扩散系统。本文中提出的数值方法是一种在演化表面上探索此类反应动力学的工具。

许多有趣的工作都是关于增长平面区域的。然而，在生物的应用中，很自然地考虑在不断发展的空间三维区域的表面上形成图案，例如我们参考⁶，和⁷等，其中研究了曲率
对固定球体和生长锥体的图案形成的影响。另一项初步研究是在这篇文章 ⁸ 中，它提出了一种数值方法，用于解决球形表面图案形成问题，并应用于实体肿瘤生长。他们的数值方法基于参数化球体上的表面，然后使用线方法。

下面提到的方法基于表面有限元方法，它是有限元方法⁵的自然延伸，并且能够处理复杂的几何形状和形状⁹¹⁰[^艾克斯和艾略特 2008]¹¹¹²¹³。这个想法是对表面进行三角剖分，并使用基于三角剖分的分段线性表面有限元空间来逼近偏微分方程组。在这种情况下，三角剖分的顶点以规定的或受某些进化规律支配的速度移动。为了做到这一点，我们需要在表面给出一个适当的守恒定律。当以适当的变分形式给出时，演化表面有限元方法利用了该守恒定律的特殊特征。

首先，我们给出演化表面上的反应扩散方程。在这个框架内，我们给出了表面梯度的符号表达，这是演化表面上的第一原理反应扩散系统推导出来的关键。我们考虑经过充分研究的活化剂耗尽衬底模型[^吉尔和迈因哈特 1972]¹⁴¹⁵，它也被称为 Brusselator 模型。同样，在此框架中可以轻松处理任何其他合理的反应动力学。我们还讨论了表面生长模型的建模。

其次，我们看看一个演化表面有限元方法，该方法应用于改变形状的演化表面上的反应扩散系统。变分公式的一个特点是不会出现某些几何量，例如曲面的平均曲率和法线。然后通过有限元方法以自然的方式利用这一点。在这种方法中，在导出反应扩散系统或利用表面有限元方法时不需要转换或参数化。在空间中离散会产生一个非自治常微分方程组，该方程组在时间上是离散的。特别是因为这个应用涉及半线性抛物方程，我们线性化非线性动力系统[^马兹瓦缪斯]，以获得线性代数系统。然后使用 GMRes 算法¹⁶求解。

最后，在固定和进化的表面上呈现出图灵模式图案。我们展示了表面有限元方法的普遍适用性，该方法具有处理生物模型可能导致的各向异性生长的能力。结果证实，结合区域增长增强了模式选择过程。该方法也适用于将表面演化与表面上的反应扩散系统耦合的模型。这在再下一个章节中有说明，它应用于模拟实体瘤的生长。

部分结论性意见在最后一个章节。其他也考虑表面偏微分方程数值解的文章的例子包括¹⁷¹⁸¹⁹²⁰²¹。这篇文章²² 对计算演化表面的不同方法的回顾。

演化表面上反应扩散方程的推导

曲面梯度

我们假设

\Gamma

可以全局表示为函数

d

的零水平集，即存在开集

\mathcal{U} \subset \mathbb{R}^{3}

和函数

\in C^{2}(\mathcal{U})

，使得

\mathcal{U} \cap \Gamma=\{x \in \mathcal{U}: d(x)=0\}, \quad \text { and } \quad \nabla d(x) \neq 0, \quad \forall x \in \mathcal{U} \cap \Gamma

我们定义单位法向量场为，

\boldsymbol{v}(x)=\frac{\nabla d(x)}{|\nabla d(x)|}

我们定义切向梯度为，

\nabla_{\Gamma} \eta(x)=\nabla \eta(x)-\nabla \eta(x) \cdot \boldsymbol{v}(x) \boldsymbol{v}(x), \quad x \in \Gamma

它只依赖于

\eta

在曲面上的值。它具有三个分量，表示为：

\nabla_{\Gamma} \eta(x)=\left(\underline{D}_{1} \eta(x), \underline{D}_{2} \eta(x), \underline{D}_{3} \eta(x)\right)

Laplace-Beltrami 算子就定义为切向梯度的切向散度：

\Delta_{\Gamma} \eta(x)=\nabla_{\Gamma} \cdot \nabla_{\Gamma} \eta(x)=\sum_{i=1}^{3} \underline{D}_{i} \underline{D}_{i} \eta(x)

演化表面上的反应扩散系统

让

\Gamma(t)

是嵌入到三维空间

\Omega(t)

的二维曲面。曲面上每一点的速度表示为，

\boldsymbol{v}:=V \boldsymbol{v}+\boldsymbol{v}_{T}

\boldsymbol{u}=\left\{u_{i}\right\}_{i=1}^{m}

是个矢量场，

\mathcal{R}(t)

是曲面的某一部分，由物质守恒定律，

\frac{d}{d t} \int_{\mathcal{R}(t)} u_{i}=-\int_{\partial \mathcal{R}(t)} \boldsymbol{q}_{i} \cdot \boldsymbol{\mu}+\int_{\mathcal{R}(t)} f_{i}(\boldsymbol{u})

这里边的

\boldsymbol{q}_{i}

表示曲面通量，

f_{i}(\boldsymbol{u})

表示曲面内的净生产率。

\frac{d}{d t} \int_{\mathcal{R}(t)} u_{i}=-\int_{\partial \mathcal{R}(t)} \boldsymbol{q}_{i} \cdot \boldsymbol{\mu}+\int_{\mathcal{R}(t)} f_{i}(\boldsymbol{u})

通过分部积分公式，

\left(\nabla_{\Gamma} \cdot \boldsymbol{w}, v\right)_{\Gamma}=\left(\mathbf{n}_{\partial \Gamma} \cdot \boldsymbol{w}, v\right)_{\partial \Gamma}-\left(\boldsymbol{w}, \nabla_{\Gamma} v\right)_{\Gamma}+\left(\boldsymbol{w}, H \boldsymbol{n}_{\Gamma} v\right)_{\Gamma}

我们得到，

\int_{\partial \mathcal{R}(t)} \boldsymbol{q} \cdot \boldsymbol{v}=\int_{\mathcal{R}(t)} \nabla_{\Gamma} \cdot \boldsymbol{q}+\int_{\mathcal{R}(t)} \boldsymbol{q} \cdot \boldsymbol{v} H=\int_{\mathcal{R}(t)} \nabla_{\Gamma} \cdot \boldsymbol{q}

同样，由经典的雷诺输运公式，

\frac{\mathrm{d}}{\mathrm{d} t} \int_{\Omega(t)} f \mathrm{dx}=\int_{\Omega(t)} \dot{f}+f \operatorname{div} \mathbf{v} \mathrm{dx}

我们有，

\frac{d}{d t} \int_{\mathcal{R}(t)} \eta=\int_{\mathcal{R}(t)} \partial^{\bullet} \eta+\eta \nabla_{\Gamma} \cdot \boldsymbol{v}

联立上述结果，我们有，

\int_{\mathcal{R}(t)} \partial^{\bullet} u_{i}+u_{i} \nabla_{\Gamma} \cdot \boldsymbol{v}+\nabla_{\Gamma} \cdot \boldsymbol{q}_{i}=\int_{\mathcal{R}(t)} f_{i}(\boldsymbol{u})

由

\mathcal{R}(t)

的任意性，我们得到，

\partial^{\bullet} u_{i}+u_{i} \nabla_{\Gamma} \cdot \boldsymbol{v}+\nabla_{\Gamma} \cdot \boldsymbol{q}_{i}=f_{i}(\boldsymbol{u}), \quad \text { on } \Gamma(t)

我们假设

\mathcal{D}

是一个对角的扩散张量矩阵，即

\mathcal{D}_{i j}=d_{i} \delta_{i j}

，并且假设化学物质根据菲克定律扩散，即

\boldsymbol{q}_{i}=-\mathcal{D}_{i j} \nabla_{\Gamma} u_{j}

那么，方程就变成了，

\partial^{\bullet} u_{i}+u_{i} \nabla_{\Gamma} \cdot \boldsymbol{v}=\nabla_{\Gamma}\left(\mathcal{D}_{i j} \nabla_{\Gamma} u_{j}\right)+f_{i}(\boldsymbol{u}), \quad \text { on } \Gamma(t)

写成向量的形式，演化曲面上的反应扩散方程系统可以写为，

\partial^{\bullet} \boldsymbol{u}+\boldsymbol{u} \nabla_{\Gamma} \cdot \boldsymbol{v}=D \Delta_{\Gamma} \boldsymbol{u}+\boldsymbol{f}(\boldsymbol{u})

采用的初值条件是，

\boldsymbol{u}(\cdot, 0)=\boldsymbol{u}_{0}(\cdot) \quad \text { on } \Gamma(0)

注意到，这里的速度场可能是依赖于曲面本身，也可能是依赖于反应扩散方程组的解。

模型说明

为了简单起见，我们考虑一个非常简单的模型，叫做 Brusselator 模型（活化剂耗尽的衬底模型）。假设

\boldsymbol{u}=(u, w)

是二维的，表示两种化学物质的浓度。

\boldsymbol{f}(\boldsymbol{u})=\left(f_{1}, f_{2}\right)^{T}=\left(\gamma\left(a-u+u^{2} w\right), \gamma\left(b-u^{2} w\right)\right)^{T}

对应于

\boldsymbol{f}(\boldsymbol{u})=\mathbf{0}

，有唯一的解，

\quad w=b /(a+b)^{2}

这刚好对应了不动曲面的唯一稳态。对于给定的模型，可能存在一组参数，其值属于图灵空间，对于这些参数，固定取域上的模型由于空间均匀结构的扩散驱动不稳定性而展现出空间模式的形成。区域增长增加了生物形态发生对的范围，与固定区域相比，它们更有可能诱发图灵模式。

曲面有限元方法

下面考虑在演化曲面上，使用曲面有限元方法求解反应扩散系统。

变分形式

对于，

\partial^{\bullet} u_{i}+u_{i} \nabla_{\Gamma} \cdot \boldsymbol{v}=d_{i} \Delta_{\Gamma} u_{i}+f_{i}(\boldsymbol{u})

使用标准的有限元变分，可以得到，

\begin{aligned} \int_{\Gamma(t)} f_{i}(\boldsymbol{u}) \varphi &=\int_{\Gamma(t)} \partial^{\boldsymbol{}} u_{i} \varphi+u_{i} \varphi \nabla_{\Gamma} \cdot v-d_{i} \int_{\Gamma(t)} \varphi \Delta_{\Gamma} u_{i} \\ &=\int_{\Gamma(t)} \partial^{\bullet} u_{i} \varphi+u_{i} \varphi \nabla_{\Gamma} \cdot v+d_{i} \int_{\Gamma(t)} \nabla_{\Gamma} u_{i} \cdot \nabla_{\Gamma} \varphi-\int_{\partial \Gamma(t)} \varphi \nabla_{\Gamma} u_{i} \cdot \mu \end{aligned}

不妨假设曲面是没有边界的，那么上面的最后一项就可以丢掉了，进一步利用输运定理，

\begin{aligned} \int_{\Gamma(t)} f_{i}(\boldsymbol{u}) \varphi &=\int_{\Gamma(t)} \partial^{\bullet} u_{i} \varphi+u_{i} \varphi \nabla_{\Gamma} \cdot \boldsymbol{v}+d_{i} \int_{\Gamma(t)} \nabla_{\Gamma} u_{i} \cdot \nabla_{\Gamma} \varphi \\ &=\int_{\Gamma(t)} \partial^{\bullet}\left(u_{i} \varphi\right)-u_{i} \partial^{\bullet} \varphi+u_{i} \varphi \nabla_{\Gamma} \cdot \boldsymbol{v}+d_{i} \int_{\Gamma(t)} \nabla_{\Gamma} u_{i} \cdot \nabla_{\Gamma} \varphi \\ &=\frac{d}{d t} \int_{\Gamma(t)} u_{i} \varphi-\int_{\Gamma(t)} u_{i} \partial^{\bullet} \varphi+d_{i} \int_{\Gamma(t)} \nabla_{\Gamma} \varphi \cdot \nabla_{\Gamma} u_{i} \end{aligned}

变分形式就变成了，寻找

u_{i} \in H^{1}(\Gamma(t))

，使得

\begin{aligned} &\frac{d}{d t} \int_{\Gamma(t)} u_{i} \varphi-\int_{\Gamma(t)} u_{i} \partial^{\bullet} \varphi+d_{i} \int_{\Gamma(t)} \nabla_{\Gamma} \varphi \cdot \nabla_{\Gamma} u_{i} \\ &=\int_{\Gamma(t)} f_{i}(\boldsymbol{u}) \varphi, \quad \forall \varphi \in H^{1}(\Gamma(t)) \end{aligned}

演化曲面有限元方法

对于参数化的有限元方法，我们选择三角逼近面片的顶点的运动来刻画曲面的运动，

\dot{X}_{j}(t)=v\left(X_{j}(t), t\right)

有限元空间表示如下，

S_{h}(t)=\left\{\phi \in C^{0}\left(\Gamma_{h}(t)\right):\left.\phi\right|_{T_{k}} \text { is linear affine for each } T_{k} \in \mathcal{T}_{h}(t)\right\}

我们用

\left\{\chi_{j}(\cdot, t)\right\}_{j=1}^{N}

来表示移动的节点基函数。

我们定义离散的物质速度，将其表示为节点速度的线性插值，即

\boldsymbol{v}_{h}=\sum_{j=1}^{N} \dot{X}_{j}(t) \chi_{j}

并且定义离散的物质导数，

\partial_{h}^{\bullet} \phi=\phi_{t}+v_{h} \cdot \nabla \phi

使用标准代数离散，我们可以得到半离散的代数系统，

\frac{d}{d t}\left(\mathcal{M}(t) \boldsymbol{\alpha}_{i}\right)+d_{i} \mathcal{S}(t) \boldsymbol{\alpha}_{i}=\boldsymbol{F}_{i}(t)

其中，

\mathcal{M}(t)_{j k}=\int_{\Gamma_{h}(t)} \chi_{j} \chi_{k}

\mathcal{S}(t)_{j k}=\int_{\Gamma_{h}(t)} \nabla_{\Gamma_{h}} \chi_{j} \cdot \nabla_{\Gamma_{h}} \chi_{k}

\boldsymbol{F}_{i j}=\int_{\Gamma_{h}(t)} f_{i}(\boldsymbol{U}) \chi_{j}

求解这样一个代数系统，本质上就得到了原问题的离散逼近。

时间离散

考虑最简单的两种物质浓度，

\boldsymbol{u}=(u, w)

，它的时间离散表示为

\left(U^{n}, W^{n}\right)

，令

d_{1}=1, d_{2}=d

，则时间离散可以写为，

\left\{\begin{array}{l} \frac{1}{\tau} \int_{\Gamma_{h}^{n+1}} U^{n+1} \chi_{j}^{n+1}+\int_{\Gamma_{h}^{n+1}} \nabla_{\Gamma_{h}^{n}} U^{n+1} \cdot \nabla_{\Gamma_{h}^{n}} \chi_{j}^{n+1} \\ =\frac{1}{\tau} \int_{\Gamma^{n}} U^{n} \chi_{j}^{n}+\int_{\Gamma_{h}^{n+1}}^{n} f_{1}\left(U^{n+1}, W^{n+1}\right) \chi_{j}^{n+1} \\ \frac{1}{\tau} \int_{\Gamma_{h}^{n+1}} W^{n+1} \chi_{j}^{n+1}+d \int_{\Gamma_{h}^{n+1}} \nabla_{\Gamma_{h}^{n}} W^{n+1} \cdot \nabla_{\Gamma_{h}^{n}} \chi_{j}^{n+1} \\ =\frac{1}{\tau} \int_{\Gamma^{n}} W^{n} \chi_{j}^{n}+\int_{\Gamma_{h}^{n+1}} f_{2}\left(U^{n+1}, W^{n+1}\right) \chi_{j}^{n+1} \end{array}\right.

将

f

的表达式代入，再利用简单的线性化技巧，

\left(U^{n+1}\right)^{2} \approx U^{n} U^{n+1}

最后得到一个线性化系统，

\left\{\begin{array}{c} \left(\frac{1}{\tau}+\gamma\right) \int_{\Gamma_{h}^{n+1}} U^{n+1} \chi_{j}^{n+1}+\int_{\Gamma_{h}^{n+1}} \nabla_{\Gamma_{h}^{n}} U^{n+1} \cdot \nabla_{\Gamma_{h}^{n}} \chi_{j}^{n+1} \\ -\gamma \int_{\Gamma_{h}^{n+1}} U^{n} W^{n} U^{n+1} \chi_{j}^{n+1}=\frac{1}{\tau} \int_{\Gamma^{n}} U^{n} \chi_{j}^{n}+\gamma a \int_{\Gamma_{h}^{n}} \chi_{j}^{n}, \\ \begin{array}{c} \frac{1}{\tau} \int_{\Gamma_{h}^{n+1}} W^{n+1} \chi_{j}^{n+1}+d \int_{\Gamma}^{n+1} \nabla_{\Gamma_{h}^{n}} W^{n+1} \cdot \nabla_{\Gamma_{h}^{n}} \chi_{j}^{n+1} \\ +\gamma_{\Gamma}^{n+1}\left(U^{n+1}\right)^{2} W^{n+1} \chi_{j}^{n+1}=\frac{1}{\tau} \int_{\Gamma^{n}} W^{n} \chi_{j}^{n}+\gamma b \int_{\Gamma_{h}^{n}} \chi_{j}^{n} \end{array} \end{array}\right.

代数系统的形式为，

\left\{\begin{array}{l} \left(\left(\frac{1}{\tau}+\gamma\right) \mathcal{M}^{n+1}+\mathcal{S}^{n+1}-\gamma \mathcal{M}_{1}^{n+1}\right) \boldsymbol{U}^{n+1}=\frac{1}{\tau} \mathcal{M}^{n} \boldsymbol{U}^{n}+\boldsymbol{F}_{1}^{n} \\ \left(\frac{1}{\tau} \mathcal{M}^{n+1}+d \mathcal{S}^{n+1}+\gamma \mathcal{M}_{2}^{n+1}\right) \boldsymbol{W}^{n+1}=\frac{1}{\tau} \mathcal{M}^{n} \boldsymbol{W}^{n}+\boldsymbol{F}_{2}^{n} \end{array}\right.

其中，

\begin{aligned} \mathcal{M}_{i j}^{n} &=\int_{\Gamma_{h}^{n}} \chi_{i}^{n} \chi_{j}^{n}, \quad \mathcal{M}_{1 i j}^{n}=\int_{\Gamma_{h}^{n}} U^{n} W^{n} \chi_{i}^{n} \chi_{j}^{n}, \quad \mathcal{M}_{2 i j}^{n}=\int_{\Gamma_{h}^{n}}\left(U^{n+1}\right)^{2} \chi_{i}^{n} \chi_{j}^{n} \\ \mathcal{S}_{i j}^{n} &=\int_{\Gamma_{h}^{n}} \nabla_{\Gamma_{h}^{n}} \chi_{i}^{n} \cdot \nabla_{\Gamma_{h}^{n}} \chi_{j}^{n}, \quad \boldsymbol{F}_{1 i}^{n}=\gamma a \int_{\Gamma_{h}^{n}} \chi_{i}^{n}, \quad \boldsymbol{F}_{2 i}^{n}=\gamma b \int_{\Gamma_{h}^{n}} \chi_{i}^{n} \end{aligned}

该方法需要对初始曲面进行三角剖分。例如，球体按以下方式进行三角剖分。使用均匀放置在球体上的 6 个点手动获得非常粗略的三角剖分。然后通过将每个三角形细分。然后将新顶点投影到球体表面。重复这种细化，直到获得可接受的最大网格尺寸。

肿瘤增长

考虑肿瘤生长模型来说明 ESFEM 在耦合曲面演化与曲面反应扩散的系统中的应用。

模型建立

该模型的曲面演化是化学驱动型的。这里

u

表示生长促进因子，

w

表示生长抑制因子。此外，我们认为表面演化受的法向速度控制，表达如下，

V=\delta u-\epsilon H

其中

\delta

是生长速率，

\epsilon

可以被认为是模拟肿瘤表面张力的参数。我们也可以将

\epsilon

视为一个小的数学上的参数，它规范了进化规律。这很有用，因为很容易看出，即使在法线方向上表面速度恒定的简单情况下，光滑表面也可能在短时间内形成拐角。ESFEM 的一个优点是我们只需要三角形顶点的位置来进行计算，从而导致该方法具有吸引力的灵活性。

数值方法

数值方法分成两步。
第一步，把曲面的线性逼近写成基函数的线性组合，组合系数为节点的位置，然后，使用标准的有限元方法，根据当前逼近曲面顶点的位置，得到下一步逼近曲面的顶点的位置。

\left(\frac{1}{\tau} \mathcal{M}^{n}+\epsilon \mathcal{S}^{n}\right) \boldsymbol{X}^{n+1}=\frac{1}{\tau} \mathcal{M} \boldsymbol{X}^{n}+\boldsymbol{G}^{n}

\boldsymbol{G}_{i}^{n}=\int_{\Gamma^{n}} \delta\left(\boldsymbol{U}_{h}^{n}\right)_{1} \boldsymbol{v}^{n} \chi_{i}

其中，

\boldsymbol{X}^{n}

表示节点坐标。
第二步，根据得到的下一步的分片曲面的顶点位置，我们可以得到新的基函数，利用基函数计算刚度矩阵和质量矩阵，进而执行之前提到的离散格式：

数值结果

网格自适应和重新划分网格

原则上一段时间后，由于三角形的顶点根据速度规律演化，演化表面的三角剖分可能不适用于有限元模拟。例如，它们可以出现非常小的角度从而出现的某种意义上退化。另一方面，三角剖分可能无法解析浓度快速变化的子区域中的解的特征。由于夹断和自相交，表面的拓扑结构也可能发生变化。

上面提出的模拟中，我们要保证初始三角剖分足够细足够均匀足以进行计算。可以通过以自适应方式细化和粗化三角剖分来执行更有效的计算，以反映曲面偏微分方程的解的行为并避免三角剖分的退化。

[^近藤和浅井 1995]:Kondo S, Asai R (1995) A reaction-diffusion wave on the skin of the marine anglefish, Pomacanthus. Nature 376:765–768

[^艾克斯和艾略特 2008]:Eilks C, Elliott CM (2008) Numerical simulation of dealloying by surface dissolution via the evolving surface finite element method. J Comp Phys 227:9727–9741

[^吉尔和迈因哈特 1972]:Gierer A, Meinhardt H (1972) A theory of biological pattern formation. Kybernetik 12:30–39

图灵斑图的 Trace 方法

\mathbf{w}=(\delta u-\epsilon H) \mathbf{n}

\frac{\partial u_{i}}{\partial t}+\mathbf{w} \cdot \nabla u_{i}+u_{i} \nabla_{\Gamma} \cdot \mathbf{w}=d_{i} \Delta_{\Gamma} u_{i}+f_{i}(\mathbf{u}), \quad i=1,2

\mathbf{f}(\mathbf{u})=\left(f_{1}, f_{2}\right)^{T}=\left(\gamma\left(a-u+u^{2} w\right), \gamma\left(b-u^{2} w\right)\right)^{T}

\phi_{t}+\left(\delta u-\epsilon \operatorname{div} \frac{\nabla \phi}{|\nabla \phi|}\right)|\nabla \phi|=0

d(x)=\left\{\begin{array}{rr} -\operatorname{dist}(x, \Gamma), & \text { x belong to interior of } \Gamma \\ \operatorname{dist}(x, \Gamma), & \text { otherwise } \end{array}\right.

\phi_{t}-\epsilon \Delta \phi=-\delta u

\int_{\Omega} \phi_{t} \psi+\epsilon \int_{\Omega} \nabla \phi \cdot \nabla \psi=-\delta \int_{\Omega} u \psi+\epsilon \int_{\partial \Omega} \nabla \phi \cdot \overline{\mathbf{n}} \psi

H^{1, e}=\left\{\left.u\right|_{\Gamma} \in H^{1}(\Gamma) \mid \nabla u \cdot \nabla d=0\right\}

\begin{gathered} \int_{\Gamma}\left(\frac{\partial u_{i}}{\partial t} v+d_{i} \nabla_{\Gamma} u_{i} \cdot \nabla_{\Gamma} v+(\mathbf{w} \cdot \mathbf{n}) u_{i}\right) \mathrm{d} \mathbf{s}=\int_{\Gamma} f_{i}(\mathbf{u}) v \mathrm{~d} \mathbf{s}, \quad \forall v \in H^{1, e} \\ \int_{\Gamma} u_{t} v+d_{1} \int_{\Gamma} \nabla_{\Gamma} u \cdot \nabla_{\Gamma} v+\int_{\Gamma}(\delta u-\epsilon H) u v=\gamma \int_{\Gamma}\left(a-u+u^{2} w\right) v, \quad \forall v \in H^{1, e} \\ \int_{\Gamma} w_{t} v+d_{2} \int_{\Gamma} \nabla_{\Gamma} w \cdot \nabla_{\Gamma} v+\int_{\Gamma}(\delta u-\epsilon H) w v=\gamma \int_{\Gamma}\left(b-u^{2} w\right) v, \quad \forall v \in H^{1, e} \end{gathered}

\begin{array}{r} \int_{\Gamma^{n+1}}\left(1+\frac{1}{\Delta t}-\varepsilon H\right) u^{n+1} v+d_{1} \int_{\Gamma^{n+1}} \nabla_{\Gamma} u^{n+1} \cdot \nabla_{\Gamma} v \\ =-\delta \int_{\Gamma^{n+1}}\left(u^{n}\right)^{2} v+a \gamma \int_{\Gamma^{n+1}} v+\gamma \int_{\Gamma^{n+1}}\left(u^{n}\right)^{2} w^{n} v+\int_{\Gamma^{n+1}} \frac{u^{n}}{\Delta t} v \\ \int_{\Gamma^{n+1}}\left[\left(u^{n+1}\right)^{2}+\frac{1}{\Delta t}+\delta u^{n+1}-\varepsilon H\right] w^{n+1} v+d_{2} \int_{\Gamma^{n+1}} \nabla_{\Gamma} w^{n+1} \cdot \nabla_{\Gamma} v \\ =b \int_{\Gamma^{n+1}} v+\int_{\Gamma^{n+1}} \frac{w^{n}}{\Delta t} v \end{array}

\begin{gathered} {\left[\left(1+\frac{1}{\Delta t}\right) \mathcal{M}^{n+1}-\varepsilon \mathcal{M}_{1}^{n+1}+d_{1} \mathcal{S}^{n+1}+\rho_{n+1} \mathcal{G}^{n+1}\right] \vec{u}_{h}^{n+1}=\mathcal{F}_{1}^{n+1}} \\ \mathcal{M}_{i j}^{n}=\int_{\Gamma_{h}^{n}} \chi_{i} \chi_{j} \\ \mathcal{S}_{i j}^{n}=\int_{\Gamma_{h}^{n}} \nabla_{\Gamma} \chi_{i} \nabla_{\Gamma} \chi_{j} \\ \mathcal{M}_{1 i j}^{n}=\int_{\Gamma_{h}^{n}} H \chi_{i} \chi_{j} \\ \mathcal{F}_{2 i}^{n+1}=\int_{\Gamma_{h}^{n+1}}\left(b+\frac{w_{h}^{n}}{\Delta t}\right) \chi_{i} \\ \left.\left.\mathcal{F}_{1 i}^{n+1}=\int_{\Gamma_{h}^{n+1}}^{n+1}=\int_{\mathcal{O}}\left(u_{h}^{n}\right)^{2}+a \gamma+\gamma\left(u_{h}^{n}\right)^{2} w_{h}^{n}+\frac{u_{h}^{n}}{\Delta t}\right] \chi_{i}^{n+1} \nabla_{\Gamma} \chi_{i}\right)\left(\mathbf{n}_{h}^{n+1} \nabla_{\Gamma} \chi_{j}\right) \\ \left.\frac{1}{\Delta t} \mathcal{M}^{n+1}-\varepsilon \mathcal{M}_{1}^{n+1}+\mathcal{M}_{2}^{n+1}+d_{2} \mathcal{S}^{n+1}+\rho_{n+1} \mathcal{G}^{n+1}\right] \vec{w}_{h}^{n+1}=\mathcal{F}_{2}^{n+1} \end{gathered}

Turing A (1952) The chemical basis of morphogenesis. Phil Trans R Soc Lond B 237:37–72 ??
Murray JD (2002) Mathematical biology I and II, 3rd edn. Springer, Berlin ??
Crampin EJ, Gaffney EA, Maini PK (1999) Reaction and diffusion on growing domains: scenarios for robust pattern formation. Bull Math Biol 61:1093–1120 ??
Maini PK, Painter KJ, Chau HNP (1997) Spatial pattern formation in chemical and biological systems. Faraday Trans 93:3601–3610 ??
Madzvamuse A, Maini PK, Wathen AJ (2003) A moving grid finite element method applied to a model biological pattern generator. J Comp Phys 190:478–500 ?? ??
Aragón JL, Barrio RA, Varea C (1999) Turing patterns on a sphere. Phys Rev E 60:4588–4592 ??
Barrio RA, Maini PK, Padilla P, Plaza RG, Sánchez-Garduno F (2004) The effect of growth and curvature on pattern formation. J Dyn Diff Equ 4:1093–1121 ??
Chaplain MAJ, Ganesh M, Graham IG (2001) Spatio-temporal pattern formation on spherical surfaces: numerical simulation and application to solid tumour growth. J Math Biol 42:387–423 ??
Dziuk G, Elliott CM (2007) Finite elements on evolving surfaces. IMA J Num Anal 27:262–292 ??
Dziuk G, Elliott CM (2007) Surface finite elements for parabolic equations. J Comp Math 25:430–439 ??
Barrett JW, Garcke H, Nürnberg (2008) Parametric approximation of Willmore flow and related geometric evolution equations. SIAM J Sci Comput 31:225–253 ??
Barreira R (2009) Numerical solution of non-linear partial differential equations on triangulated surfaces, D.Phil Thesis, University of Sussex ??
Elliott CM, Stinner B (2010) Modeling and computation of two phase geometric biomembranes using surface finite elements. J Comp Phys 229:6585–6612 ??
Prigogine I, Lefever R (1968) Symmetry breaking instabilities in dissipative systems. II. J Chem Phys 48:1695–1700 ??
Schnakenberg J (1979) Simple chemical reaction systems with limit cycle behaviour. J Theor Biol 81:389–400 ??
Golub GH, Van Loan CF (1996) Matrix Computations. JHU Press, Baltimore ??
Greer J, Bertozzi AL, Sapiro G (2006) Fourth order partial differential equations on general geometries. J Comput Phys 216:216–246 ??
Calhoun DA, Helzel C (2009) A finite volume method for solving parabolic equations on logically cartesian curved surface meshes. SIAM J Sci Comput 6:4066–4099 ??
Dziuk G, Elliott CM (2010) An Eulerian approach to transport and diffusion on evolving surfaces. Comput Vis Sci 13:17–28 ??
Lefevre J, Mangin J-F (2010) A reaction-diffusion model of the human brain development. PLoS Comput Biol 6:e1000749 ??
Elliott CM, Stinner B, Styles VM (2010) Numerical computation of advection and diffusion on evolv-ing diffuse interfaces. IMA J Num Anal. Advance Access published on 11 May 2010. ??
Deckelnick KP, Dziuk G, Elliott CM (2005) Computation of Geometric PDEs and Mean Curvature Flow. Acta Numerica 14:139–232 ??