PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of plosonePLoS OneView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
 
PLoS One. 2017; 12(3): e0172262.
Published online 2017 March 8. doi:  10.1371/journal.pone.0172262
PMCID: PMC5342208

Large-scale effects of migration and conflict in pre-agricultural groups: Insights from a dynamic model

John P. Hart, Editor

Abstract

The debate on the causes of conflict in human societies has deep roots. In particular, the extent of conflict in hunter-gatherer groups remains unclear. Some authors suggest that large-scale violence only arose with the spreading of agriculture and the building of complex societies. To shed light on this issue, we developed a model based on operatorial techniques simulating population-resource dynamics within a two-dimensional lattice, with humans and natural resources interacting in each cell of the lattice. The model outcomes under different conditions were compared with recently available demographic data for prehistoric South America. Only under conditions that include migration among cells and conflict was the model able to consistently reproduce the empirical data at a continental scale. We argue that the interplay between resource competition, migration, and conflict drove the population dynamics of South America after the colonization phase and before the introduction of agriculture. The relation between population and resources indeed emerged as a key factor leading to migration and conflict once the carrying capacity of the environment has been reached.

Introduction

The debate on the causes of conflict in human societies has deep philosophical roots [1], with recent works even proposing it as one of the keys to understand human evolution [2]. However, the actual extent of the conflict in small hunter-gatherer groups remains unclear [3], with some authors suggesting that large-scale warfare and mass-killing only arose with the spreading of agriculture and the building of complex societies, as before resources were too sparse to make their conquest or defense meaningful [36]. Evidence supporting this view points out to the fact that, at least in some regions, the Neolithic revolution was accompanied by “unprecedented levels” of violence [7]. Even if evidence of widespread warfare also exists for previous periods [811], the question of whether the rate of conflict in pre-agricultural societies was sufficient to significantly affect the evolution of human populations on scales larger than the single group remains open [12].

In order to shed light on this issue, we developed a model based on operatorial techniques which simulates population-resource dynamics in a pre-agricultural context at the continental scale. The model outcomes under different conditions were compared with recently available demographic data for prehistoric South America [13]. Using as a proxy for population size the number of occupied archaeological sites and the probability density of summed calibrated radiocarbon dates (SCPDs), the authors reconstructed the spatio-temporal patterns of human population growth in South America, from 14 k to 2 k years ago. This period encompassed three phases: a first phase of rapid expansion from 14 k to 9 k years ago, which represents the continent colonization, followed by approximatively 4 k years of a dynamical equilibrium around a carrying capacity, with no net growth. Finally a new growth phase occurred, due to the rising of agriculture which allowed to sustain a larger population [13, 14].

Our research focuses on the pre-agriculture period and more specifically on the second phase, when the carrying capacity was reached: a condition which is simulated in our model. The demographic curve on that phase is characterized by sharp oscillations of large amplitude (see figure 3 in [13]), a behavior that resembles the dynamics of invasive species. Similar dynamics have indeed been observed in the field. However, the spatial scale of such observations is generally much smaller than the continental one [15]. Consequently, considering the whole South America as a single invaded site does not seem appropriate. From an ecological point of view, it could be better conceptualized as a mosaic of different ecosystems, in some cases separated by geographic barriers. Its population can hence be regarded as a meta-population split into different sub-populations, each of them occupying a patch of the mosaic [16]. This kind of representation appears especially suitable considering hunter-gatherers groups of the mid-Holocene, who normally lived within well defined territories whose size depended on the local ecological conditions [9, 12]. While in each patch of the mosaic we can expect oscillating dynamics between human and resource populations, there is no reason why these oscillations should be in phase on a larger scale, and empirical research confirmed that they often are not [17]. As a consequence, local oscillations should compensate each other, leading to a relatively smooth curve for the global population density.

