GTBIM：一种基于GraphSage与Transformer的多层网络动态层加权影响最大化方法

《Chaos, Solitons & Fractals》：GTBIM: A GraphSage–transformer-based approach for influence maximization with dynamic layer weighting in multilayer networks

　　安德鲁·莫罗佐夫|卡利安·曼纳|马莱·巴纳吉|谢尔盖·彼得罗夫斯基
英国莱斯特大学计算与数学科学学院，LE1 7RH

摘要
对骚乱和抗议活动的数学建模正逐渐成为理解社会动荡动态的强大工具，其最终目标是确保人类社会的可持续发展。本文采用反应-扩散建模框架来分析社会抗议活动的时空模式。我们的模型包含两个变量：抗议者数量以及抗议活动所造成的累积破坏程度。该系统既通过解析方法进行研究，也通过一维和二维物理空间中的大规模数值模拟进行探讨。抗议模型呈现出多种动态状态，包括规则模式和混沌模式。特别地，我们发现了一种有趣的现象：在相同的模型参数下，两种具有不同空间结构的混沌吸引子可以共存。我们还揭示了模型参数空间的复杂结构，确定了不同动态状态可能出现的参数范围。

1. 引言
长期以来，人们对通过街头行动表现出来的社会现象——如示威、游行、街头犯罪、骚乱（常伴随抢劫行为）等——一直很感兴趣[1]、[2]、[3]。无论这些行为是暴力的还是非暴力的，它们都可能由各种原因引发，并可能因不同事件而触发。例如，街头行动常常是公众对不受欢迎政府决策的反对表现[4]、[5]、[6]，也可能是对严重社会不公的反应[7]，或是由于极度不公平的社会制度所致（如[8]、[9]）。尽管街头行动的具体原因各不相同，但它们往往并非仅仅由某个特定时间和地点的单一事件引起。根据公众失望程度的不同，街头行动可能持续数天甚至数周，有时甚至数月或数年。在某些情况下，这类行为可能会变得普遍化（如巴西或美国大型城市中的街头犯罪），或者促使国家政治格局发生变化，比如催生新的政党或社会运动[2]。最终，它们甚至可能引发革命，从而导致整个政治体系的彻底变革。同样，街头行动的空间范围也可能相当广泛，覆盖整个大城市甚至全国（如[1]、[10]）。

尤其是在持续时间较长的街头行动中，预测其持续时间和强度是一个重要且极具实践意义的难题（例如，在特定日期参与行动的人数，参见[11]），因为这两个因素决定了对公共财产、私人财产以及整体经济的破坏程度。行动的持续时间和相应的时间模式各不相同，从几天到几年不等，参与者数量也可能从单调下降[12]、[13]到持续波动[14]、[15]不等。影响街头行动随时间变化的因素有很多，尤其是社会规范的差异性[16]、[17]、行为反应的变动性[18]、警方执法的影响[3]、[13]，以及昼夜规律性（在街头犯罪和抢劫活动中尤为明显，参见[14]、[15]）。

虽然街头行动的时间变化是一种相对容易被理解的常见现象，但其空间强度的变化虽然经常被观察到，见图1，却不太明显，因此更具研究价值。已有学者提出了一些解释这一现象的因素或机制。其中之一是骚乱者与警察之间的“猫捉老鼠”游戏，即闹事群体不得不根据警方的行动不断改变位置（参见[19]）。另一种更为常见的机制则是空间异质性的影响，这种影响可以以多种方式体现出来。例如，在较大的空间尺度上（如一个大区域或全国范围内），人口密度较高的地区更可能出现更多抗议者[10]。在较小的空间尺度上（如大型城市内），抗议者或骚乱者聚集的地方可能与交通便利性有关（如城市交通系统的布局，参见[20]），或者与容易遭到抢劫的目标物如商店的分布有关[21]。

因此，毫无疑问，空间异质性会在很大程度上影响街头行动的空间分布模式。不过，我们在这里假设，可能存在另一种不太明显的通用机制，也能导致抗议者在空间上的非均匀分布。实际上，当涉及大量人类时，尤其是在足够大的空间和时间尺度上，他们的行为规律可以用生态系统中的种群动态定律来描述[22]、[23]、[24]。这又使得一些已知的模式形成机制变得重要起来[25]。因此，我们可以预期，由这些机制产生的空间和时空模式，例如扩散不稳定性以及时空混沌，可能在一定程度上决定了街头行动中抗议者的分布情况。

需要指出的是，可以通过多种研究方法来分析街头行动在空间和时间上的动态特征，例如对抗议者数量数据的统计分析[11]，包括地理统计学方法[29]、[30]，以及对参与者社会特征的回归分析[31]。近年来，越来越多的研究者开始运用机器学习技术来研究空间中的复杂犯罪网络[32]、[33]。理论研究则通常采用数学建模和非线性动力学方法（如[19]、[34]、[35]）。贝雷斯蒂基及其合作者可能是最早使用反应-扩散框架来描述抗议和骚乱的时空动态的人[36]、[37]，不过他们的研究主要关注骚乱通过“骚乱波”在空间中的传播速度，因此没有涉及空间分布模式的问题（尽管在他们的一些数据中可以观察到这种现象，如[37]中的图2d和图3）。这一研究方向在[38]中得到了进一步探索（也可参见[39]），在那里人们发现了抗议活动的自崩溃波现象。最近还有文献对利用物理学方法研究社会现象——包括犯罪的空间传播——进行了深入综述，见[40]、[41]。

下载：高分辨率图像（722KB）
下载：全尺寸图像

图1. 来自历史数据的各类社会动荡（A)-(E)以及犯罪活动（F）的空间分布实例（详情见下文）。这些示例展示了不同空间尺度上零星分布的热点区域。所有案例的联合图形左侧都标有南北方向的箭头。(A) 2019年10月18日至22日，智利圣地亚哥发生动荡期间的骚乱强度分布（改编自[20]）。(B) 2018年11月12日至2019年3月2日期间，法国“黄背心”运动的地理分布（改编自[10]）。(C) 2023年6月6日，法国巴黎“纳赫尔·梅尔祖克”骚乱活动的空间分布。橙色表示媒体报道的骚乱事件，蓝色则表示社交媒体上报道的夜间骚乱事件（改编自[26]）。(D) 2022年，伊朗德黑兰城市中街头抗议活动的空间分布（改编自[27]）。(E) 2011年8月，英国大伦敦地区骚乱期间与骚乱相关的犯罪活动空间分布（改编自[21]）。(F) 2001年6月，美国加利福尼亚州长滩地区住宅盗窃案的高发区域（改编自[28]）。

