{'fraction_non_alphanumeric': 0.03098752989071362, 'fraction_numerical': 0.013695381514199323, 'mean_word_length': 4.891489113983002, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 3, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Global hunter-gatherer population densities constrained by influence of seasonality on diet composition

The dependence of hunter-gatherers on local net primary production (NPP) to provide food played a major role in shaping long-term human population dynamics. Observations of contemporary hunter-gatherers have shown an overall correlation between population density and annual NPP, but with a thousand-fold variation in population density per unit NPP that remains unexplained. Here we build a process-based hunter-gatherer population model embedded within a global terrestrial biosphere model, which explicitly addresses the extraction of NPP through dynamically-allocated hunting and gathering activities. The emergent results reveal a strong, previously unrecognized effect of seasonality on population density via diet composition, whereby hunter-gatherers consume high fractions of meat in regions where growing seasons are short, leading to greatly reduced population density due to trophic inefficiency. This seasonal carnivory bottleneck largely explains the wide variation in population density per unit NPP, and questions the prevailing usage of annual NPP as the proxy of carrying capacity for ancient humans. Our process-based approach has the potential to greatly refine our understanding of dynamical responses of ancient human populations to past environmental changes.


Introduction
Hunter-gatherer populations, that subsist on hunting, gathering and fishing, rely intimately on the biotic fabrics of their local environments [1][2][3] . Although only a handful of hunter gatherer populations remain, over 300 societies of contemporary hunter-gatherers have been documented sporadically over the past two centuries 1 . Due to the assumed similarities between their lifestyles and those of our foraging ancestors, contemporary hunter-gatherers have provided many insights on our species' past [4][5][6][7][8] , informing reconstructions of population changes on evolutionary timescales 5,6 , the timing and rates of global human dispersals 4,7 , and the drivers behind the origin of agriculture 8 . Yet, these contemporary societies do not provide direct analogues of ancient foragers, as they had more complex technologies 9 and had experienced interactions with and pressures from neighbouring agricultural and industrial societies 10,11 , such as acquiring supplemental agricultural food 2,12 or encountering novel pathogens near the time of documentation, which may have modulated the population density-ecosystem relationships. Thus, many of the observational data may provide distorted views on the pre-existing state, casting doubt on statistical models and demographic parameters that are directly fitted to the contemporary data [5][6][7][8] . Taking full advantage of the insights from contemporary hunter-gatherers requires a mechanistic, process-oriented understanding of how environmental factors influence the distribution and abundance of hunter-gatherer populations.
Because hunter-gatherers acquire food directly from their surrounding environment, it has long been thought that their population density should be closely linked to the productivity of their local ecosystems 1,2 , and numerous studies have used net primary production (NPP) as the main predictor for ancient hunter-gatherer density changes 4,7,13 . However, although contemporary hunter-gatherer population density is positively correlated with NPP, NPP alone explains less than 30% of the variability of population density (Fig.  1a). The fraction of NPP consumed by hunter-gatherers (Φ NPP , see Methods), which also indicates population density per unit NPP, varies by three orders of magnitude across the recorded groups ( Fig. 1b and Supplementary Figure 1e). A portion of this scatter could reflect the confounding historical influences, but it has also been assumed that the edible proportion of primary production differs among biomes 2 . However, our analysis does not show systematic differences in Φ NPP between biome categories, except for a lower Φ NPP in boreal forest and a marginally significantly higher Φ NPP in Mediterranean forest ( Fig. 1b and Supplementary Figure 2b). Nor does Φ NPP show a consistent change across the NPP gradient (Supplementary Figure 2c). The wide spread in Φ NPP thus remains unexplained. The weak statistical dependence of population density on NPP has prompted a recourse to other explanatory variables, including biodiversity 3 , pathogen stress 3 , precipitation seasonality 8 , climate variability 14 , and social complexity-related variables such as food storage-dependence 14 . It is however difficult to tease apart the causality among these inter correlated variables using statistical methods alone.
Here we build a mechanistic understanding of hunter-gatherer populations with a global, activity-based process model for which population density is an emergent feature. Our model operates within the framework of an Earth system model (ESM) (Fig. 2). During the past decades, the land components of ESMs have integrated vegetation modules to resolve key processes of biogeochemical cycling of carbon, energy and water 15 , and have recently begun to incorporate large mammalian herbivores based on metabolic and demographic equations 16,17 . Our model directly couples human population dynamics to such a land model 17 , simulating human time allocation to hunting versus gathering in response to the interaction between humans and food resources (see Methods). The model uses explicit formulations of daily carbon/energy flows among vegetation, herbivores and humans, the outcome of which determines human reproduction and mortality rates. These fluxes depend on a time allocation scheme with two simple assumptions: first, total foraging time increases or decreases depending on the level of fullness of the previous day, and second, gathering vs. hunting time depends on the relative abundance of plant vs. animal food and on an underlying a priori preference for meat. The model is resolved on a daily time-step, thus capturing the seasonal cycle in high detail.

Modelled global distribution of hunter-gatherers
Simulated steady-state global population density and foraging times under present-day climate are shown in Fig. 3. The model simulates high densities in most regions with high NPP and closely matches the observed relationship of density versus NPP derived from contemporary hunter-gatherer data (Fig. 3b). Simulated foraging times are generally 3.5-5.5 hours per person per day throughout the world, close to the average of available observations from a few hunter-gatherer groups 2,18,19 (ca. 4 hr. day -1 ind. -1 ). In the model, slightly more time is spent hunting than gathering in the tropics, whereas most of the foraging time is devoted to hunting at high latitudes (Fig. 3cde).
The simulated population density agrees reasonably well with the ethnographic data (Pearson correlation coefficient r=0.58, Supplementary Figure 3ab), and there is also a strong correlation between simulated and ethnographically recorded diet compositions, though less of the observed variance is explained (r=0.40 for the meat fraction of the diet, Supplementary Figure 3cd). Model-data discrepancies may be partly due to the lack of fishing in the model, which is a major source of food in many coastal societies (Supplementary Figure 1d). In addition, the model treats all plant or animal food equally, thus ignoring the wide variety such as the different energy content, protein richness, and relative abundance of each food item within both broad categories 2,[20][21][22] . Any bias in the vegetation and herbivore dynamics model 17 could also propagate to the human model, contributing to the model-data mismatch (note that the S2 experiment has corrected part of the biases in the vegetation inputs to the FORGE model, see Methods). Furthermore, discrepancies could arise from climate shifts and non-steady state factors given that the model simulates equilibrium densities forced by present-day climate (while the ethnographic records were collected during the past two hundred years 1 ) whereas some of the populations were under positive growth rate when being studied 23 . There are also significant uncertainties in the diet composition data, as the estimation methods for the amounts of gathered, hunted, and fished food are often inconsistent, sometimes measured in different units (weight or calorie) and occasionally based on the ethnographers' impressions 2 , which may further add to the model-data scatter.
The modelled global total hunter-gatherer population is 17 million, which is at the high end of the estimates of prehistoric (pre-agriculture) population derived by extrapolating national historical records, which range from 1-20 million 24 . This can be partly explained by the fact that our model was calibrated with data from contemporary hunter-gatherers, who presumably had more advanced technologies and access to non-local foods than early-Holocene foragers 25 (see Supplementary Discussion 2 for the test runs regarding the technology-dependent foraging efficiency parameter). Geographical differences also contribute to our estimate being higher than previous upscalings of contemporary hunter gatherer populations, which include regression models (10~12 million, refs. 3,26 ) and a proportional projection method (7 million, ref. 1 ). For the former, the difference is primarily attributable to a higher density in Africa and southern Asia in our model (Supplementary Figure 4), where the populations were suppressed in ref. 3 due to a high pathogen stress as a predictor in their regression (see Supplementary Discussion 1 for more discussion about pathogen stress). For the latter, we find that the global estimate with the same projection method but a different upscaling biome map ranges from 10 to 17 million, contingent on how an average population density is calculated for each biome (see Supplementary Discussion 3). Therefore, we consider our result of 17 million, simulated by a process-based model, to be a reasonable, independent estimate of potential global hunter-gatherer population size under modern climate with advanced foraging technology, assuming foraging as the only subsistence type.

Linkage among seasonality, diet composition and population density
As in the observations, simulated Φ NPP (fraction of NPP consumed by hunter-gatherers) displays a huge spatial variation with a thousand-fold range (Fig. 4b). The spatial distribution of Φ NPP follows an inverse relationship with the meat fraction of the diet (compare Figs. 4a and 4b). The negative correlation is clear in the scatter plot, both for grid cells where the observational data are located (Fig. 4d) and for all populated grid cells in the model (Supplementary Figure 5), showing a drop of Φ NPP by almost an order of magnitude as meat fraction increases by 60%. Remarkably, this emergent property of the model is corroborated by ethnographic data, which falls over a comparable range (compare Figs. 4c and 4d). Indeed, in both the observations and the model, a higher meat fraction of the diet is associated with a lower population density under the same level of NPP (Supplementary Figure 6). It is possible that the use of aquatic resources, which are not included explicitly in the model, could increase the density of some populations relative to terrestrial NPP (and thus Φ NPP ) and thus confound the model-observation comparison. However, we checked that the observed correlation between Φ NPP and meat fraction holds when the groups with fish-rich diets are excluded (Supplementary Figure 7). Therefore, although aquatic resources are important for hunter-gatherers in some locations (Supplementary Figure  1d), the association between Φ NPP and meat fraction appears to be a robust feature of the observations. Notably, even the groups with intermediate (30%≤fishing<50%) or high (≥50%) fishing percentages display significant negative correlations between Φ NPP and meat fraction as well (Supplementary Figure 7), lending confidence to this overall trend.
The meat fraction of the modelled hunter-gatherer diet is mainly driven by a scarcity of vegetation edible by humans during winter or dry seasons in the model. In regions with relatively small seasonal variations in NPP, foragers can subsist mainly on gathering throughout the year due to a substantially higher abundance of plant food than animal food at equilibrium (Fig. 5a). By contrast, in regions with long non-growing seasons, edible plant biomass becomes depleted and foragers depend on hunting to subsist through the season of scarcity. As a result, the production rate of animals imposes a strong constraint on human population density in regions with short growing seasons, regardless of how abundant plant food may be during the growing season (Fig. 5b). Since relatively little energy flows from primary to secondary production due to trophic inefficiency (Fig. 5), the Φ NPP is expected to be lower for humans subsisting more on animals. Therefore, the apparent negative correlation between Φ NPP and meat fraction ( Fig. 4d and Supplementary Figure 5) is caused both directly by the effect of shifts in humans' trophic position, and indirectly by the covariation between meat fraction and seasonality, the latter being able to influence Φ NPP through the limiting seasonal minima.
With the process-based model, we quantify the contribution of the two effects, by conducting two experiments in which humans only subsist on either plants or animals (via allocating all foraging time in either gathering or hunting, hereafter "onlyGather" and "onlyHunt"). The difference in the resulted Φ NPP between the two runs could be regarded as the direct impact of the trophic position. In the results of onlyGather, humans only persist in regions with long growing seasons (> ca.180 days), but the total population size in these populated areas is 9.5 times higher than that in the onlyHunt experiment for the same area (Supplementary Figure 8). Correspondingly, Φ NPP is mostly higher in onlyGather than in onlyHunt (albeit the opposite occurs in a few grid cells with relatively shorter growing seasons where a full reliance on hunting slightly increases population density), and log 10 (Φ NPP ) increases by 0.85 (or, Φ NPP increases by seven times) on average in onlyGather compared to onlyHunt (Supplementary Figure 8e). Therefore, the direct effect of diet composition via trophic energetics contributes about half to the apparent slope of Φ NPP versus meat fraction.
The modelled contrast of seasonal cycles in diet composition between long versus short growing season regions (Fig. 5) is qualitatively supported by observations in the Ache (in tropical forest) 20 , and in the Hiwi (in savanna) 27 and Bushmen groups (in Kalahari desert) 28 . The former, living in places with relatively high NPP all year round, exhibits little seasonal variance in both meat and plant intake, whereas the latter two, living in highly seasonal environments, exhibit a highly varied share of plant food across seasons. An in-depth quantitative evaluation is however difficult due to the limited number of hunter-gatherer groups with detailed, seasonally based diet records, and because the model cannot capture the highly diverse phenology of various plant and animal species at the local ecosystems.
Given the strong role of meat fraction in controlling population density predicted by the process-based model, we carried out a statistical analysis of contemporary hunter gatherer data in order to test if meat fraction emerges as a significant predictor (see Supplementary Discussion 1). Our analysis shows that, indeed, a multivariable stepwise linear regression with nine empirical predictors (NPP, tree and grass coverage, biodiversity, growing season length, percentages of food derived from hunting, gathering and fishing, and absolute latitude) identifies the hunting percentage as the strongest explanatory variable for population density, followed by NPP and grass cover, while the other variables are non significant. Hunting percentage itself can be explained by growing season length, followed by grass coverage (see Supplementary Discussion 1). Whereas it has been previously noted that the fraction of animal food tends to increase at higher latitudes 25,29 , our result suggests that it is mainly the short growing season that explains this gradient, as the effect of latitude itself is non-significant when growing season length is considered. A structural equation model, which considers both direct and indirect effects of the predictors, also confirms the dominant effect of hunting percentage on population density in the observational dataset, and supports the link that seasonality influences population density via the hunting percentage (Supplementary Figure 9 and Supplementary Table 1). Relationships between Φ NPP and other variables, apart from hunting percentage (Fig. 4c), were also examined (Supplementary Figure 10), showing a significant positive correlation between Φ NPP and growing season length (r=0.37, p<0.001), which however becomes non-significant after controlling for hunting percentage in the partial correlation analysis (Supplementary Table  2).

Discussion
This study built a process-based hunter-gatherer dynamics model, coupled to a global terrestrial biosphere model, to explore the influence of environmental factors on hunter gatherer population density. The model explicitly simulates micro-scale processes including daily carbon/energy flows among plants, animals and humans, combined with a dynamic time allocation for hunting and gathering activities with simple assumptions. The emergent macroscopic relationships from the model are well supported by ethnographic observations of contemporary hunter-gatherers, revealing a causal mechanism whereby short growing seasons drive high fractions of meat in the diet, leading to greatly reduced population density per unit NPP largely as a result of trophic inefficiency.
The empirical negative association between hunting fraction and population density has been reported or implied in previous studies 14,30 , but it was assumed that population density was the causal factor, rather than vice versa: a higher population density was suggested to lead to less frequent residential moves, more food storage, and thus greater dependence on the stationary and higher productive food of plants and coastal resources 14 . Our process-based coupled human-ecosystem model suggests a more parsimonious explanation, given that it reproduces the observed relationships as a straightforward consequence of bioenergetics. Reinforcing mechanisms between demography and diet via cultural factors may still be relevant, as suggested by ref. 14 and by anthropological studies on diet from both ecological/ biological and sociocultural perspectives 31 . In particular, food storage and sharing have received much attention, due to their close association with societal characteristics and ability to mitigate food scarcity at daily to seasonal timescales [32][33][34] . It may be interesting to include food storage and sharing in our model framework in the future to quantitatively explore the interplay between environment, demography, diet, and sociocultural aspects.
The integral links between seasonality, diet composition and population density have strong implications for studies that explore ancient human populations. Although seasonal cycles in the availability of different types of food have left traces in fossil teeth of extinct hominins 35 and been detected in the gut microbiome composition of the Hadza foragers in East Africa 36 , seasonality has not been included in prior modelling studies 4,6,7,13 that used only annual NPP and/or mean paleoclimate variables, which would have biased expectations for the spatio-temporal pattern of early human distribution and migration. Furthermore, the lengthening of growing seasons following the last deglaciation, combined with the end-Pleistocene megafaunal extinctions that occurred in different regions at different timing and rates 37 , might have driven long-term changes in diet composition towards less meat dependence. Our results suggest that a shift from meat-dominated to plant-dominated diet could have boosted population much more than implied by NPP changes alone. Given that demography might be at the core of cultural evolution 38 , which further shaped human societies 8,39 , it is crucial to study the trajectory of past population growth and its drivers, under the changing paleoclimate, vegetation and animal distributions.
The model developed here represents a new type of computational modeling for ancient human-ecology studies, breaking from a tradition 40 that has mainly included niche (species distribution) models based on statistical methods and without human feedbacks on the environment, and agent-based models which resolve individual behaviors and interactions but require very detailed local information that is challenging to assess and have only been applied at a local/community scale 41,42 . By incorporating hunter-gatherers in a realistic, interactive and dynamical global environment, the ESM framework helps to clarify how short growing seasons force humans to rely on meat-rich diets, reducing Φ NPP and thus the population density per unit NPP. Additional factors beyond those resolved in our current model, such as technological transitions 25,43 and long-range migrations, are sure to have played important roles in determining changes in hunter-gatherer population abundance and diet through time. With further development, the approach has the potential to provide more highly resolved pictures of the complex, multifaceted interactions between our ancestors and the Earth system.

Contemporary hunter-gatherer data
The hunter-gatherer population density and coordinates were acquired from ref. 3 , which combined the ethnographic datasets from refs. 1,2,44 (357 data points in total). For the diet composition (percentages of hunting, gathering, and fishing), we combined the two datasets of subsistence from Kelly 2 and Binford 1,45 (340 data points in total), which indicate proportions of food derived from terrestrial animals, terrestrial plants, and aquatic resources 1 , in terms of weight or calorie 2 . Average of the two datasets were used if a population is present in both sources.
To explore the potential impact of environmental variables on population density and diets, we extracted the values of the following variables from their global maps: annual NPP from the satellite-derived MOD17A3 product 46 (averaged over 2001-2010); fractional coverage of trees and grasses from a satellite-derived vegetation cover product 47 ; biodiversity index from ref. 3 , which was the average of scaled richness of animals (mammals + birds) and vascular plants; growing season length (GSL) calculated as the days in a year when daily gross primary production (GPP) exceeds 20% of the maximum daily GPP for each grid cell, wherein the daily GPP was from an upscaled global map based on FLUXNET tower sites 48 (different thresholds were tested but made little difference, e.g. GSL defined by a threshold of 10% and 20% are highly correlated, r=0.93). Due to slight differences of the land-sea mask of these global land-only products, a few coastal populations were located in grid cells with a non-valid value; in such case, the value of the nearest land pixel was used. These potential environmental predictors are used in statistical analyses of the ethnographic data as detailed in Supplementary Discussion 1.
Unlike Tallavaara et al. 3 , we did not use the climate-based Miami model to calculate NPP but instead used MODIS NPP, because the Miami model was shown to overestimate NPP in the tropics since it only considers annual mean temperature and precipitation and misses important factors like nutrient limitation in tropical forests 49 . It was argued in ref. 3 that the satellite-based NPP was unfavourable because of the recent human appropriation of NPP. However, we consider it a minor problem here because the majority of the sites have very small coverages of cropland, according to a satellite product 50 (242 out of the 357 sites have <5% cropland cover, 305 sites <20%).
We calculated the fraction of NPP consumed by hunter-gatherers (Φ NPP ), which represents the proportion of total energy available to heterotrophs that flows into hunter-gatherers, as follows: per capita daily energy consumption, averaged over eight different groups 2 , is 2480 kcal ind. -1 day -1 , which equals to 174 kgC ind. -1 yr -1 , using the conversion coefficients of 9.8 MJ kgDM -1 (see the text below Eq. 4) and 0.45 gC gDM -1 . This same value is multiplied by population density and divided by NPP to derive Φ NPP for each site.
A potential systematic bias in the ethnographic diet composition data is worth attention. In this study we regard the recorded hunting fractions as representing meat fractions of the diet. However, the hunting fractions are likely to underestimate the true consumption of animals due to the fact that early observers had devalued women's contribution in meat provision, including small animals, eggs and insects. Such contributions had been often categorized as "gathering", or simply ignored, as the focus of "hunting" had been on relatively large, mobile prey. This bias could be large for groups in arid regions, e.g., the Western Desert in Australia 51,52 , for which the recorded hunting percentage is about 30% while the true meat fractions are up to 80% when women's contribution to hunting is included 51 . The bias appears to be lower for groups in the tropics like Amazon Basin and Congo Basin 53 , where recent field observations show less than 10% of total hunted food from women 53 . Unfortunately, revisiting and verifying the secondary data on diet composition is difficult due to the disappearance of the foraging lifestyles for most of the recorded hunter-gatherer groups. We therefore conducted a sensitivity test to see how the observed negative relationship between Φ NPP and hunting fraction (Fig. 4c) might be influenced by expected uncertainties and biases in the ethnographic dietary records. As detailed in Supplementary Discussion 1, we found that the relationship between Φ NPP and hunting fraction is robust to the expected biases, with only a minor impact on the slope of the correlation over the most probable range of bias.

Description of the process-based hunter-gatherer dynamics model
Model overview-We designed the hunter-gatherer dynamics model, FORGE (FORager dynamics in Global Ecosystem model), in the framework of Earth system models. It is coupled to the ORCHIDEE (Organising Carbon and Hydrology In Dynamic Ecosystems) terrestrial biosphere model 54 , which is the land component of the IPSL (Institut Pierre Simon Laplace) ESM. ORCHIDEE simulates the energy and water balance, vegetation dynamics, and carbon cycle of land ecosystems. Inputs for ORCHIDEE include meteorological variables (air temperature, precipitation, incoming short and long wave radiation, wind speed, air humidity, and air pressure), atmospheric CO 2 concentration, land cover, topography and soil texture maps. Spatial resolution and domain of each simulation are customized depending on the inputted climate forcing. In the model, each grid cell is occupied by a suite of plant functional types (PFTs), with their fractional covers adding up to one. ORCHIDEE has recently incorporated a module of large mammalian herbivores that simulates animal population density driven by vegetation and climate 17 . The main revision of the herbivore module compared to ref. 17 is the inclusion of browsers (herbivores that eat woody plants) whereas the previous version only calculated grazers. The same parameterizations are used for browsers as grazers except that browsers feed on the leaf and fruit compartments of tree PFTs, assuming 10% of these biomass is reachable to browsers. The summation of grazer and browser density in each grid cell provides the potential animal food density for humans.
The new FORGE model (Fig. 2) couples hunter-gatherers with vegetation and herbivores in each grid cell via daily foraging activities and resultant carbon/energy flows, which then updates an energetic reserve (body fat) that determines birth rate and mortality rate of human populations. The human foraging activities are formulated based on a time allocation scheme. The feedback of human hunting on herbivore dynamics is taken into account, while the impact of humans on vegetation growth is neglected in the current model, but the model infrastructure leaves room for future developments regarding human-induced environmental changes (e.g. changes in land cover, fire regimes, etc.) and subsequent impacts on the ecosystems. Major simplifications of FORGE include that i) only an average human (body weight = 50 kg ind -1 ) is represented; ii) fishing is not implemented in the current land-only model, leading to a potential underestimation in population density of the groups for which marine resources provide much food, mostly in coastal regions (Supplementary Figure 1); iii) human migration 55 across grid cells is not accounted for, although the fact that the ethnographically recorded average migration distance per move and the total distance moved between camps in a year are ca. 27 km and 280 km, respectively 1 , means that most hunter gatherers would not move over more than the dimensions of one grid cell (on the order of 10 4 km 2 each) in an annual cycle; and (iv) the lack of representation of small animals in the model, including birds, reptiles, insects, and small mammals (<5~10 kg ind. -1 ), could lead to underestimations in available animal food and meat fraction of the diet, which coincides with the potential systematic underestimations in the ethnographic hunting fraction data mentioned above. Detailed formulations of FORGE are described below.
Daily intake-The model considers two types of subsistence activity, hunting and gathering. The contribution of each activity to the diet depends on the abundance of animal biomass and edible vegetation biomass, the time spent in hunting and gathering, as well as the technology-dependent efficiency of both activities: I V = B V t g e g a I A = B A t ℎ e ℎ a (1) where I V and I A are daily dry matter intake of vegetation and meat (kgDM day -1 ind. -1 ); B V and B A are biomass density of edible vegetation and animals (kgDM m -2 ); t g and t h are the time allocated to gathering and hunting (hr day -1 ); a is the searching area per person hour, fixed at 4000 m 2 hr -1 ind. -1 ; e g and e h are the efficiencies of gathering and hunting, representing the fraction of the corresponding biomass that is acquired when hunter-gatherers pass the searched area.
In a linear form this equation would imply that the food resources and foragers are both uniformly distributed within each grid cell and foragers encounter them randomly. In reality, excursions are undertaken from localized camps, different food items are distributed in patches, and foragers know how to target the higher energy-return items first 2 . To partly account for these aspects of sub-grid heterogeneity, we assume e g and e h to decrease as B v and B a decrease relative to human density, capturing the increasing difficulty of foraging as the best resources are sequentially depleted. We model this with an asymptotic function: where B H is human biomass density (kg m -2 ), equal to P H (population density which is prognostically simulated in the model, see Eq. 11) multiplied by mean body weight (W H , 50 kg ind -1 ); e max is the maximum fraction of the edible plant/animal biomass that can be acquired when hunter-gatherers pass the searched area, fixed at 1 and tested in a sensitivity analysis (see below); c determines the patch depletion rate; the value of c is fixed at 100, and is tested in the sensitivity analysis. Notably, the inclusion of Eq. 2 produces a dynamic similar to that found in ecological studies of per capita prey consumption by wild predators, which often decreases with predator density for a given prey density 56 .
Energy expenditure and energetic reserve-The per capita daily energy expenditure E (MJ day -1 ind. -1 ) is calculated as: where E c is the energy expenditure excluding those spent during foraging, fixed as 8.37 MJ day -1 ind. -1 (i.e. 2000 kcal day -1 ind. -1 ); E g and E h are the energy spent per hour gathering or hunting, fixed as 1.25 MJ hr -1 ind. -1 (according to ref. 57 ) for both parameters in the model. Based on daily food intake and energy expenditure, per capita energetic reserve (mainly in fat cells) is updated daily as below. This is the key variable that indicates health condition and is used to calculate birth rate and mortality rate (see below). (ref. 58 ). The gross energy value of plants and animals is around 3800-4800 for different plants (ca. 5000 for seeds) and 4400-5600 for different animal taxa (ref. 59 , in the unit kcal kgDM -1 ). We thus assume the same value for both food types for simplicity. Further considering the ca. 30% energy loss in excreta, and that the metabolisable energy content is less than the gross energy value of food, we finally set m V and m A to 9.8 MJ kgDM -1 . It should be noted that the energy density of a mammal's carcass depends on its fat content 22 , thus a fixed m A value in the current model is a simplification and does not account for variations across species and over time.
There is an upper limit of F (F max ), assumed to be 30% of body weight, i.e. 15 kg ind. -1 in the model. Accordingly, whenever F is about to exceed F max , daily intake will be reduced so that F stays at F max .  (6) where th frac is the fraction of t f allocated to hunting. The additional meat-craving term, C meat , is introduced to represent that hunter-gatherers generally have a preference for meat, which may reflect nutritional needs (e.g. fat and essential amino acids) 62 as well as cultural importance 2 . The parameter q is fixed at 2.5 but also tested (see Supplementary Table 3). In cases when f meat falls below 10% (the minimum observed hunting+fishing fraction), th frac is set to 0.95. t h and t g are then calculated as t f × th frac and t f × (1 − th frac ) respectively. This time allocation scheme is highly simplified. It does not consider, for example, seasonal variations of other necessary activities such as providing shelter, extreme weather that makes foraging activities risky or unworthy (e.g. the Hiwi in tropical savannah don't forage during the middle of the day to avoid heat 27 ), and additional social/cultural factors. Further improvements and tests of the scheme are possible in the future, given more observational data on human time allocation dynamics. Yet, it provides an unprecedented approach to include the fundamental limit of time on the ability to capture existing biomass.
Updating food resources-Not all vegetation and herbivore biomass are edible or accessible to humans. We denote the edible and/or accessible fraction of vegetation and herbivores as Φ V and Φ A , respectively, which are assumed to be fixed parameters in the current model over the globe (Φ V =0.015, Φ V =0.1, see sensitivity tests below). This is a simplification; for example, many fruits in tropical forests grow high in the canopy, making them harder to reach than the case in short-vegetation biomes 2 . Unfortunately, robust estimates of the edible/accessible fraction for major biomes worldwide are unavailable, and therefore we do not vary Φ V and Φ A between grid cells in the current model, in order to avoid arbitrarily introducing spatial relationships of human density versus NPP. Apart from fruits, hunter-gatherers can also eat roots and tuber, which are from the belowground primary production of some herbaceous plants. We thus also include part of the belowground grass NPP as potential food sources.

Europe PMC Funders Author Manuscripts
Where λ (day -1 ) is the daily turnover rate of the edible vegetation biomass; f fruit is the fraction of NPP (including both tree and grass) allocated to fruits, set to 6.5% in the current model according to a synthesis 63 of fruit fall observations across tropical and temperate forests; f below is the fraction of grass NPP allocated belowground, calculated as a function of mean annual precipitation (MAP, mm yr -1 ) following the empirical equation from ref. 64 , which captures the higher belowground allocations in dryer ecosystems. The lower limit of f below is set to 20%. The inputted NPP and λ are either from ORCHIDEE outputs or from observation-derived values (see "Model setup" section for the different runs).
The herbivore dynamics module is coupled once per year with the human dynamics module. Detailed descriptions of the herbivore module can be found in ref. 17 . In summary, herbivore population density (P A,tot , ind. km -2 ) is updated on a yearly time-step using a logistic equation (Eq. 8 in ref. 17 ). The herbivore birth and mortality rates depend on an animal fat reserve pool (kg ind. -1 ) that is updated daily according to their energy intake through grazing or browsing.
Within the human module, B A is updated on a daily time-step: At the end of each year, the accumulated number of hunted animals ∑I A P H W A × 0.28 is subtracted from the herbivore population density P A,tot . The decreased P A,tot is passed to the herbivore module to calculate a new P A,tot for the next year. B A is then updated and used in the human module for the new year: B A = Φ A P A, tot W A × 0.28 (9) where W A is the mean body weight of herbivores, set to 180 kg ind. -1 in the current model; 0.28 is the conversion between live weight and dry mass after excluding water and bones.
Note that the herbivore population dynamics is also affected by seasonality, with relatively higher animal starvation mortalities in regions with longer winters/dry seasons 17 . However, temporally, simulated meat density is relatively stable throughout the year in the model (see Fig. 5), because the animal population density is updated only once per year. The annual update is justified by the fact that the animals represented have lifespans of ~25 years 17 , so that seasonal fluctuations in population numbers are strongly damped. Unlike humans, herbivores can feed on plant litter (dead grasses and fallen leaves/fruits in the model) to live through non-growing seasons, providing vital food for hunter-gatherers to survive the long winters/dry seasons in the model.
It should be emphasized that the modelled feedbacks of human foraging are very different for vegetation and animal density. Gathering is assumed to have negligible impacts on primary productivity, since the B v that is harvested represents the edible part of plants, which generally do not contribute to plant growth, and for which harvesting can even aid in reproduction (e.g. through the spreading of seeds 65 ). In contrast, hunting can directly reduce the annual production of herbivores by reducing P A,tot . However this effect would be small at low rates of hunting, for which competition among herbivores places the limit on P A,tot .
Human population dynamics-The population birth rate (r birth , yr -1 ) is assumed to depend on the average body condition, represented by the daily-varying energetic reserve F. Thus, r birth is calculated every day and averaged over the year to be used in Eq. (11), where r a = 0.1, r b = 15, r c = 0.5. This sigmoidal function is similar to the equation used to calculate herbivore birth rate 17 , with the parameters modified so as to give a maximum birth rate of 10%, close to the recorded highest crude birth rate during the past two centuries (ca. 6%, ref 66 ) (see sensitivity tests of the three parameters below).
For mortality, we consider two processes. First, a background mortality rate (M b ) which is the inverse of lifespan, fixed as where P H (ind. km -2 ) is the human population density for each grid cell. The P H is initialized as P 0 =10 -3 ind. km -2 , a lower value than the minimum density ever recorded (0.002 ind. km -2 , ref. 1,2 ). Whenever P H falls below P 0 , P H is reset to P 0 .
Once each year, the reproduction energy cost (313 kcal day -1 , or 478 MJ for a year, according to ref. 67 ) is subtracted from the energetic reserve F in order to account for the energy input to birth: Model setup-To derive the input variables for FORGE, we firstly ran ORCHIDEE for the globe at 2° spatial resolution and for 200 years as a spin-up to reach the equilibrium of vegetation production and biomass under present-day climate conditions. The climate forcing was the CRU-NCEP reanalysis dataset 68  The last 10 years' outputs from ORCHIDEE are used to provide the inputs to FORGE. These include the inputs required by the herbivore module: fractional cover of PFTs (aggregated into three types: grass, tree, bare ground), carbon influx rates (i.e. NPP allocated to the edible plant tissues, and influx to the edible litter pool) and decay rates of the edible pools, and annual mean temperature; as well as inputs required by the human module: NPP, decay rate of the fruit compartment, and annual precipitation (used in Eq. 7). The input variables are at a daily time step. FORGE was then run for 300 years to reach the equilibrium of population density of both herbivores and hunter-gatherers. The last 50 years were averaged and presented as the results.
Supplementary Figure 11 shows the ORCHIDEE-simulated NPP, in comparison with MODIS NPP. ORCHIDEE can reproduce the overall pattern of the satellite-derived NPP, with a reasonable agreement with MODIS NPP at the sites (Pearson correlation coefficient r=0.67). However, ORCHIDEE underestimated NPP in tropical forests and arid ecosystems including some classical hunter-gatherer areas like the Great Basin in North America, Kalahari in southern Africa, and interior Australia (Supplementary Figure 11). The underestimation of annual NPP in the dry regions is partly due to an underestimation of growing season length by ORCHIDEE (Supplementary Figure 12).
In order to test the impact of the bias of ORCHIDEE-simulated NPP, and to minimize this bias that would propagate to the FORGE-simulated hunter-gatherer densities, we conducted three sets of simulations, with different input files while identical FORGE model: S0: inputs to FORGE are directly from ORCHIDEE outputs; S1: compared to S0, the daily NPP is multiplied by a scaling factor so that the annual NPP (average for 2001-2010) equals to MODIS NPP for each grid cell; S2: compared to S0, the daily NPP is replaced by the observation-based values, that is, annual MODIS NPP interpolated into daily values according to an observation-derived daily GPP product 48 . Note that the herbivore module in FORGE separates grazers feeding on grass NPP and browsers feeding on tree NPP, and here we assume the same NPP values (per unit PFT area) for tree and grass PFTs in the same grid cell. Besides, the ORCHIDEE-simulated daily decay rates of edible plant tissues (λ) are replaced by a constant of 0.01 day -1 during growing season (when NPP daily >NPP daily,max ) and 0.05 day -1 during non-growing season (when NPP daily <=NPP daily,max ), in order to avoid an inconsistent timing of senescence (a period of high decay rates of leaves and fruits) with the observed seasonality in NPP. The λ values are chosen so that the integrated annual decay rates are generally comparable to those simulated by ORCHIDEE (used in S0 and S1).
Modelled population densities and their relationships between NPP and meat fraction of the diet from the S0 and S1 runs are shown in Supplementary Figures 13 and 14. Among S0, S1, and S2, the emergent relationships between population density and NPP, and Φ NPP and meat fraction, are similar (Supplementary Figure 14), except for a steeper slope of log-population density versus NPP in S0 due to the lack of data points under the high end of NPP (as ORCHIDEE underestimates NPP in tropical forests, see Supplementary Figure 11). Regarding the global patterns of population density (Supplementary Figure  13), however, both S0 and S1 predict an unrealistic absence of hunter-gatherers in many arid ecosystems including those in western Plains in North America, Kalahari, and interior Australia, the classical hunter-gatherer territories. These regions are populated in the S2 results (Fig. 3a). The difference between S1 and S2 (for which the annual total NPP in the inputs are identical) therefore indicates the critical role of a substantially long growing season length in sustaining hunter-gatherers in the less productive ecosystems. Considering a prevailing underestimation in growing season length for the arid ecosystems in state-of-the art terrestrial biosphere models (Supplementary Figure 15), our results highlight the need for a more realistic simulation of plant phenology in order to improve the simulation of contemporary and ancient hunter-gatherers. Unless otherwise specified, the model results shown are from the S2 experiment.

Sensitivity tests of parameters
We consider several parameters in FORGE, listed in Supplementary Table 3, to be highly uncertain. We therefore conducted sensitivity tests using Sobol's method 69 (a variance-based sensitivity analysis), which decomposes the variance of model output into fractions that can be attributed to different inputs. Sobol's method has the advantage of measuring sensitivity across the whole input space, as well as accounting for nonlinearity and parameter interactions, so that the total-order index (S T ) indicates the importance of each parameter considering both its main effect (first-order sensitivity index, S 1 ) and the contribution of its interaction with other parameters.
Sobol's method requires thousands of runs, which is computationally expensive for global experiments. We thus carried out the tests at two sites with contrasting characteristics. First, a temperate forest with a long growing season and a low hunting fraction in the diet at 42N, 123W in. Second, a boreal forest with a short growing season and a high hunting fraction at 56N, 69W (Supplementary Figure 16). For each site, we conducted 9000 runs in which the parameters were sampled within their ranges listed in Supplementary  Table 3, using Saltelli's sampling scheme which is more efficient than random sampling. Supplementary Figure 17 shows how the varied parameter values affect the corresponding equations. All the sensitivity tests were done using the Python package "SALib" (https:// salib.readthedocs.io/en/latest/index.html).
The resulting S T and S 1 (Supplementary Figure 18) indicate different sensitivities to the parameters between the two sites. For the hunting-dominated site, population density is most sensitive to Φ A , followed by e max and c, whereas for the gathering-dominated site, population density is sensitive to the four parameters Φ V , Φ A , e max , and c. The three parameters used to calculate birth rate make negligible contributions to the variance in modelled population. Indeed, birth rate impacts the time it takes to reach equilibrium population under a given set of environmental conditions, but has a negligible effect on the equilibrium value itself. Nor is the population density sensitive to the parameter q, which determines the meat craving response.
Supplementary Figure 19 further shows the quantitative response of model results to the parameters. As expected, a higher Φ V (or Φ A ) value increases population density at the gathering-dominated (or hunting-dominated) site, while a higher c value decreases the population at both sites. A higher e max value increases population density at both sites, but the beneficial effect diminishes at high population densities for the hunting-dominated site, probably a result of overhunting. Meat fraction in the diet is mainly determined by the relative values of Φ V and Φ A (Supplementary Figure 18). Compared to population density and meat fraction, foraging time is relatively more equally sensitive to all the tested parameters (Supplementary Figure 18), but its results in the 9000 runs are relatively centralized across the tested parameter space for both sites (Supplementary Figure 19).
For the standard configuration, we set Φ V =0.015, Φ A =0.1, e max =1, and c=200 so as to match the average value of population density and meat fraction across the sites. Note that changing these globally-constant parameter values increases or decreases population everywhere, but has minor impacts on the modelled relationships with environmental variables (Supplementary Figure 20).

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.  Figure 2). The white circles represent median values and the thick (or thin) black bars the inter-quantile (or 5 th -95 th ) ranges. Except for a lower Φ NPP in boreal forest and a slightly higher Φ NPP in Mediterranean forest, Φ NPP does not differ significantly among the other five biome types (ANOVA test p >0.05). FORGE simulates daily energy (carbon) intake from plants and animals and energy expenditure (black arrows); these energy in/out-fluxes update an energetic reserve pool (body fat) every day, which indicates the average health state and then impacts birth rate and mortality rate, determining changes in hunter-gatherer population density every year. The daily intake rate depends on food abundance and on the time spent in hunting and gathering, calculated with time allocation algorithms (see Methods for details).   Model results were averaged over all grid cells of either gathering-dominated (modelled annual mean hunted food <40%, see Fig. 4a) or hunting-dominated (hunted food ≥40%) regions, at dynamic equilibrium. Mean population density of the two regions is 0.19 and 0.04 ind./km 2 , and mean meat fraction of the diet is 22% and 55%, respectively. To synchronize the seasonal cycles between the two hemispheres, outputs of the southern hemisphere were shifted by 6 months. The numbers on the links indicate percentages of annual NPP that flow to herbivores, humans, and litter (dead plants) and soil pools. Note that for legibility, the link widths are not strictly proportional to the magnitude of flows. Spring: MAM; Summer: JJA; Autumn: SON; Winter: DJF.