This intuition was further confirmed by our model, where, consistently with the mosaic concept, South America is represented as a two-dimensional lattice, with humans and resources interacting in each cell of the lattice. The model produced demographic curves that indeed were oscillating at the local level (namely the single cell) but smooth at the global one (namely the whole lattice). The reconstructed curve for the population density in South America instead exhibits high peaks and dramatic drops, from a density of SCPDs less than 0.13 to almost 0.29 around a mean carrying capacity K ≈ 0.185, which cannot be entirely attributed to the calibration process [13]. This suggests that some coordinating mechanism must have been at work to keep in phase the local oscillations.

As a further step, we proposed a possible coordinating mechanism and implemented it into the model. More specifically, we argued that the interplay between resource competition, migration and conflict drove the population dynamics in South America after the colonization phase and before the introduction of agriculture. Unlike the no-migration case, under this second scenario, the model was able to fit remarkably well the reconstructed curve, strongly supporting the correctness of our interpretation. The rest of the paper is organized as follows. The next section illustrates our model, the subsequent one presents the results, and the last one discusses them and derives their general implications.

Methods

Model overview

We modeled South America as a two-dimensional lattice formed by L2 different cells. Each cell includes two interacting populations 𝒫a and 𝒫b, respectively humans and natural resources. Being interested in a pre-agricultural situation, the latter population represents the only sustenance for humans. Starting from some initial conditions, populations evolve following a dynamics similar to the one of a predator-prey system, which from an ecological point of view is the most appropriate to represent a user-resource system with users depending on renewable biogenic resources [18, 19].

We explored two scenarios. In the first, no migration is allowed; in the second, humans are allowed to migrate to neighboring cells, while resources cannot diffuse across cells. Moreover, migration from a cell α to a cell β is assumed to depend on a variable Kα, which reflects the local (i.e., in α) ratio between the densities of resources and humans. More specifically, humans migrate at a rate that increases for decreasing Kα, simulating the fact that migration from α increases when local resources are low. After humans migrate from α to β, resources in α tend to recover, while human density in β, after an initial obvious increment, can subsequently decline. This produces a reduction of the total human populations in α and β, which we interpret as the outcome of increased mortality due to competition and/or conflict over local resources.

In order to reproduce population dynamics, we employed operatorial techniques commonly used in quantum mechanics and already adopted in the analysis of population dynamics and migration [2022]. The details are presented in the following subsections.

One-cell case

Let 𝒮 be a system composed by two populations: the users 𝒫a (i.e. the humans), and the natural resources 𝒫b. The system 𝒮 can assume four basic states, while any other possible state can be obtained as a combination of these: (i) [var phi]0,0, i.e. the vacuum of the system, corresponding to an almost complete absence of both humans and resources; (ii) [var phi]0,1 (resp. [var phi]1,0) corresponding to a low (resp. high) density of humans and abundance (resp. lack) of resources; (iii) [var phi]1,1 abundance of both human and resources.

We associate the four basic states of the system with four independent mutually orthogonal vectors in a four-dimensional Hilbert space endowed with the scalar product left angle bracket[center dot], [center dot]right angle bracket, a procedure that follows the general framework described in details in Bagarello [20] and used in several contexts.

A convenient way to construct the vectors [var phi]j, k with j, k = 0,1 makes use of two fermionic operators, a and b, which together with their adjoints, a and b, satisfy the canonical anticommutation relation (CAR):

{aa} = {bb} = 1 14a2b2 = 0,  {ab} = 0, 
(1)

where {x, y} = xy + yx is the anticommutator of x and y, a can be both a or a, and 1 14 is the identity operator in the Hilbert space ℋ = ℂ4. We first introduce the vacuum vector [var phi]0,0 which satisfies a[var phi]0,0 = b[var phi]0,0 = 0. The existence of a non zero vector [var phi]0,0 with this property is a standard result in quantum mechanics [23]. The other vectors [var phi]j, k can be constructed from [var phi]0,0 in the following way:

φ1,0: = aφ0,0φ0,1: = bφ0,0φ1,1: = abφ0,0.