在本文中，我们运用数学建模的思想和方法来研究抗议者空间分布现象。具体而言，我们考虑了一个由反应-扩散方程描述的通用模型，该模型将抗议者的局部密度与其街头行动所造成的破坏密度联系在一起。我们证明，该模型能够预测多种模式形成场景，包括图灵模式和非图灵模式，后者在空间上可能是规则的或不规则的，在时间上可能是周期性的或混沌的。有趣的是，我们的模型还揭示了两种不同的时空混沌吸引子共存的可能性——这是一种在文献中很少被提及的特殊现象。

本文的结构如下：下一节将详细介绍所采用的数学模型。第3节则研究模型的简化版（无空间维度版本），以便在参数平面上揭示其完整的分岔结构。第4节则针对具有空间显式性的系统展开研究，展示在不同模式形成场景下所出现的丰富多样的时空模式，涵盖了1维和2维空间情况（分别为4.1节一维空间情况，4.2节二维空间模式）。第4.3节则通过弱非线性分析方法对参数平面的精细结构进行研究。第5节将对我们的研究成果进行总结与讨论。我们认为，这些研究为理解街头行动在时间和空间上的动态特征提供了全新的视角，尤其有助于提高对其可预测性及管理方法的认知。

2. 模型方程
我们的数学模型是在以往街头抗议模型基础上的扩展（参见[12]、[38]）。该系统包含两个组成部分：积极参与街头行动的人数U（即抗议者数量），以及街头行动所造成的累积破坏量V（包括直接和间接破坏）。直接破坏的例子包括故意破坏财物、抢劫商店、焚烧汽车等，而间接破坏则会影响到经济和城市正常运转（如公共交通中断、当地商店关闭等）。我们假设，破坏产生的速度与参与街头行动的人数成正比。与一些先前的研究[12]、[38]不同，这里我们假设政府会不断修复这些破坏，模型中修复的速度与破坏量V成正比。模型中的U和V这两个分量都会按照随机布朗运动在空间中扩散，这一过程在方程中用经典的扩散项来表示。需要注意的是，该模型还可以进一步扩展，引入更复杂的空间扩散方式，例如趋化作用或/和交叉扩散[28]，但我们在此不予考虑。关于U和V的时空动态，相应的模型方程如下：
(1) ?U/?t = DUΔU + F(U) ? G(U,V)，
(2) ?V/?t = DVΔV + βU ? γV，
其中Δ表示拉普拉斯（扩散）算子，DU和DV分别是对应的扩散系数。F(U)表示新抗议者的加入率，而G(U,V)则表示由于抗议者失望而退出运动的比率。在V的方程中，βU和γV分别表示抗议者造成的破坏的产生和消除过程，β和γ均为正系数。沿袭以往的研究[12]、[38]，我们对函数F(U)和G(U,V)给出如下参数化形式：
(3) εF(U) = εU + aU2h2 + U2，G(U,V) = cUV，
其中ε、a、h、c均为正系数。从F(U)的函数形式可以看出，当U值较低或较高时，加入率接近线性函数，而在中等抗议者数量时则呈现S形增长趋势。这里我们假设，累积的破坏程度是抑制抗议者加入的主要因素，因此没有在F(U)中明确考虑大U值时的饱和效应。此外，与以往研究不同，为简化问题，我们在这里省略了加入率中的常数项，假设F(0)=0（参见[12]、[38]）。需要说明的是，我们总是可以通过重新调整破坏量的单位，使得c≡1，因此，我们进一步假设G(U,V)=UV，而不失一般性。同样，通过重新调整空间尺度，我们也可以假设破坏的扩散系数DV=1。这样，DU就可以理解为抗议者扩散系数与所造成破坏扩散系数之间的比值。

将F(U)和G(U,V)的参数化形式代入方程(4)–(5)后，可以得到以下系统：
(4) ε?U/?t = DUΔU + εU ? UV + aU2h2 + U2，
(5) ?V/?t = ΔV + βU ? γV。
在上述系统中，我们假设初始条件为非负值（初始条件选择的相关说明见下文），边界条件则为零通量。空间域定义为Ω?Rn（其中n=1,2），其边界?Ω为分段光滑的边界。在大多数模拟中，我们考虑的是矩形空间域（针对二维空间情况）。另一个关键问题在于如何合理选择模型参数?a、h、β、γ、?、DU的取值范围。必须指出，要精确确定这些数值是非常困难的（甚至是不可能的）。在本文中，我们对部分参数采用粗略估计，而对其他参数则在其较大范围内变动，以此验证结果的敏感性。尤其是在大多数模拟中，我们固定?a=20000、h=110、β=20、?=1，然后将γ和DU作为主要控制参数进行变化。这些参数值的选取是基于我们对该无扩散模型的数值模拟结果，该模拟得出了2018–2020年法国黄背心运动中抗议者数量峰值及其抗议周期持续时间等具有合理特征的动态模式（为此我们使用了非空间系统，相关图表见补充材料SM1）。在此情况下，U的单位为数千名抗议者。对于a和h的值，我们可以设定一个约束条件，即a/h2的量级应接近1（例如，若该比值过小，则实际中的加速效应可忽略不计）。我们还认为参数h应在100左右（以千人为单位）；这一参数决定了抗议者数量的一个临界值，超过该值后招募速度的上升趋势会停止，F(U)则开始呈现下降趋势。最后需要强调的是，我们还简要验证了结果对参数β的依赖性。

为对公式(4)–(5)进行数值模拟，我们采用显式计算方案对模型进行离散化处理——扩散部分使用中心差分法，反应部分使用欧拉法，并设置无通量边界条件。此外，我们还适当选择了空间和时间网格的大小，以避免出现数值伪影。为确保所选网格足够精细，我们将时间和空间步长减半，检查使用更细网格的模拟是否会影响结果。在我们的模拟中，无论是于一维还是二维空间情形下，我们都考虑以下两种初始条件：(i) U和V的初始空间分布是通过对方程组共存稳态E?=(U?,V?)进行小幅随机扰动得到的（该系统始终存在唯一的正共存状态，详见下文）。这些随机扰动源自不相关的高斯（白）噪声，且每个空间网格中的扰动是相互独立的。(ii) U和V的初始空间分布为U和V的静态非均匀分布，这类分布通常在DU值较低时出现，被称为图灵模式（详见后续章节）。为获得图灵模式，我们先从类型(i)的初始条件出发，然后进行长时间模拟以得到渐近解（关于构建静态图灵模式所用的参数的更多细节见第4节）。

使用MALTAB平台编写的相应计算代码可通过https://github.com/KalyanManna1/Sample-codes-for-social-protest-modelling获取。

3. 均质（非空间）情形下的抗议活动动态

在探讨包含扩散项的完整空间模型(4)–(5)之前，我们先简单研究一下抗议者和破坏行为的空间分布均质的均质情形。这种情形出现在空间位移的尺度与区域的实际（地理）尺寸相当时。我们将这种情形称为非空间模型。

