John David Weissmann，苏黎世大学Axel Timmermann，釜山大学Marcia S. Ponce de Leon，苏黎世大学Christoph P.E. Zollikofer，苏黎世大学
编译来源:
Vahdati, A. R., Weissmann, J. D., Timmermann, A., de León, M. S. P., & Zollikofer, C. P. (2019). Drivers of Late Pleistocene human survival and dispersal: an agent-based modeling and machine learning approach. Quaternary Science Reviews, 221, 105867.
本文作者之一：Ali R. Vahdati
一、前言 目前几乎达成共识的是，智人（Homo sapiens）在非洲经历了漫长的进化过程，且我们的物种的非洲人群在晚更新世时期开始向全球扩散，最终实现了我们星球上所有可居住区域的人口定居（Groucutt et al.，2015；Galway-Witham, & Stringer，2018）。然而，关于人类如何在时间和空间上展开全球扩散的具体过程，仍然是一个激烈研究的问题。是否只有一次“走出非洲”事件，还是有多次扩散波，或是一个持续的人口流动？人类是否沿着特定的路线迁移，或是沿着环境梯度扩散？有哪些因素调控了非洲人口的扩散？ 多种方法已经被用来重建这些扩散事件，包括古人类学、考古学以及基因和表型多样性模式的分析，这为理解人类空间时间人口扩散的复杂性提供了重要见解。例如，非洲以外的到达时间、人口瓶颈和创始人效应，本地人口的混血和替代现象（Liu et al.，2015；Bae et al.，2017）。此外，关于智人与古代人群（如尼安德特人和丹尼索瓦人）的混血的研究，也揭示了更长时间段和更大空间范围内的基因交流（Dannemann, & Racimo，2018；Wolf, & Akey，2018）。 这一领域的研究迅速发展，几乎每天都会产生新的扩散模式和动态证据（Harvati et al.，2019）。然而，关于这些扩散背后的原因和机制，我们仍在初步探索阶段。多种因素被提出可能影响这些扩散过程，如气候变化（Abbate, & Sagri，2012），海平面变化（Armitage et al.，2011），构造作用（Bailey, & King，2011），生态变化，动物群更替率（Palombo，2013），食物资源（Dennell，2003），人口统计（Dennell，2003）以及认知与技术文化（Tattersall，2009）。这些因素的组合可能共同促成了人类从非洲的扩散，但重建每个因素的具体影响仍是一项挑战。 问题在于，直接实验测试过去的扩散过程是不可能的，且现有经验数据稀缺且不完整。因此，计算机模型被用作“体内”实验平台，模拟不同的扩散情景，评估在现有经验数据下的可能性（Eriksson et al.，2012）。然而，计算机模拟测试替代情景存在局限性。首先，仅能明确测试少部分情景（Reyes-Centeno et al.，2015）。其次，是否最符合经验数据的模型情景是最可能的情景仍不确定。为此，使用近似贝叶斯计算（ABC）方法可以测试成千上万种模型参数组合，评估最可能的参数设置（Marjoram, & Tavare，2006）。这表明，气候、生态系统生产力以及资源利用率在过去的扩散过程中起了重要作用（Eriksson et al.，2012）。 为了解决这些问题，我们提出了一个新的扩散模拟模型，结合了个体和人口级别的参数，以及随时间变化的地理、海平面和气候条件。模拟的目的是探索不同扩散驱动因素的综合效应，而非复制历史扩散事件。我们之前已验证了该模型与理论预期的一致性（Callegari et al.，2013）。 图1展示了我们分析的不同阶段。我们使用基于二十面体细分的球形网格，每个节点与邻居相隔25到33公里。每个网格节点包含生物性（食物可得性）和非生物性（地理特征）属性。食物可得性通过净初级生产力（NPP）衡量，基于模拟的过去区域气温和降水量估算（Timmermann, & Friedrich，2016）。地理特征包括海拔（Amante, & Eakins，2018）、海平面（Waelbroeck et al.，2002）、冰盖（Ganopolski et al.，2010）和现今的河流（Made with Natural Earth，2018）。这些属性决定了每个节点的承载能力。我们使用坡度函数将节点的NPP转化为承载能力K，并根据地理特征修改K。例如，个体不能在冰覆盖或过高的地区生存，而河流和海岸线则增加了节点的承载能力。模拟从85千年前南部非洲的140个个体开始（“预热”阶段），大致与智人从非洲第二阶段的扩散相吻合（Groucutt et al.，2018）。在每个时间步，智能体可以执行不同的行动，如移动、交配和繁殖。智能体的移动方向是一个随机选择，依据邻近节点的海拔和承载能力加权，且彼此独立。 图1 研究设计。通过四种模拟情景（不同节点和种群参数设置）探究晚更新世人类灭绝与扩散的驱动机制。基于网格化模拟计算种群灭绝与迁徙抵达时间，并采用机器学习分析参数空间与模拟结果的关联，用于：结果可预测性评估；敏感性分析；相似抵达时间的参数组合识别。 除了网格节点的属性（即NPP和地理特征），我们还分析了个体移动概率（pm）和人口出生率（b0）对生存和扩散的影响。出生和死亡概率的计算方法见材料与方法。我们模拟了四个复杂度层次的扩散情景（图1）：一个“均匀地球”情景，个体在没有特征的地球上朝任何方向以相等的概率移动；一个“均匀大陆”情景，个体在大陆形状的陆地上以相等的概率移动；一个“恒定NPP”情景，个体生活在时间不变的地形和NPP条件下；一个“动态NPP”情景，个体生活在随时间变化的地形和NPP条件下（Waelbroeck et al.，2002）。由于模拟计算代价高昂，每个85,000年的模拟需要2到5小时，且由于智能体模型的随机性，模拟需要多个重复。为了解决这些问题，我们使用了机器学习（Hastie et al.，2009）。我们构建了数学替代模型，通过对数据子集的测试，确保它们识别出数据的普遍模式，而不受参数值的影响。然后，我们使用替代模型进行敏感性分析，量化每个参数及其交互效应对结果的影响。 通过这种方法，我们回答了两个问题：a.在给定参数下，局部人口是否能生存？b.如果能，人口到达指定区域需要多长时间。我们构建了决策树，揭示了参数间复杂的相互作用，并分析了每个参数如何影响生存概率或到达时间。 二、研究方法 1、计算出生率和死亡率我们为单个个体确定了出生概率pb和死亡概率pd，其中出生概率为pb = b0，死亡概率为pd = d0 + (b0 - d0) * N / K。这里，N是人口规模，K是承载能力，b0和d0是人口规模为零时的出生率和死亡率。 2、确定区域承载能力我们使用坡度函数来确定每个单元的承载能力：若NPP < NPPmax，则K = NPP * Kmax / NPPmax；否则，K = Kmax。 3、Sobol敏感性分析我们使用Sobol敏感性分析方法（Sobol，2001）分析模型输入的变动对机器学习模型输出的敏感性。简而言之，敏感性分析衡量了模型输出的不确定性，这些不确定性可以通过改变模型输入来解释。我们使用Python包SALib实现了Sobol算法（Herman, & Usher，2017）。 分析过程如下：首先，我们指定了每个参数的下限和上限，它们与模拟中使用的最小值和最大值相同。接着，我们使用SALib包中的saltelli.sample函数生成参数值样本，这些样本产生了Sobol序列的参数值组合。所有敏感性分析使用的参数值为N = 3000（生成样本的数量），并设置calc_second_order = True，以计算更高阶的交互效应。根据模拟情景（均匀地球、均匀大陆、恒定NPP和动态NPP），生成了24,000、24,000、30,000和30,000个参数设置，并使用机器学习替代模型预测其行为。需要注意的是，模拟模型的参数比我们在这里测试的要多。我们将这些参数保持不变，因为测试所有参数的影响会过于耗费计算资源。 4、机器学习模型我们使用Python的Scikit-learn库（Pedregosa et al.，2011）构建了所有机器学习模型。对于所有分类和回归估计，我们使用了以下算法：K-近邻（KNN）、逻辑回归、岭回归、随机森林和梯度提升回归树。 KNN是一种基于距离的算法：它通过计算新数据点与训练样本的距离来预测结果。KNN是一个易于理解且快速的算法。如果问题是分类问题，KNN使用多数投票系统决定新数据的类别；如果是回归问题，则取新数据邻近训练样本的均值。 逻辑回归是一种分类算法，使用sigmoid函数为每个新数据点分配概率，并据此分类。 岭回归类似于简单线性回归，但通过最小化模型系数，使其尽可能接近零，这种正则化有助于调整不同特征的重要性，并避免过拟合。 随机森林和梯度提升算法是决策树的集成方法。决策树往往会过拟合数据（即具有高方差），而树的集成方法则有助于利用多棵决策树的能力，减少方差。随机森林是由多棵决策树组成的，每棵树与其他树略有不同。每棵树只对部分训练数据做出较好的预测。随机森林通过向所有树询问新数据来做出预测，分类时使用多数投票，回归时使用均值。梯度提升也使用多棵决策树来做预测，但它通过使每棵新树尽量减少前一棵树的误差来构建树。 三、研究结果 1、人口生存驱动因素我们首先训练了机器学习算法，以根据人口参数预测生存率，这些算法的预测准确性和特异性介于97%和99%之间。然后，使用这些模型对数万个随机参数组合进行了生存预测，并分析了参数变化对模型结果的敏感性。结果（图3）显示，出生率b0是决定人口生存的最重要参数，参数之间的交互作用对生存的影响超过了任何单一参数。参数交互作用发生在多个参数的综合效应非加性时（即一个参数的效应随着其他参数的变化而变化，Berrington de Gonzalez, & Cox，2007）。在四种模拟情景中，参数交互作用分别占生存率变异的18%、19%、10%和5%。即使在最简单的均匀地球模拟中，人口生存仍依赖于复杂的参数交互（图2）。 图2 本决策树基于24,000次模拟结果，展示均质地球情景下种群存活的关键参数阈值。树结构按参数重要性自上而下排列，左右分支分别表示高于/低于阈值时种群的灭绝（左，橙）或存活（右，蓝）趋势。各节点显示：父节点阈值条件下的模拟占比（首行）、灭绝/存活比例（颜色）及子节点阈值标准（末行）。 图3 种群存活参数重要性分析图。展示四种模型情景(颜色图例)中各参数对种群存活的影响：横轴为模型参数，纵轴为一阶Sobol指数(反映参数独立作用)。末位柱状图表示交互作用效应。不同模型参数数量不同(柱状图数量差异)，误差线为95%置信区间。 2、人口扩散驱动因素以下参数变化能增加生存概率：更高的b0、更低的移动概率pm和更高的NPP。Kmax（最大承载能力）的变化对生存概率影响较小。我们选择的当前范围是因为Kmax低于40（每平方公里0.052个个体）会导致许多灭绝，而Kmax超过50（每平方公里0.065个个体）则远离经验估计（Zahid et al.，2016）。b0较高时，生存概率增加，因为它提高了低于承载能力时的出生率。随着pm增加，人口迅速扩展至更大区域，减少了个体的交配机会，从而降低了生存概率。更高的全球NPP通过增加区域的承载能力提高生存概率。 在两个较简单的情景中，参数交互对人口扩散的影响没有生存率那么显著（图4与图3比较）。这表明，在这两个简单情景中，每个参数对到达时间的影响在其他参数值变化时几乎保持不变，而其他参数值不会改变某个参数对到达时间的贡献。然而，在恒定NPP和动态NPP模型中，参数交互可以解释到达时间的大方差（图4）。这表明引入环境异质性后，模型行为变得非线性。b0仍然是决定到达时间的最大因素，甚至比移动概率pm更为重要。唯一例外是从中东到东南亚的扩散，其中pm更为重要。全球NPP和pm在均匀地球和均匀大陆模型中对到达时间的影响程度相似。然而，在恒定NPP和动态NPP模型中，NPPmax和Kmax的影响较小。pm是几乎所有模型中决定到达时间的主要因素（图4）。为了使人口迅速扩散，pm必须足够高，但不能过高，否则会导致灭绝。b0和NPP的重要性在于它们影响人口规模，较大的人口更有可能成功扩散到新区域。 图4 抵达时间参数重要性分析图。展示四种模型情景：a.均质地球(分析非洲之角单区域抵达时间)；b.均质大陆；c.恒定NPP情景；d.动态NPP情景(阿拉伯半岛和东南亚的抵达源自中东)。仅纳入机器学习代理模型可解释90%以上模拟结果的区域。各子图横轴为模型参数，纵轴为一阶Sobol指数(反映参数对采样区域抵达时间的独立影响)，末位柱状图表示交互作用效应。 将人口扩散限制于陆地（均匀大陆情景）并不会显著改变参数对扩散的影响（图3与图4比较）。仅加入地理特征或地理特征和气候波动会增加模拟结果的随机性。具体而言，在恒定NPP和动态NPP情景中，到达时间变得大多不可预测：模型参数本身难以解释模型结果，因为景观和/或气候特征使得到达时间具有随机性。我们还研究了撒哈拉沙漠是否是主要的随机性来源，导致到达时间不可预测。为此，我们构建了机器学习模型，以预测从中东到其他非洲以外地区的到达时间。在恒定NPP情景中，新起始区域仅有助于预测到阿拉伯半岛和东南亚的到达时间，但对于其他地区的预测依然较差。在动态NPP情景中，从中东的到达时间更难预测。 3、等效性为了评估是否存在某些参数值导致相似的到达时间（等效性），我们计算了在动态NPP情景下，导致到达时间相似（误差在1000年以内）的参数空间占比。我们预测了60,000个参数设置的人口生存情况，并确定了存活人群的到达时间。到达时间分布的峰值显示，超过30%的参数设置可以导致相似的到达时间。这些导致相似到达时间的参数设置可能彼此差异很大，差异高达最大可能差异的92%。 我们估计的伊比利亚半岛到达时间，依据60,000次模拟的参数组合，在831到73,831年前之间。到达时间分布显示了三个明显的峰值，其中一个约在55-35千年前，与考古学研究中估计的伊比利亚半岛首次到达时间（约45千年前）大致重合（Cortes-Sánchez et al.，2019；de la Peña，2019）。我们获得了所有导致相似到达时间的不同参数设置，最大误差为1000年，位于50-40千年前区间内。这些参数设置的差异最大可达最大可能差异的66%。例如，人口A的参数为b0=0.35，NPPmax=0.84，Kmax=40.32，pm=0.07，预计会在约44千年前到达伊比利亚半岛，这与人口B的参数（b0=0.32，NPPmax=1.17，Kmax=49.48，pm=0.19）相似。 这些参数差异转化为两个群体之间的生物学差异：人口A每个女性生命周期的平均婴儿数为9.45，而人口B为8.64。A群体的每个个体年均移动2.1公里，而B群体个体的平均移动距离为5.7公里。A群体和B群体的人口密度分别为每千平方公里52.4和64.3个个体。NPPmax和Kmax之间的群体差异表明，人口B在高NPP环境中的利用效率优于人口A。 四、讨论 近年来，越来越多的研究认识到彻底记录、验证和测试模型框架的重要性，我们的研究遵循这一方向（Grimm et al.，2006；Stillman et al.，2015）。与传统的基于“踏脚石”模型的方法（Kimura，1953）不同，我们采用了全球大规模的智能体模型（ABM），并在不同的模型复杂性水平上进行分析。与数学模型（如常用的反应-扩散模型）相比，ABM具有多个优势：能够模拟系统的涌现特性，这些特性来源于个体间的相互作用，并且本质上不可预测（Mitchell，2009）；能够模拟离散现象，融入随机性；以及可以模拟大量异质化个体在异质环境中的相互作用（Macal, & North，2005）。因此，即使是看似简单的系统，也可以表现出复杂的特性，如在均匀地球情景下，参数交互对人口生存和灭绝的复杂影响（图2）。 我们的ABM模型可以用于测试关于人类扩散的各种假设。例如，更高的NPP可能意味着更好的土地利用和狩猎策略，NPP值也可以用于模拟猎物物种的分布。增加的移动概率代表更高的适应能力来应对栖息地迁移（Rolland，2010）。NPP水域（NPPwater）可以代表人类如何利用淡水资源，以及在淡水存在下，区域植被如何改善（Parton et al.，2015）。此外，其他模型参数如分娩间隔、怀孕的最小和最大年龄、人口生育/死亡率，以及个体在高海拔地区的生存能力等，可以用来测试更多关于生育行为和生理学对扩散的影响的假设。 另一创新是使用机器学习分析模拟的参数敏感性，进一步理解ABM的复杂输出。机器学习分析模拟结果的几个优点包括：减少不同参数值组合测试的计算成本，突破了参数空间探索的实际限制；减少了所需的模拟重复次数；揭示了何时输出变异性能通过模型参数解释，何时其他因素（如景观特征和气候变化）影响结果。这使得参数的重要性不仅是系统的全球静态属性，而是空间和时间的函数，从而揭示了历史偶然性在最复杂情景中的作用。 敏感性分析表明，在简单情景中，交互作用对生存变异的贡献较大（均匀地球和均匀大陆情景中分别占18%和19%），但对到达时间的影响较小（4%）。在复杂情景中，根据地区不同，交互作用对到达时间的影响有较大差异。例如，在恒定NPP情景下，从南部非洲到非洲之角及从中东到东南亚的扩散中，交互作用对到达时间影响较小。而在动态NPP情景中，巴塔哥尼亚的路径被气候条件阻塞，只有最后部分的路径对参数效应产生影响。 对模拟输出的可预测性分析表明，气候和地理特征是决定到达时间的最重要因素，超过个体/人口特征（Timmermann, & Friedrich，2016）。在恒定NPP的异质环境中，人口穿越撒哈拉沙漠等通道的过程取决于个体的随机行为，而不太依赖于人口参数。因此，这种模拟中的到达时间不可预测。加入波动气候后，扩散变得更加可预测，因为气候波动完全封锁了重要的通道，所有人群都会被困在这些通道后。当气候打开一个通道时，不同的人群从新的起点开始扩散，使得到达时间更多地成为气候变化的函数，从而变得更加可预测。 从经验数据推断人类扩散情景仍是一个争议问题，研究建议在200千年至50千年之间可能存在单次或多次扩散事件（Reyes-Centeno et al.，2015；Hershkovitz et al.，2018）。基因数据提示人类从非洲的扩散发生在130千年至45千年之间（Scally, & Durbin，2012；Fernandes et al.，2012），一些考古数据表明扩散发生在125千年前（Armitage et al.，2011），而其他数据则认为第一次扩散发生在60千至50千年前（Mellars et al.，2013）。化石记录和基因推断各有局限性，化石记录存在采样偏差，无法确定化石是否为谱系的首次证据，基因推断模型的假设可能不现实，放宽假设会导致不同结果（Chikhi et al.，2010）。化石记录表明，现代人类在130千年至80千年前已存在于黎凡特和阿拉伯半岛（Grun, & Stringer，1991；Grun et al.，2005），并可能在55千年前重新进入黎凡特（Hershkovitz et al.，2015）。西欧和西伯利亚的首次定居大约发生在45千年前（Fu et al.，2014；Benazzi et al.，2011）。 没有单一的“最佳”参数设置能够完全符合这些经验数据。原因有三：首先，参数间的相互依赖性强（图3和图4）；其次，许多参数设置，即使差异很大（最大差异可达92%），也能导致相似的结果（等效性问题）；第三，类似的参数设置可能因巨大的随机性（历史偶然性）导致不同结果。这表明，任何旨在推测人类到达时间的合理参数值的模型，都必须测试一系列的参数值。报告一个单一的最佳参数组合，可能会误导结论。我们再次强调，我们的模型并非为了拟合经验数据，而是展示了不同参数空间的可能性。进一步的经验数据可以缩小人类扩散和人口动态的潜在参数空间，从而减少计算模型的不确定性。 我们验证并分析了一个复杂的人类扩散计算模型。我们展示了在简单景观中，到达时间是模型参数的函数，且模型行为中的随机性较小。然而，加入地理特征和气候波动后，随机性显著增加，模型行为变得难以预测。预测能力的局限性表明，历史偶然性在人类从非洲的扩散中发挥了重要作用。即使模型表现出可预测性，不同参数设置也可能导致相似的到达时间，这表明模型验证在得出关键参数结论前至关重要。此研究框架可用于评估人类扩散的替代理论，或探索为什么现代人类能胜过尼安德特人，模型两种物种在不同参数组合下的竞争。