The set φ = {φj,kjk = 0, 1} is an orthonormal basis for , so that a general vector Ψ of 𝒮 can be expanded as Ψ=j,k=01cj,kφj,k, with j,k=01|cj,k|2=1. By using a standard argument, we can interpret Ψ as a state where the probability to find 𝒮 in the state [var phi]j, k is given by |cj, k|2.

We now introduce the number operators n^(a)=aa and n^(b)=bb of 𝒫a and 𝒫b. Then:

n^(a)φn(a),n(b)=n(a)φn(a),n(b),n^(b)φn(a),n(b)=n(b)φn(a),n(b),
(2)

where n(a), n(b) can be either 0 or 1. It holds:

n^(a)Ψ=k=01c1,kφ1,k,n^(b)Ψ=j=01cj,1φj,1,

and the mean values, or densities, n(a), (b) of the operators n^(a),(b) over the state Ψ are defined as n(a),(b)=Ψ,n^(a),(b)Ψ, leading to the following expressions:

n(a)=k=01|c1,k|2,n(b)=j=01|cj,1|2.
(3)

From a biological point of view, n(a) (resp. n(b)) is interpreted as the density for the population 𝒫a (resp. 𝒫b) associated with the state Ψ. Both n(a), n(b) are real numbers between 0 and 1. For instance, if at t = 0 n(a) = 0, it means that Ψ is only a combination of [var phi]0,0 and [var phi]0,1 corresponding to the (almost complete) absence of humans (zero density for 𝒫a). Conversely, if n(a) = 1 then Ψ is a combination of [var phi]1,0 and [var phi]1,1 corresponding to the maximum density for humans.

2D case

A square lattice divided in L2 cells is introduced in the 2D case. In each cell the two populations 𝒫a and 𝒫b interact as in the previous case.

In each cell α = 1, (...), L2, we introduce the fermionic operators aα, aα and the related number operators n^α(a)=aαaα for 𝒫a and bα, bα and n^α(b)=bαbα for 𝒫b. As before, we suppose that these operators satisfy the standard CAR anticommutation rules:

{aα,aβ}={bα,bβ}=δαβ11,α,β
(4)

{aα,aβ}={bα,bβ}={aα,bβ}=0,αβ.
(5)

Moreover aα2=bα2=0 holds for all α. In Eq (4) 1 1 is the identity operator on the Hilbert space , with dim(ℋ) = 4L2, endowed with the scalar product left angle bracket[center dot], [center dot]right angle bracket.

Extending what we have done before, we first introduce the vacuum vector of the system φ0a,0b, where 0a=0b=(0,0,,0) are two L2-dimensional vectors. The vacuum is annihilated by all the operators aα, bα, i.e.,

aαφ0a,0b=bαφ0a,0b=0,α=1,..,L2.

We then construct the states of the basis of by acting with the operators aα,bα over φ0a,0b,

φma,mb=(a1)ma(1)(aL2)ma(L2)(b1)mb(1)(bL2)mb(L2)φma,mb,
(6)

where ma=(ma(1),,ma(L2)), mb=(mb(1),,mb(L2)) are all the possible L2-dimensional vectors whose entries are only 0 or 1. In fact, aα2=bα2=0, any other choice would destroy the state. Notice that we can construct at most 2L2 vectors ma and 2L2 vectors mb, so that we can build 4L2 possible pairs (ma,mb).

The set φ of all the vectors obtained by this construction forms an orthonormal basis of . Moreover:

n^α(a)φma,mb=ma(α)φma,mb,
(7)

n^α(b)φma,mb=mb(α)φma,mb,
(8)

for all α = 1, (...), L2. For more details on CAR, we refer to Roman [23].

The vectors φma,mb can now be interpreted similarly to [var phi]j, k in the previous section. The vacuum φ0a,0b, for instance, describes a situation where all the lattice cells hold few of both humans and resources. The vector φma,mb with ma=(1,0,0,0) and mb=(0,0,,0,1) describes the same situation, except that there is a large amount of humans in the first cell and of resources in the last one.