对于任何参数组合，非空间模型都存在一个平凡平衡点E0=(0,0)，表示没有抗议活动及与之相关的破坏行为。该模型还可能存在一个非平凡平衡点E?=(U?,V?)，此时U?和V?的值由以下代数方程组确定：(6)??βU?3??γU?2?(aγ?βh2)U???γh2=0，(7)V?=βγU?。根据笛卡尔符号法则，当aγ>βh2时，方程(6)有唯一正解；而当aγ< /><βh2时则可能还存在两个额外的共存平衡点。需注意，该模型不允许出现抗议者数量为正但破坏程度为零的平衡状态，这是因为该模型假设破坏的产生与现有破坏程度无关。 平衡点的（局部）稳定性可通过在这些状态下计算的雅可比矩阵来确定。对于平凡平衡点e0=(0,0)，我们有?J(E0)=?0β?γ。可以看出，J(E0)的特征值为?λ1=?>0和λ2=?γ<0，因此E0是一个鞍点，即始终处于不稳定状态。而在共存平衡点E?=(U?,V?)处，雅可比矩阵的表达式为J(E?)=aU?(h2?U?2)(h2+U?2)2?U?β?γ。对应的特征方程为λ2?trace(J(E?))λ+det(J(E?))=0，其中trace(J(E?))=aU?(h2?U?2)(h2+U?2)2?γ，det(J(E?))=βU?+aγU?(U?2?h2)(h2+U?2)2。因此，当trace(J(E?))<0且det(J(E?))>0时，非平凡平衡点E?=(U?,V?)会呈现局部渐近稳定状态[42]。

下载：下载高分辨率图像（641KB）
下载：下载全尺寸图像

图2. 非空间抗议模型中可能出现的不同动态模式，以及这些模式与模型关键参数的关联。(A) 非空间模型在(γ,β)平面上的双参数分岔图。图中编号区域含义如下：区域1——全局稳定平衡态；区域2——全局吸引的极限环；区域3——两个共存的吸引子：一个局部稳定的极限环和一个局部稳定的平衡态。蓝色曲线为霍普夫分岔曲线，灰色虚线曲线为周期极限点，即不稳定周期环与稳定周期环的交汇处。BT点和GH点分别代表博格丹诺夫-塔肯斯分岔点和广义霍普夫分岔点。图(B)–(D)展示了分岔图区域1–3对应的相图示例。在这些图中，绿色点和红色点分别表示稳定平衡态和不稳定平衡态，绿色曲线和红色曲线则分别表示稳定极限环和不稳定极限环。E?为共存平衡态，E0为平凡平衡态。参数值为：(B) β=20，γ=22；(C) β=20，γ=18；(D) β=20，γ=20。其他参数值为??=1，a=20000，h=110。

下载：下载高分辨率图像（312KB）
下载：下载全尺寸图像

图3. 在不存在空间扩散（A，即非空间模型）以及存在一定范围γ值（即抗议者造成的破坏的人均清除率）的空间扩散（B）条件下，抗议模型所处的稳定性特征。 (A) 分岔图显示了非空间模型中平衡态E?的稳定性随γ变化而发生的变化规律。绿色曲线表示模型的吸引子（即稳定状态或极限环），红色曲线则表示排斥子。(B) 在(γ,DU)参数空间中观察到的空间模型中的图灵不稳定性。图灵分岔曲线用实心黑色曲线表示。其他参数与图(A)中的相同。黑色虚线和点划线垂直直线分别对应图(A)中的霍普夫分岔阈值和鞍结分岔阈值。颜色含义如下：绿色——稳定区域；粉色——图灵区域；灰色——霍普夫区域；黄色——霍普夫-图灵区域；青色——具有双稳性的稳定区域；白色——具有双稳性的图灵区域。其他参数值为??=1，a=20000，h=110，β=20。（关于图中颜色标注的含义，可参阅本文的网络版本。）

当改变某个模型参数时，非空间模型会在平衡点E?=(U?,V?)处发生霍普夫分岔。在本文中，当满足以下条件时，随着参数γ的变化就会发生霍普夫分岔[42]：trace(J(E?))|γ=γH=0，d[trace(J(E?))]dγ/γ=γH≠0，且det(J(E?))|γ=γH>0。我们通过数值方法验证了上述条件，并在参数空间中绘制出了霍普夫分岔曲线。由于模型参数的不同，第一个李雅普诺夫系数的符号可能为正也可能为负，因此霍普夫分岔可能是临界分岔也可能是亚临界分岔（详见下文）。

我们通过固定其他模型参数，在(γ,β)参数空间中为非空间模型绘制了一系列分岔图。图2A展示了一个此类分岔图的示例。与该图对应的相图则呈现于图2B–D中。从图2A可以看出，对于给定的??、a和h值，可能存在三种不同的参数区域。在区域1中，共存平衡态是全局吸引子：从所有正初始条件出发，系统轨迹都会趋向于E?状态。在参数区域2中，全局吸引子为极限环，因此抗议者数量以及相应的破坏程度会随时间周期性波动（此类时间动态的示例见SM1）。最后，在参数位于区域3时，我们会观察到双稳性现象：根据初始条件的不同，系统轨迹最终可能落在平衡态E?处，也可能落在极限环上——这两个吸引子的吸引域由那个不稳定的极限环决定（见图2D中的红色曲线）。

我们还展示了仅改变γ值时上述(γ,β)分岔图的横截面示意图（见图3A）。当β保持为20时，较小的γ值（即破坏清除速率较慢）会导致系统处于稳定平衡态。随着γ值逐渐增大，会出现超临界霍普夫分岔（在γ=γH(1)处），此时U?的振荡幅度会迅速增大，远远超过霍普夫分岔前稳定平衡态下的U?值。继续增大γ值后，平衡态会重新趋于稳定，但这是通过亚临界霍普夫分岔实现的（在γ=γH(2)处），因此仍有可能出现大幅度的周期性波动，具体取决于初始条件。最后，当破坏清除速率进一步增大时，不稳定的极限环与稳定的极限环会在γ=γSN处发生碰撞（即极限环的鞍结分岔）。值得注意的是，尽管在大γ值和小γ值情况下系统的行为在定性上相似（即当γ>γSN和γ<γh(1)时），但在高破坏清除速率下，平衡态下的抗议者数量会显著增多。这说明快速的破坏清除会导致持续性的抗议活动或骚乱，参与人数也会大幅增加。

4. 抗议活动的空间动态特性

4.1 一维空间情形