A generic state Ψ(0) of the system can be written as a linear combination of the elements in φ:

Ψ(0)=ma,mbcma,mbφma,mb,
(9)

where cma,mb are complex scalars such that ma,mb|cma,mb|2=1.

From Eqs (7) and (8), the φma,mb are eigenstates of the number operators and the densities of 𝒫a and 𝒫b in α are simply their related eigenvalues. This is true at the initial time t = 0. At a later time, we need to compute the mean values over the state Ψ(t) describing the system 𝒮 at time t:

nα(a)(t)=Ψ(t),n^α(a)Ψ(t)=aαΨ(t)2,
(10)

nα(b)(t)=Ψ(t),n^α(b)Ψ(t)=bαΨ(t)2,
(11)

As already stated (see also [20]), these mean values are phenomenologically interpreted as densities of humans and resources in the cell α.

Evolution of the system

To determine the densities Eqs (10) and (11), we need first to compute the time evolution of the state of the system Ψ(t). If Ψ(0)=Ψ0=ma,mbcma,mbφma,mb, then Ψ(t) can be written as Ψ(t)=ma,mbcma,mb(t)φma,mb.

The time-depending functions cma,mb(t) are obtained through the Schrödinger equation iΨ(t)t=HΨ(t), where H is the Hamiltonian operator describing the dynamics of humans and resources. To make the model more realistic, we assume that H explicitly depends on time and on the local densities of 𝒫a and 𝒫b. Using the orthogonality conditions of the basis vectors φma,mb, we obtain the following ordinary differential equation (ODE) system for the coefficients cma,mb(t)

icma,mb(t)t=φma,mb,HΨ(t),
(12)

which, when solved, returns the densities of humans and resources in α:

nα(a)(t)=ma,mbma(α)=1|cma,mb(t)|2,nα(b)(t)=ma,mbmb(α)=1|cma,mb(t)|2.
(13)

Note that the above general construction implies that both nα(a)(0)Ψ(0)2=1 and nα(b)(0)Ψ(0)2=1 at t = 0. As we will see later (see Eq (20)), this condition is preserved for all time, even in the current setting where the Hamiltonian is not purely quadratic and time independent.

It remains now to explicitly define the self—adjoint Hamiltonian of the system. The human-resource interaction is ruled in each cell α by the operator Hα while the migration is ruled by the operator HM, so that the full Hamiltonian is given by: H = (∑α Hα)+HM, where the last term is clearly zero in the first scenario, i.e., in absence of migration. A natural choice for Hα and HM is:

Hα=ωαa(t)aαaα+ωαb(t)bαbα+λα(t)aαbα+bαaα,
(14)

HM=αβpα,βγα,β(t)aαaβ,
(15)

where pα, β ≠ 0 if α and β are neighboring cells, and pα, β = 0 otherwise.

The meaning of the terms Hα, HM has been explained in details and discussed in [21] and [20]. Here we simply recall that the terms in ωαa(t)aαaα and ωαb(t)bαbα are related to a sort of (time-dependent) inertia of the populations within the system, because increasing the values ωαa,b(t) leads to a more static behavior of the associated populations keeping the densities close to their initial values. The terms λα(t)(aαbα+bαaα) are the interaction parts of the Hamiltonian, where the functions λα(t) measure the strength of the interaction between 𝒫a and 𝒫b in α. This interaction leads to a sort of predator-prey dynamics, with the term aαbα that increases the density of 𝒫a whereas the density of 𝒫b decreases. The hermitian conjugate term aαbα induces an opposite effect. The terms aαaβ in HM rule the migration of humans between α and β, because aαaβ increases the density of humans in β while decreasing the density in α. Of course, such a migration is only possible for non zero values of pα, β.

Finally, we have to define the analytic expression for the functions ωαa,b(t), λα(t), and γα, β(t). When migration is allowed, we assume that humans migrate from a cell α to a neighboring one only if the local human/resource ratio Kα(t)=nα(b)(t)/nα(a)(t) is low enough, with migration that intensifies when Kα(t) decreases. Moreover, natural resources tend to proliferate in a given cell in case of low human density, namely if Kα(t) is large, otherwise they tend to exhaust.

In accordance to such assumptions, we chose the following functions:

ωαa(t)=σαexp-(Kα(t)-1/τα)2Kα(t)1/2,
(16)

ωαb(t)=σαexp-(Kα(t)/τα)2Kα(t)-11/2,
(17)

λα(t)=ωαa(t)+ωαb(t)+μα,
(18)

γα,β(t)=ωαb(t)+ωβb(t),
(19)

where σα, τα, μα are real parameters. From Eqs (16)–(19), it follows that the higher the value of Kα(t), the higher (resp. the smaller) the value of ωαa(t) (resp. ωαb(t)). As a consequence, when Kα(t) is high the migratory effects are mainly given by the contribution ωβb(t) in γα, β(t). For instance, when in β the ratio resources/humans is low then the migration occurs mainly from β to α. The function λα(t) in Eq (18) increases for high values of ωαa(t) or ωαb(t), i.e., when Kα(t) is either high or low. This implies that, when human density is low in a cell, it grows at the expense of the resources and vice versa.

The fixed term μα has been introduced to guarantee a minimal predator-prey dynamics even when both ωαa(t) an ωαb(t) are close to zero. Finally, the values of σα, τα tune how strongly Kα(t) affects the system dynamics. In particular, by decreasing τα we limit the interaction effect only to the case of very high or very low values of Kα(t), while increasing σα strengthens it.

Since H = H a general argument shows that

n(a)(t)+n(b)(t)=αnα(a)(t)+nα(b)(t)=n(a)(0)+n(b)(0),
(20)

so that, if the global population density increases, the global resources density decreases and vice-versa. The quantity K=n(a)(0)+n(b)(0)2 can be interpreted as a sort of carrying capacity, since the various terms in the Hamiltonian work to reduce n(a)(t) (and to increase n(b)(t)) when n(a)(t) > K, and vice-versa.

The values of the parameters used in our numerical simulations to reproduce the empirical data were σα = 12.5, τα = 0.35, μα = 0.25, pα, β = 1, ∀α, β. We further fixed K = 0.5. The coefficients in Eq (9) were chosen so that the initial human density nα(a)(0) in each cell randomly differs at most by 10% from K/L2, and with the further global condition n(a)(0) = K. The rationale for this choice was to reproduce the phase when the carrying capacity was reached in the continent. The resources densities were set to nα(b)(0)=(2K/L2-nα(a)(0)).

Notice that ωαa(t) or ωαb(t) can theoretically blow-up when nα(a)=0 or nα(b)=0, i.e., when human or resource local density become zero. However, these conditions never occurred in the simulations, as our initial conditions ensured non zero local densities in each cell. Indeed, even if at some time a local density approaches zero, the interaction and migration terms in H work to avoid any further decreasing. The Matlab code used for the simulation is presented as supporting information.

Results

No-migration scenario

We implemented the model described above with a number of cells varying from 9 to 1000. Although in each cell we observed oscillating dynamics, when considering the global population density—i.e., the density of the population in the whole lattice—the curve became smoother as oscillations compensated each other (Fig 1a). In addition, we observed that this smoothing effect increased with the number of cells included in the model (Fig 1b). As a result, the global density curve became dramatically different from the one for South America between 9 and 5.5 k years BP (see Fig 3b in [13]) suggesting that some coordination mechanism, able to keep in phase the local oscillations, was at work.

Fig 1
Model dynamics in absence of migration.

Migration scenario

In the second scenario, migration was allowed, i.e., humans could move from their cell to a neighboring one, with a possible subsequent decline of population due to increased competition over natural resources (see the Methods section). The model parameters were set to fit the empirical data presented in [13]. Due to computational constraints, simulations of the second scenario are shown here only for an L2 lattice with L = 3. Slightly larger lattices were also considered but no significant difference was observed.