我们首先通过研究一维物理空间来探讨抗议活动的空间动态特性。这种情形可用于描述那些在某一空间方向上延伸、而在垂直方向上较窄的城市地区的抗议活动。同时需要强调的是，了解一维情形下的空间动态特性对于更好地理解同一模型在二维情况下的动态模式至关重要。我们首先对均匀静态状态e?=(u?,v?)进行了局部稳定性分析。通过标准方法（详细内容见附录a），我们发现当存在扩散作用时，该状态可能会失去稳定性，从而出现图灵模式。图3b总结了随着参数γ（即破坏清除速率）变化时图灵不稳定性方面的研究结果。在所研究的γ参数范围内，存在六个不同的区域：稳定区域（存在扩散时e?仍保持局部稳定）；图灵区域（存在扩散时e?失去稳定性）；霍普夫区域（在没有扩散的系统中e?本身就不稳定）；霍普夫-图灵区域（同时存在时间不稳定性和空间不稳定性）；具有双稳性的稳定区域（存在扩散时e?仍保持局部稳定，但均匀系统中还存在另一个吸引子，即极限环）；以及具有双稳性的图灵区域（存在扩散时e?失去稳定性，且均匀系统中也存在另一个吸引子，即极限环）。图3b中不同区域的局部不稳定性特征通过散射关系图进行展示，这些图分别见于sm2中。需要特别指出的是，只有当破坏传播的速度远快于抗议者移动速度时（即du?dv=1），才会出现图灵不稳定性。

随后，我们对系统进行了更为全面的数值研究，旨在寻找远离均匀状态e?的时空动态模式。我们的研究结果（结合局部稳定性分析）汇总在图4的(γ,du)图中，图中不同的符号代表了在一维空间模型中出现的不同动态模式。相关空间模式的示例可见图5、图6以及补充材料sm3-sm4中。在抗议活动模型中观察到的所有空间非均匀模式都可以归为三类：静态的（且空间上具有规则性）；非静态的周期性模式；以及时空混沌模式。

下载：下载高分辨率图像（916kb）
下载：下载全尺寸图像

图4. 4. 抗议活动的空间动态特性 4.1 一维空间情形 我们首先通过研究一维物理空间来探讨抗议活动的空间动态特性。这种情形可用于描述那些在某一空间方向上延伸、而在垂直方向上较窄的城市地区的抗议活动。同时需要强调的是，了解一维情形下的空间动态特性对于更好地理解同一模型在二维情况下的动态模式至关重要。我们首先对均匀静态状态e?=(U?,V?)进行了局部稳定性分析。通过标准方法（详细内容见附录A），我们发现当存在扩散作用时，该状态可能会失去稳定性，从而出现图灵模式。图3B总结了随着参数γ（即破坏清除速率）变化时图灵不稳定性方面的研究结果。在所研究的γ参数范围内，存在六个不同的区域：稳定区域（存在扩散时E?仍保持局部稳定）；图灵区域（存在扩散时E?失去稳定性）；霍普夫区域（在没有扩散的系统中E?本身就不稳定）；霍普夫-图灵区域（同时存在时间不稳定性和空间不稳定性）；具有双稳性的稳定区域（存在扩散时E?仍保持局部稳定，但均匀系统中还存在另一个吸引子，即极限环）；以及具有双稳性的图灵区域（存在扩散时E?失去稳定性，且均匀系统中也存在另一个吸引子，即极限环）。图3B中不同区域的局部不稳定性特征通过散射关系图进行展示，这些图分别见于SM2中。需要特别指出的是，只有当破坏传播的速度远快于抗议者移动速度时（即DU?DV=1），才会出现图灵不稳定性。 随后，我们对系统进行了更为全面的数值研究，旨在寻找远离均匀状态e?的时空动态模式。我们的研究结果（结合局部稳定性分析）汇总在图4的(γ,du)图中，图中不同的符号代表了在一维空间模型中出现的不同动态模式。相关空间模式的示例可见图5、图6以及补充材料sm3-sm4中。在抗议活动模型中观察到的所有空间非均匀模式都可以归为三类：静态的（且空间上具有规则性）；非静态的周期性模式；以及时空混沌模式。 下载：下载高分辨率图像（916kb） 下载：下载全尺寸图像>
4. 抗议活动的空间动态特性

4.1 一维空间情形

我们首先通过研究一维物理空间来探讨抗议活动的空间动态特性。这种情形可用于描述那些在某一空间方向上延伸、而在垂直方向上较窄的城市地区的抗议活动。同时需要强调的是，了解一维情形下的空间动态特性对于更好地理解同一模型在二维情况下的动态模式至关重要。我们首先对均匀静态状态e?=(u?,v?)进行了局部稳定性分析。通过标准方法（详细内容见附录a），我们发现当存在扩散作用时，该状态可能会失去稳定性，从而出现图灵模式。图3b总结了随着参数γ（即破坏清除速率）变化时图灵不稳定性方面的研究结果。在所研究的γ参数范围内，存在六个不同的区域：稳定区域（存在扩散时e?仍保持局部稳定）；图灵区域（存在扩散时e?失去稳定性）；霍普夫区域（在没有扩散的系统中e?本身就不稳定）；霍普夫-图灵区域（同时存在时间不稳定性和空间不稳定性）；具有双稳性的稳定区域（存在扩散时e?仍保持局部稳定，但均匀系统中还存在另一个吸引子，即极限环）；以及具有双稳性的图灵区域（存在扩散时e?失去稳定性，且均匀系统中也存在另一个吸引子，即极限环）。图3b中不同区域的局部不稳定性特征通过散射关系图进行展示，这些图分别见于sm2中。需要特别指出的是，只有当破坏传播的速度远快于抗议者移动速度时（即du?dv=1），才会出现图灵不稳定性。

随后，我们对系统进行了更为全面的数值研究，旨在寻找远离均匀状态e?的时空动态模式。我们的研究结果（结合局部稳定性分析）汇总在图4的(γ,du)图中，图中不同的符号代表了在一维空间模型中出现的不同动态模式。相关空间模式的示例可见图5、图6以及补充材料sm3-sm4中。在抗议活动模型中观察到的所有空间非均匀模式都可以归为三类：静态的（且空间上具有规则性）；非静态的周期性模式；以及时空混沌模式。

下载：下载高分辨率图像（916kb）
下载：下载全尺寸图像