At the scale of the single cell, we observed oscillating dynamics as in the previous scenario, albeit sharper. However, unlike the previous case, oscillations did not compensate each others at the global level and the global population density curve exhibited peaks and drops consistent with the ones empirically estimated for South America (Fig 2). Specifically, we compared the amplitude with respect to the carrying capacity of the oscillations in the empirical curve and in the one produced by the model. This comparison highlighted that under the migration scenario the curves have a very similar behavior: in our model the global population density oscillations fell in the the [0.66K, 1.50K] interval, while in Fig 3b by Goldberg et al. they fell in the [0.70K, 1.57K] interval [13]. In addition, we considered the most relevant peaks and valleys (i.e., the ones almost reaching the interval limits) of both curves, and their number turned to be very similar.

Fig 2
Model dynamics in a nine-cell system with and without migration.

In order to test the model robustness, we performed a sensitivity analysis using a wide range of parameter configurations. Under no condition we were able to reproduce the empirically estimated oscillations in absence of migration and conflict, which hence emerged in our model as a crucial mechanism able both to exacerbate the local oscillations and to keep them synchronized at the global level.

How the various parameters affect the Hamiltonian Eq (15) is well understood from [20], and similar conclusions hold for the present model, with some differences due to the nonlinear parameters that have been introduced here (see the Sensitivity analysis section in the Supporting information). For instance, increasing the value of ωαa(t) produces a stronger inertia for humans in cell α, with their density exhibiting fewer and smaller oscillations. Changing the parameters pα, β, which tunes migration, we obtained different oscillatory behaviors of the population density. In particular, increasing pα, β induced the formation of sharp oscillations in both the global and local density curves, even if this is especially evident at the local (single cell) level (Fig 3). Further details on the sensitivity analysis are included as Supporting information.

Fig 3
Model dynamics for different values of the migration parameters pα, β.

Conclusions

Recent data suggested that pre-agricultural human populations were subjected to significant oscillations on scales larger than the single group and reaching the continental level. Although part of the empirically-estimated oscillations for South America may be an effect of the calibration process, statistical analysis proved them to be significantly different from a random dynamics [13]. Our operator-based model showed that the simple aggregation of local fluctuations, in absence of a more general coordination mechanism, could not reproduce the observed dynamics. A mechanism that our model instead proved to be sufficient to explain the observed dynamics was the synergistic work of migration and conflict.

The picture that can be derived is that, in case of local resource overuse, mobile hunter-gathered groups tried to move to neighboring territories. If these were already occupied, as in South America after the colonization phase, migration likely led to conflict. Both scarcity and conflict hence had the potential to spread in growing circles, with a snowball effect that could reach even the continental scale in a relatively short time.

A major advantage of our interpretation lies on the fact that it is consistent with what we know about hunter-gatherer groups which, as shown by anthropological research, presented significant levels of violence fundamentally caused by competition over resources [9, 10]. In particular, our results are consistent with recent findings showing that reduced environmental productivity is a strong predictor of lethal aggression among prehistoric groups in Central California [24].

The main limit of our work lies in the fact that it has only been checked against a specific phase within a single dataset. New datasets could be used in the future for further testing. North America and Australia look as especially interesting cases as they both share with South America a story of colonization by humans of previously unoccupied territories. Moreover, at the current stage of the research, we only modelled a period of about 4000 years, which is significantly shorter than the time-range of Goldberg and colleagues’ dataset. In future steps, both the colonization phase and the expansion of the carrying capacity linked to the development of agriculture could be included in the model.

Despite these limitations, the current version of our model was already able to shed light on the long enduring debate about the roots of conflict for humans societies [12]. It has formerly been proposed that, due to their low density and the few possessions they owned, conflict among hunter-gatherer groups was limited, if not absent [36]. This “peaceful noble savage” view has already been challenged elsewhere [811, 24], and our model showed that it does perhaps not apply to prehistoric South America as well. More than the lifestyle (hunter-gatherer vs. more complex societies), the relation between population and resources appears to be a key factor eventually leading to migration and conflict once the carrying capacity of a given environment has been reached.