图4.>图5. 一维抗议模型（4）–（5）中出现的动态（即非稳态）模式示例。这些示例包括混沌和规则的动力学模式（详见下文）。左列展示了时空混沌的示例（DU=0.1）；右列显示了规则周期性模式（DU=0.15）；两幅图的其余参数相同：a=20000，h=110，β=20，γ=10。这两种模式都对应于图4(A)中的图谱中的图灵区域。在两列中，上排表示时间范围t∈[10000,10050]内抗议人群规模U的时空动态；中间排展示了t=10000时各组成群体（U和V）的空间分布；下排则绘制了时间区间t∈[8000,10000]内空间平均密度〈U〉和〈V〉的相图。在整个γ范围内，当DU值较小时会出现稳定（且空间上规则）的模式。对于这类模式，U的空间分布表现为局部出现大量抗议者的尖锐峰值，而这些峰值之间则是骚乱强度较低的区域。有趣的是，U和V的最大值可以比均匀稳定状态E?的数值高得多。因此，模型中引入扩散效应会在某些地理位置形成骚乱的“热点”。在完全均匀的系统中，如此大量的抗议者是难以被察觉的。下载：下载高分辨率图像（761KB）下载：下载全尺寸图像

图6. 一维抗议空间模型中混沌多稳定性的示例：从不同的初始条件出发，系统会趋向于两个共存的混沌吸引子之一。上排：从两种不同初始条件出发，空间平均密度〈U〉和〈V〉的相图。在(A)-(B)两幅图中，黑色曲线对应于从γ=0.05时得到的渐近非均匀稳定模式开始的模拟，而蓝色曲线则是从略微扰动的均匀稳定状态（U?,V?）（在两幅图中均用实心红色圆圈标出）开始的模拟。A幅图中的相图包含了瞬态动态，即t∈[0,10000]；B幅图则仅展示了最终的混沌吸引子，省略了初始的瞬态动态，时间为t∈[8000,10000]。下排：从不同初始条件出发，时间范围t∈[10000,10050]内抗议人群规模U的时空动态模拟结果：C幅图和D幅图分别对应于(A)-(B)幅图中的蓝色曲线和黑色曲线。模型参数为DU=0.2，??=1，a=20000，h=110，β=20，γ=13。空间域的大小为Ω=[0,10]。（关于本图例中颜色含义的解释，可参阅本文的网络版本。）模型中的非稳态（动态）模式可以是周期性的，也可以是混沌的。图5给出了周期性和混沌性时空模式的两个示例。我们还在同一图中展示了空间平均密度〈U〉和〈V〉的相图（见下排）。需注意，我们判断模式具有混沌性质的依据是相近轨迹的指数级发散（为简化起见，此处未展示相应图表）。我们发现，从周期性动态向时空混沌过渡的典型途径是通过周期倍增分岔。关于通过周期倍增实现混沌过渡的示例，可见于补充材料SM4中。我们的模拟表明，即使空间分布相对简单，例如仅由一个区域构成，平均密度的复杂混沌动态也有可能出现（详情见SM4中的图S4）。需要指出的是，此前在其他反应-扩散模型中也已有研究报道过物种空间分布规则的系统存在时间混沌现象[43]。我们还发现，混沌模式（无论是空间分布规则还是不规则）大多出现在图4中的霍普夫域或图灵-霍普夫域内，并不需要满足图灵条件DU?DV=1。具有扩散效应的骚乱模型的一个特别有趣的特点是，不同的时空混沌吸引子可以共存，即系统拥有多个混沌吸引子。图6展示了一个示例，根据初始条件的不同（详情见图注），轨迹会趋向于两个吸引子之一。其中一个吸引子的特征是抗议者数量和破坏程度更高（在〈U〉 - 〈V〉平面上，它位于由实心圆圈标出的不稳定平衡点附近）。另一个吸引子的U和V值则较低。两个吸引子对应的抗议者空间分布在定性上相似，但在空间规则性程度上存在差异，如图6 C、D所示。需要注意的是，对应的无空间结构系统在给定参数下并不表现出多稳定性。我们还发现，趋向于某个混沌吸引子可能会引发漫长的瞬态阶段，如图6A中的黑色轨迹所示，而对于另一个吸引子，轨迹则会很快到达其附近区域（蓝色轨迹）。因此，在合理的时间内（例如30–50年），只能达到一个混沌吸引子。另一种多稳定性的例子是，在相同的参数设定下，不同的稳定状态可以共存，这体现在图4B中，它是图4A主图的放大部分。特别是在图灵稳定区域（即图灵曲线之上），可能出现稳定的非均匀（规则）模式。因此，在这种情况下，均匀状态E?仅具有局部稳定性：对该平衡点的超临界扰动会导致其不稳定，最终的吸引子将是模型组分在空间上交替出现高密度和低密度抗议者的规则周期性空间分布（参见SM3）。需注意，图4B中的图谱是通过第2节中描述的两种不同类型的初始条件构建的。最后，我们研究了空间域大小对骚乱时空动态的影响。我们发现，一般来说，空间域的物理尺寸增大会导致空间不规则性增加，进而促进时空混沌的产生。这一观察结果与早期关于低维系统中混沌现象的一些研究[44]、[45]相符。我们发现在更大的空间环境中，出现纯时间周期性模式的可能性较小。图7给出了一个示例，其中空间域的大小在保持图6相同参数的情况下增加了10倍。结果，图6中较小空间域所观察到的周期性模式（右列）消失了，取而代之的是混沌动态。重要的是，在较大的空间域中，空间平均值〈U〉和〈V〉的振荡幅度有所减小。从图7可以看出，振荡发生在〈U〉 - 〈V〉平面上的不稳定平衡点E?周围。这可以通过以下事实来解释：整个空间被分割成多个独立的子空间，每个子空间内都存在大幅度的混沌振荡。对大量独立混沌振荡器进行平均后，整体方差降低，即偏离E?给出的平均值的程度减小。下载：下载高分辨率图像（1010KB）下载：下载全尺寸图像

图7. 空间域大小对一维模型中抗议模式的影响。该图展示了空间域为Ω=[0,100]时的抗议动态，其大小是图5中所展示的空间的十倍。上排表示时间范围t∈[10000,10050]内抗议人群规模U的时空动态；中间排展示了t=10000时各组成群体（U和V）的空间分布；下排则展示了时间区间t∈[8000,10000]内空间平均密度〈U〉和〈V〉的相图。所有模拟都是在使用与图5相同的参数情况下进行的（空间域大小除外）。下载：下载高分辨率图像（1MB）下载：下载全尺寸图像

图8. 二维物理空间中抗议空间模型（4）–（5）可能出现的状态的参数结构与分类。(A) 为空间域Ω=[0,10]×[0,10]构建的(γ,DU)参数图。(B) 为(A)图中东南角的放大视图。其他参数为??=1，a=20000，h=110，β=20。黑色虚线和点划线垂直直线分别表示霍普夫分岔阈值和鞍点分岔阈值。不同的模式由以下彩色符号表示：- 稳定的空间均匀分布；- 热点；- 冷点；- 热点与条带的混合；- 冷点与条带的混合；- 迷宫式结构；- 动态热点；- 动态热点与条带的混合；- 动态迷宫式结构；- 相互作用的螺旋混沌；- 其他混沌模式；- 动态冷点与条带的混合；- 动态冷点；- 周期性模式（包括空间均匀和非均匀类型）。4.2 二维空间模式图8展示了二维抗议空间模型的参数图。可以看出，增加第二个空间维度会带来多种在一维情况下无法观察到的新模式。另一方面，两种空间情形下参数图的整体结构仍有一定的相似性（参见图4）。我们可以将所有发现的二维非均匀模式大致分为四类：稳定的空间结构；接近稳定结构的动态模式；空间不规则性高的混沌动态；以及周期性动态。二维模型中大多数模式的示例见SM5（见图S5）。需要指出的是，在上述分类中，我们有意没有明确界定那些接近稳定结构但具有混沌性质的动态模式与具有高空间不规则性的“真正”混沌动态之间的确切界限。该模型还存在漫长的瞬态阶段，在这些阶段中，趋向于某个吸引子（例如非均匀稳定模式）的过程可能持续很长时间，远远超过实际骚乱动态的时间尺度。我们首先介绍在DU值较低时观察到的二维稳定模式。这类模式共有五种，其中四种如图9所示（该图展示了初始瞬态消失后的最终抗议者数量空间分布）。这五种类型分别是：热点、冷点、热点与条带的混合、冷点与条带的混合，以及迷宫式结构。从骚乱建模的角度来看，热点模式描述的是骚乱活动仅在某些地点发生，而大部分区域抗议者数量较少的情况。冷点模式则呈现出相反的空间分布，即大多数地区都被骚乱者占据。空间分布可能是条带和斑点（冷点或热点）的复杂混合体，尤其是环状条带，一些示例见图9 C和SM5中的图S5。下载：下载高分辨率图像（2MB）下载：下载全尺寸图像