The progressive introduction of agriculture in South America allowed for an increase in the carrying capacity, leading to a new population expansion phase between 5.5 and 2 k years ago [13]. However, also the technological expansion of the possibilities of the natural environment has its limits, and new resource conflicts are likely to arise when they are reached. As planetary boundaries are likely to become a major concern in the near future [25, 26], with a possible increase of environment-related conflict [2729], the lessons that can be derived from the modeling of prehistoric large-scale demographic dynamics could provide valuable insights also for today’s world.

Supporting information

S1 Appendix

Sensitivity analysis.

This appendix presents an extended sensitivity analysis on the model.

(PDF)

S2 Appendix

Model code.

This appendix presents the Matlab code used for the numerical simulation.

(PDF)

S1 Fig

Model dynamics for different values of μα.

(a) Global (whole nine-cell lattice) human population density in the migration scenario. (b) Global human population density in the no-migration scenario. (c) Local (single cell) human population density in the migration scenario, with the central cell selected as example. (d) Local human population density in the no-migration scenario.

(TIF)

S2 Fig

Model dynamics for different values of τα.

(a) Global (whole nine-cell lattice) human population density in the migration scenario. (b) Global human population density in the no-migration scenario. (c) Local (single cell) human population density in the migration scenario, with the central cell selected as example. (d) Local human population density in the no-migration scenario.

(TIF)

S3 Fig

Model dynamics for different values of the migration parameter τα.

(a) Global human population density in the migration scenario. (b) Global human population density in the no-migration scenario.

(TIF)

Acknowledgments

The authors gratefully acknowledge support, comments and suggestions by Torsten Lindström, Andrei Khrennikov, Carlo Corvaja, Tobias Haller, Thomas Campagnaro, Mark W Allen, and two anonymous referees. A special thank goes to Christopher High for his linguistic revision. FB and FG acknowledge technical support from the University of Palermo and GNFM.

Funding Statement

The author(s) received no specific funding for this work.

Data Availability

Data Availability

The code to produce all relevant data is provided as supporting information.

References

1. Lawler A. The Battle Over Violence. Science. 2012;336(6083):829–830. doi: 10.1126/science.336.6083.829 [PubMed]
2. Bowles S. Warriors, Levelers, and the Role of Conflict in Human Social Evolution. Science. 2012;336(6083):876–879. doi: 10.1126/science.1217336 [PubMed]
3. Lawler A. Civilization’s Double-Edged Sword. Science. 2012;336(6083):832–833. doi: 10.1126/science.336.6083.832 [PubMed]
4. Fry DP, editor. War, peace, and human nature: the convergence of evolutionary and cultural views. Oxford: Oxford University Press; 2013.
5. O’Connell RL. Ride of the Second Horseman: The Birth and Death of War. Oxford: Oxford University Press; 1997.
6. Helbling J. The dynamics of war and alliance among the Yanomami. Sociologus. 1999;49:103–118.
7. Meyer C, Lohr C, Gronenborn D, Alt KW. The massacre mass grave of Schöneck-Kilianstädten reveals new insights into collective violence in Early Neolithic Central Europe. Proceedings of the National Academy of Sciences. 2015;. doi: 10.1073/pnas.1504365112 [PubMed]
8. Allen MW, Jones TL, editors. Violence and warfare among hunter-gatherers. Walnut Creek, CA: Left Coast Press; 2014.
9. Gat A. War in human civilization. Oxford: Oxford University Press; 2006.
10. Guilaine J, Zammit J. The origins of war: Violence in prehistory. Malden: Blackwell Publishing; 2005.
11. Sala N, Arsuaga JL, Pantoja-Pérez A, Pablos A, Martínez I, Quam RM, et al. Lethal Interpersonal Violence in the Middle Pleistocene. PLoS One. 2015;10(5):e0126589 doi: 10.1371/journal.pone.0126589 [PMC free article] [PubMed]
12. Bowles S. Did Warfare Among Ancestral Hunter-Gatherers Affect the Evolution of Human Social Behaviors? Science. 2009;324(5932):1293–1298. doi: 10.1126/science.1168112 [PubMed]
13. Goldberg A, Mychajliw AM, Hadly EA. Post-invasion demography of prehistoric humans in South America. Nature. 2016;532(7598):232–235. doi: 10.1038/nature17176 [PubMed]
14. Boserup E. The conditions of agricultural growth: The economics of agrarian change under population pressure. New Brunswick: Transaction Publishers; 2005.
15. Arim M, Abades SR, Neill PE, Lima M, Marquet PA. Spread dynamics of invasive species. Proceedings of the National Academy of Sciences. 2006;103(2):374–378. doi: 10.1073/pnas.0504272102 [PubMed]
16. Wiens JA, Stenseth NC, Horne BV, Ims RA. Ecological Mechanisms and Landscape Ecology. Oikos. 1993;66(3):369 doi: 10.2307/3544931
17. Hanski I. Single-species metapopulation dynamics: concepts, models and observations. Biological Journal of the Linnean Society. 1991;42(1–2):17–38. doi: 10.1111/j.1095-8312.1991.tb00549.x
18. Bravo G, Tamburino L. Are two resources really better than one? Some unexpected results of the availability of substitutes. Journal of Environmental Management. 2011;92(11):2865–2874. doi: 10.1016/j.jenvman.2011.06.010 [PubMed]
19. Perman R, Ma Y, McGilvray J, Common M. Natural resource and environmental economics. Harlow: Pearson Education; 2003.
20. Bagarello F. Quantum Dynamics for Classical Systems: With Applications of the Number Operator. New York: J. Wiley and Sons; 2012.
21. Bagarello F, Oliveri F. A phenomenological operator description of interactions between populations with applications to migration. Math Models Methods Appl Sci. 2013;23(03):471–492. doi: 10.1142/S0218202512500534
22. Bagarello F, Gargano F, Oliveri F. A phenomenological operator description of dynamics of crowds: escape strategies. Applied Mathematical Modelling. 2015;39(8):2276–2294. doi: 10.1016/j.apm.2014.10.038
23. Roman P. Advanced Quantum Mechanics. New York: Addison-Wesley; 1965.
24. Allen MW, Bettinger RL, Codding BF, Jones TL, Schwitalla AW. Resource scarcity drives lethal aggression among prehistoric hunter-gatherers in central California. Proceedings of the National Academy of Sciences. 2016;113(43):12120–12125. doi: 10.1073/pnas.1607996113 [PubMed]
25. Rockström J, Steffen W, Noone K, Persson Åsa, Chapin FS, Lambin EF, et al. A safe operating space for humanity. Nature. 2009;461:472–475. doi: 10.1038/461472a [PubMed]
26. Steffen W, Richardson K, Rockström J, Cornell SE, Fetzer I, Bennett EM, et al. Planetary boundaries: Guiding human development on a changing planet. Science. 2015;347(6223):736 doi: 10.1126/science.1259855 [PubMed]
27. Frank AB, Collins MG, Levin SA, Lo AW, Ramo J, Dieckmann U, et al. Dealing with femtorisks in international relations. Proceedings of the National Academy of Sciences. 2014;111(49):17356–17362. doi: 10.1073/pnas.1400229111 [PubMed]
28. Homer-Dixon TF. Environment, scarcity, and violence. Princeton: Princeton University Press; 1999.
29. Schleussner CF, Donges JF, Donner RV, Schellnhuber HJ. Armed-conflict risks enhanced by climate-related disasters in ethnically fractionalized countries. Proceedings of the National Academy of Sciences. 2016;113(33):9216–9221. doi: 10.1073/pnas.1601611113 [PubMed]

Articles from PLoS ONE are provided here courtesy of Public Library of Science