图9. 二维抗议空间模型中观察到的不同类型稳定模式的示例。所有子图均展示了t=1000时抗议者数量（U）的分布。(A) 热点（在图8中用标记）：γ=10，DU=0.1；(B) 冷点（在图8中用标记）：γ=21，DU=0.02；(C) 热点斑块与条带的混合（在图8中用标记）：γ=10，DU=0.02；(D) 迷宫式结构（在图8中用标记）：γ=16，DU=0.1。其他参数值分别为??=1，a=20000，h=110，β=20。当扩散系数DU从较低值增加到较高值时，静态的非均匀模式会转变为动态模式。这类动态模式的空间结构与静态模式类似（例如热点），但会随时间发生变化。在SM5（图S5）中展示了接近静态非均匀结构的动态模式示例。而空间平均值〈U〉和〈V〉的对应图表则显示，这些模式可能是混沌的，也可能是近似周期性的，相关内容见SM6（图S6）。我们还发现，上述某些动态模式可能是长期暂态现象，逐渐逼近静态吸引子（例如图S6(e)所示）：数值模拟表明这类模式的暂态时间可能非常长，可达104年或更长的弛豫时间。下载：高分辨率图像（868KB）下载：全尺寸图像图10. 在二维抗议活动空间模型中观察到的不同时空混沌模式示例。上方面板(A)-(C)为t=1000时抗议规模(U)的空间分布快照：(A)在γ=12.5、DU=0.3条件下，通过相互作用螺旋产生的混沌动力学（如图8中的标记所示）；(B)在γ=19.5、DU=0.5条件下，具有高度空间不规则性的混沌（如图8中的标记所示）；(C)在γ=13、DU=0.5条件下，空间不规则性较低的混沌。下方面板(D)-(F)为上方面板所示模式对应的空间平均密度〈U〉和〈V〉的时间变化情况（t∈[900,1000]）。其他参数仍为??=1，a=20000，h=110，β=20。当DU进一步增加时，模型会呈现混沌动力学特征，即空间和时间上都存在高度不规则的模式。这种情况主要出现在霍普夫分岔区和图灵-霍普夫分岔区。相应的混沌模式示例见图10A、B、D、E。尤其我们观察到了螺旋混沌，也被称为“螺旋湍流”，其混沌现象是由旋转螺旋的破裂所引发的（见图10A中的示例）。随着DU的进一步增加，混沌动力学的不规则程度会降低（见图10C、F），这是因为斑块的特征尺寸会增大，达到与空间域大小相当的水平。在霍普夫分岔区和图灵-霍普夫分岔区，DU进一步增加时，系统会出现周期性振荡，且异质性程度很低。最后，我们发现了双稳态现象，这与我们之前在一维空间模型中的研究结果一致。这一现象在图9 B中有部分展示，显示了靠近图灵曲线的分岔区域内的静态模式（为简化起见，图中未同时标出表示局部稳定均匀状态的符号）。关于系统在图灵曲线附近行为的更多信息将在下一节中介绍。我们还简要验证了空间域大小对系统空间动力学的影响，发现与一维情况类似，空间域大小的增加通常会促进混沌动力学的发展（为简洁起见，此处未附相应图表）。4.3 二维模型的弱非线性分析在上一节中，我们通过数值模拟揭示了二维暴动模型中的多种空间模式。在这里，我们运用分析工具来证明某些特定类型模式的存在，尤其是那些在接近图灵不稳定曲线参数值时出现的模式。为此，我们采用弱非线性分析方法，参照[46]中的步骤推导振幅方程。我们将(4)–(5)式中的反应项分别记为?F1(U,V)=?U?UV+aU2h2+U2以及F2(U,V)=βU?γV。接着，我们在平衡点E?=(U?,V?)附近展开泰勒级数至三阶项，用偏离E?的值，即U?U?和V?V?来重新表述该系统。为简化符号，后续我们将继续使用U和V来表示对平衡点的偏离量。由此可得：(8)?U?t=DUΔU+∑m+n=13AmnUmVn,(9)?V?t=ΔV+∑m+n=13BmnUmVn，其中Amn=1m!n!?m+nF1?Um?VnE?，Bmn=1m!n!?m+nF2?Um?VnE?，且m,n≥0。通过简单计算可得A10=aU?(h2?U?2)(h2+U?2)2，A01=?U?，B10=β，B01=?γ，A20=2ah2(h2?3U?2)(h2+U?2)3，A11=?1，A02=0，B20=0，B11=0，B02=0，A30=24ah2U?(U?2?h2)(h2+U?2)4，A21=0，A12=0，A03=0，B30=0，B21=0，B12=0，B03=0。上述系统可用紧凑的向量符号W=(U,V)重新表示，即：(10)?W?t=LW+H，其中线性算子L和非线性项H定义为L=A10+DUΔA01B10B01+Δ，H=∑m+n=23AmnUmVn∑m+n=23BmnUmVn=2ah2(h2?3U?2)(h2+U?2)3U2?UV+24ah2U?(U?2?h2)(h2+U?2)4U30。设DUc代表图灵不稳定阈值，也就是图灵不稳定曲线。因此，当DU< /><μ2=0时该状态稳定，而当μ>μ2=0时则变得不稳定。(b)条带状模式对应ρ1=μΨ1≠0，ρ2=ρ3=0。当μ>μ3=h02Ψ1(Ψ2?Ψ1)2时该模式稳定，否则不稳定。(c)六边形模式表现为ρ1=ρ2=ρ3=|h0|±h02+4(Ψ1+2Ψ2)μ2(Ψ1+2Ψ2)?ρH±，且φ=0或π，只要μ>μ1=?h024(Ψ1+2Ψ2)该模式就存在。其中ρH+仅在μ<μ4=(2ψ1+ψ2)h02(ψ2?ψ1)2时稳定，而另一种解ρh?始终不稳定。(d)混合模式表现为ρ1=|h0|ψ2?ψ1，ρ2=ρ3=μ?ψ1ρ12ψ1+ψ2，且要求ψ2>Ψ1、μ>Ψ1ρ12，这类模式始终不稳定。我们运用上述标准来研究图灵不稳定曲线附近暴动模型中的静态模式。作为示例，我们考察了γ∈[19.6,22]范围内的图灵曲线附近区域（即位于霍普夫分岔区和霍普夫-图灵分岔区右侧的区域）。基于弱非线性分析的预测结果见图11。需注意，图11中的其他模型参数与图8中的主参数图相同。我们由此得到了三个区域：R1：μ< />< />< />< />< />< />< />< />< />< />< />< />< />< /><μ是无法成立的。所有区域ri对应的振幅ρ值见图8a，而图11b则展示了(γ,du)坐标系下的ri区域。我们对模型(4)–(5)进行的直接数值模拟总体上证实了通过弱非线性分析得到的预测结果。具体而言，我们发现以e?附近的小幅度扰动作为初始条件时，系统会趋于空间均匀状态；而采用正向延续技术（即从du较小时得到的图灵模式开始模拟）则会导致hπ模式（冷点模式）的出现。因此，图11b的结果与图8b在图灵分岔阈值duc附近的结论是一致的。下载：高分辨率图像（356kb）下载：全尺寸图像图11. 通过弱非线性分析方法，研究远离均匀平衡态e?的二维抗议活动模型中的空间模式。（a）为(μ,ρ)坐标系下振幅方程(26)–(28)的稳定性图，此处μ=(duc?du)/duc，duc为图灵不稳定阈值，参数设定为??=1，a=20000，h=110，β=20，γ=21，duc=0.0223。图中实线与虚线分别代表稳定和不稳定模式状态，s、hπ和h0分别表示条带状模式、φ=π时的六边形模式以及φ=0时的六边形模式（详细内容见正文）。垂直直线标出了μ1=?0.1231、μ2=0、μ3=7.6287和μ4=26.7721这几个数值。（b）为基于弱非线性分析得到的模型(4)–(5)在(γ,du)参数空间中的分岔图，参数μi的数学表达式已在正文中给出，参数域ri的定义也在图中标明。通过研究γ较小的图灵分岔区域也能得出类似结论（为简洁起见此处未展示结果）。另一方面，当γ∈[10,11.5]时（即仍处于霍普夫-图灵分岔区附近），数值模拟表明，当du接近图灵不稳定曲线时，模型在图灵分岔区域内会呈现动态模式（见图8(a)）。这是一种有些反直觉的现象，可以用无扩散系统中霍普夫分岔的向后效应来解释。这一例子也显示出弱非线性分析作为分析方法的潜在缺陷：模型中强烈的非线性效应使得我们无法使用仅在图灵分岔阈值duc附近才有效的近似方法。在其他反应-扩散系统中也报告了类似弱非线性分析的局限性[48]，上述差异的具体机制仍有待进一步研究。5. 讨论与结论在抗议或暴动等公共街头活动中，来自不同来源、针对不同事件的数据表明，抗议者在空间上的分布往往具有非均匀性，即抗议者数量较多的区域与数量较少的区域交替出现。这种非均匀性存在于不同的空间尺度上：从小范围（如大城市的某个街区）到大范围（整个国家），详见图1。通常人们认为，抗议者分布的非均匀性仅仅是空间非均匀性的结果。然而在本文中，我们提出了一种假设：即抗议活动强度的高低差异（可用抗议参与者的空间密度作为替代指标）也可能是一种自组织现象，因而可能是由系统的内部不稳定性所导致的。在我们研究了抗议活动动态的通用时空模型后，即方程(4)–(5)，证明了确实如此。所提出的框架的优势主要有以下几点。首先，我们所设计的抗议模型十分简洁，仅需要两个动态变量（抗议者数量u和累积破坏程度v），并且参数数量也很少。因此，可以说这是生成抗议/暴动活动非均匀时空动态的最简模型。模型的简洁性使我们能够全面分析各种参数值、初始条件以及不同规模的物理域ω。此外，我们还通过分析方法证明了某些非均匀状态的存在（详细内容见第4节），这也是该方法的另一个优势。最后，利用该模型还可以提出一些简单的抗议管理策略，比如通过调整破坏/垃圾清理速率来实现。我们的模型展现了丰富多样的空间模式，包括不同大小的静态或动态斑块、条带状以及迷宫状模式。这种自组织的模式形成会对系统的时空动态产生直接影响（例如通过抗议者密度的空间平均值来体现）。虽然无空间模型的动态（即不含扩散项的方程(4)–(5)）最多只能呈现周期性，但包含空间结构的模型的时空动态则可能出现准周期性和混沌行为，后者实际上属于时空混沌，因为它本质上是模式形成的结果。 通过弱非线性分析方法，研究远离均匀平衡态e?的二维抗议活动模型中的空间模式。（a）为(μ,ρ)坐标系下振幅方程(26)–(28)的稳定性图，此处μ=(DUc?DU)/DUc，DUc为图灵不稳定阈值，参数设定为??=1，a=20000，h=110，β=20，γ=21，DUc=0.0223。图中实线与虚线分别代表稳定和不稳定模式状态，S、Hπ和H0分别表示条带状模式、φ=π时的六边形模式以及φ=0时的六边形模式（详细内容见正文）。垂直直线标出了μ1=?0.1231、μ2=0、μ3=7.6287和μ4=26.7721这几个数值。（B）为基于弱非线性分析得到的模型(4)–(5)在(γ,DU)参数空间中的分岔图，参数μi的数学表达式已在正文中给出，参数域Ri的定义也在图中标明。通过研究γ较小的图灵分岔区域也能得出类似结论（为简洁起见此处未展示结果）。另一方面，当γ∈[10,11.5]时（即仍处于霍普夫-图灵分岔区附近），数值模拟表明，当DU接近图灵不稳定曲线时，模型在图灵分岔区域内会呈现动态模式（见图8(A)）。这是一种有些反直觉的现象，可以用无扩散系统中霍普夫分岔的向后效应来解释。这一例子也显示出弱非线性分析作为分析方法的潜在缺陷：模型中强烈的非线性效应使得我们无法使用仅在图灵分岔阈值DUc附近才有效的近似方法。在其他反应-扩散系统中也报告了类似弱非线性分析的局限性[48]，上述差异的具体机制仍有待进一步研究。5. 讨论与结论在抗议或暴动等公共街头活动中，来自不同来源、针对不同事件的数据表明，抗议者在空间上的分布往往具有非均匀性，即抗议者数量较多的区域与数量较少的区域交替出现。这种非均匀性存在于不同的空间尺度上：从小范围（如大城市的某个街区）到大范围（整个国家），详见图1。通常人们认为，抗议者分布的非均匀性仅仅是空间非均匀性的结果。然而在本文中，我们提出了一种假设：即抗议活动强度的高低差异（可用抗议参与者的空间密度作为替代指标）也可能是一种自组织现象，因而可能是由系统的内部不稳定性所导致的。在我们研究了抗议活动动态的通用时空模型后，即方程(4)–(5)，证明了确实如此。所提出的框架的优势主要有以下几点。首先，我们所设计的抗议模型十分简洁，仅需要两个动态变量（抗议者数量u和累积破坏程度v），并且参数数量也很少。因此，可以说这是生成抗议 暴动活动非均匀时空动态的最简模型。模型的简洁性使我们能够全面分析各种参数值、初始条件以及不同规模的物理域ω。此外，我们还通过分析方法证明了某些非均匀状态的存在（详细内容见第4节），这也是该方法的另一个优势。最后，利用该模型还可以提出一些简单的抗议管理策略，比如通过调整破坏>在一定的参数范围内，观察到了两种不同混沌吸引子的共存现象（见图6）：这一奇特现象在文献中极为罕见[49]、[50]。有趣的是，从实际应用的角度来看，我们的研究为抗议活动管理提供了一些简单的方法。事实上，从我们的结果中可以明显看出，系统的性质严重依赖于参数γ，即废物/破坏的清除速率（例如由州政府或地方当局负责），详见图3、4A和8A。一个有趣且看似违反直觉的发现是，γ的增加会导致抗议者总数的增加（见图3A）。这可被视为政府与抗议者之间一种不情愿的合作。此外，当γ从较小值增加到中等值时，系统动态的复杂性也会增加，尤其是会引发时空混沌，从而使抗议活动的走向难以预测。由此可以得出的一个直接结论是，为了降低抗议强度并使其更加有序，当局不应急于清除废物和修复损坏。需要注意的是，我们模型中将抗议者的移动视为扩散过程（参见方程式(4)–(5)中的扩散项）；因此，人类的移动被视作随机的。这一假设可能看起来不现实，因为对于任何特定个体而言，其在空间中的移动通常是有特定原因的，因此几乎不可能是随机的。然而，我们认为，由于这些原因可能多种多样且往往不为观察者所知，将一定时间间隔内的人类移动汇集起来，就可以从统计学角度对其进行描述，为便于数学建模而将其视为随机现象。这是统计物理学中用于描述复杂多粒子系统的常用方法。从这个意义上说，随机性假设更多反映的是我们对“为何移动”这一原因的了解不足，而非真正的随机性。1 我们还要提到，此前在一些有影响力的研究中也曾采用将人类移动视为随机的方法[52]。我们的研究存在一些局限性，虽然我们有意忽略了这些局限性，但未来研究应予以考虑。这些局限性包括：互动场所的物理异质性、抗议者群体内部的组织结构，以及外部干扰的影响。下面我们将简要讨论这些方面。在我们的模型中，抗议者在均匀的空间中移动。而在现实中，情况很少如此。在中等尺度上（单个小型/中型城市或大城市的某个街区），城市空间以多种不同的方式呈现结构，例如街道几何形状、贫富社区交替分布、距离交通枢纽的远近等。在大尺度上（大区域或整个国家），空间结构则由城乡区域交替分布和/或交通系统的空间结构决定[20]。因此，更现实的建模方法应考虑在具有结构的空间中，例如网络中，进行抗议者移动的模拟。这将成为未来研究的重点。我们的方法将所有参与者视为相同的个体。实际上，不同的个体抗议者在许多方面显然存在差异（有时差异很大）。因此，他们对各种情况的行为反应也可能有很大不同；这又会直接影响他们参与或退出街头活动的意愿（在当前模型中，这一点通过方程式(3)中的系数来量化）。在更现实的方法中，应根据不同的行为反应对人群进行分类[16]。尽管在非空间环境中已经有一定程度的此类研究，例如参见[12]、[13]、[53]，但在具有明确空间结构的系统中，这种行为结构会产生何种影响仍不清楚，这也应成为未来研究的重点。政府（或地方当局）通常会采取措施镇压或消除街头抗议/骚乱。这些措施可能包括部署额外的安保力量、限制公众进入特定地点、逮捕闹事者、通过社交媒体发起攻击性宣传以劝阻人们参加抗议活动等[54]。此外，抗议者和政府之间也可能展开谈判，从而缓解抗议的根源问题。对于数学模型而言，这意味着某些参数不再保持恒定，因此系统中的一些稳定模式可能不再存在，此时就需要考虑非自治的抗议模型。最后，需要指出的是，模型(4)–(5)的应用并不局限于街头抗议的动态分析。实际上，它也可以被视为一种具有自我抑制作用的种群动态生态模型。在此框架下，U表示某种物种的种群密度，而V则表示该种群对环境造成的负面影响（例如通过产生长期滞留在栖息地内的废物或毒素）。此前有一些研究探讨了自我毒性对种群动态的影响[55]、[56]；有趣的是，我们在这篇论文中获得的结果与他们的发现高度一致。特别是，一些关于产毒微生物的实验室实验表明，由于毒素的作用，物种密度可能会出现振荡动态，而在相同实验条件下，无毒素的细菌则没有出现这种现象[55]。还有类似我们的建模框架被用来描述废物产生对人类种群增长的影响[56]，该研究也表明，产生废物/毒素可能导致系统出现振荡。不过据我们所知，之前的研究尚未考虑过在具有自我毒性系统中空间的作用。可以说，我们的研究有助于更好地理解这类系统的时空动态特性。CRediT作者贡献说明：Andrew Morozov：写作——初稿撰写、可视化、验证、监督、方法论、研究设计、形式分析、概念构建。Kalyan Manna：写作——初稿撰写、验证、方法论、形式分析。Malay Banerjee：监督、方法论、形式分析、概念构建。Sergei Petrovskii：写作——初稿撰写、可视化、监督、方法论、概念构建。