Home | About | Journals | Submit | Contact Us | Français |

**|**PLoS One**|**v.12(3); 2017**|**PMC5342208

Formats

Article sections

Authors

Related links

PLoS One. 2017; 12(3): e0172262.

Published online 2017 March 8. doi: 10.1371/journal.pone.0172262

PMCID: PMC5342208

Francesco Gargano,^{1} Lucia Tamburino,^{2} Fabio Bagarello,^{1,}^{3} and Giangiacomo Bravo^{4,}^{5,}^{*}

John P. Hart, Editor^{}

New York State Museum, UNITED STATES

**Conceptualization:**GB LT FB.**Data curation:**FG .**Formal analysis:**FG.**Investigation:**FG.**Methodology:**FB FG.**Software:**FG.**Validation:**FG.**Visualization:**FG.**Writing – original draft:**LT GB FG FB.**Writing – review & editing:**LT GB FG FB.

* E-mail: es.unl@ovarb.omocaignaig

Received 2016 October 3; Accepted 2017 February 2.

Copyright © 2017 Gargano et al

This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

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.

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 [3–6]. 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 [8–11], 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.

We modeled South America as a two-dimensional lattice formed by *L*^{2} 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 [20–22]. The details are presented in the following subsections.

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)*
_{0,0}, i.e. the *vacuum* of the system, corresponding to an almost complete absence of both humans and resources; *(ii)*
_{0,1} (resp. _{1,0}) corresponding to a low (resp. high) density of humans and abundance (resp. lack) of resources; *(iii)*
_{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 , , 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 _{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):

{*a*, *a*^{†}} = {*b*, *b*^{†}} = 1 1_{4}, *a*^{2} = *b*^{2} = 0, {*a*^{♯}, *b*^{♯}} = 0,

(1)

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

The set ℱ_{φ} = {*φ*_{j,k}, *j*, *k* = 0, 1} is an orthonormal basis for ℋ, so that a general vector Ψ of 𝒮 can be expanded as $\Psi ={\displaystyle \sum _{j,k=0}^{1}{c}_{j,k}{\phi}_{j,k}}$, with $\sum _{j,k=0}^{1}{\left|{c}_{j,k}\right|}^{2}=1$. By using a standard argument, we can interpret Ψ as a state where the probability to find 𝒮 in the state _{j, k} is given by |*c*_{j, k}|^{2}.

We now introduce the number operators ${\widehat{n}}^{\left(a\right)}={a}^{\u2020}a$ and ${\widehat{n}}^{\left(b\right)}={b}^{\u2020}b$ of 𝒫_{a} and 𝒫_{b}. Then:

$$\begin{array}{l}\begin{array}{ccc}{\widehat{n}}^{\left(a\right)}{\phi}_{{n}^{\left(a\right)},{n}^{\left(b\right)}}& =& {n}^{\left(a\right)}{\phi}_{{n}^{\left(a\right)},{n}^{\left(b\right)}},\end{array}\\ \begin{array}{ccc}{\widehat{n}}^{\left(b\right)}{\phi}_{{n}^{\left(a\right)},{n}^{\left(b\right)}}& =& {n}^{\left(b\right)}{\phi}_{{n}^{\left(a\right)},{n}^{\left(b\right)}},\end{array}\end{array}$$

(2)

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

$$\begin{array}{c}\hfill {\widehat{n}}^{\left(a\right)}\Psi =\sum _{k=0}^{1}{c}_{1,k}{\phi}_{1,k},\phantom{\rule{1.em}{0ex}}{\widehat{n}}^{\left(b\right)}\Psi =\sum _{j=0}^{1}{c}_{j,1}{\phi}_{j,1},\end{array}$$

and the mean values, or densities, *n*^{(a), (b)} of the operators ${\widehat{n}}^{\left(a\right),\left(b\right)}$ over the state Ψ are defined as ${n}^{\left(a\right),\left(b\right)}=\langle \Psi ,{\widehat{n}}^{\left(a\right),\left(b\right)}\Psi \rangle $, leading to the following expressions:

$$\begin{array}{c}\hfill {n}^{\left(a\right)}=\sum _{k=0}^{1}|{c}_{1,k}{|}^{2},\phantom{\rule{1.em}{0ex}}{n}^{\left(b\right)}=\sum _{j=0}^{1}{\left|{c}_{j,1}\right|}^{2}.\end{array}$$

(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 _{0,0} and _{0,1} corresponding to the (almost complete) absence of humans (zero density for 𝒫_{a}). Conversely, if *n*^{(a)} = 1 then Ψ is a combination of _{1,0} and _{1,1} corresponding to the maximum density for humans.

A square lattice divided in *L*^{2} 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, , *L*^{2}, we introduce the fermionic operators *a*_{α}, ${a}_{\alpha}^{\u2020}$ and the related number operators ${\widehat{n}}_{\alpha}^{\left(a\right)}={a}_{\alpha}^{\u2020}{a}_{\alpha}$ for 𝒫_{a} and *b*_{α}, ${b}_{\alpha}^{\u2020}$ and ${\widehat{n}}_{\alpha}^{\left(b\right)}={b}_{\alpha}^{\u2020}{b}_{\alpha}$ for 𝒫_{b}. As before, we suppose that these operators satisfy the standard CAR anticommutation rules:

$$\begin{array}{ccc}\hfill \{{a}_{\alpha},{a}_{\beta}^{\u2020}\}& =& \{{b}_{\alpha},{b}_{\beta}^{\u2020}\}={\delta}_{\alpha \beta}1\phantom{\rule{-0.40em}{0ex}}1,\phantom{\rule{1.em}{0ex}}\forall \alpha ,\beta \hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill \{{a}_{\alpha}^{\u266f},{a}_{\beta}^{\u266f}\}& =& \{{b}_{\alpha}^{\u266f},{b}_{\beta}^{\u266f}\}=\{{a}_{\alpha}^{\u266f},{b}_{\beta}^{\u266f}\}=0,\phantom{\rule{1.em}{0ex}}\forall \alpha \ne \beta .\hfill \end{array}$$

(5)

Moreover ${a}_{\alpha}^{2}={b}_{\alpha}^{2}=0$ holds for all *α*. In Eq (4)
1 1 is the identity operator on the Hilbert space ℋ, with dim(ℋ) = 4^{L2}, endowed with the scalar product , .

Extending what we have done before, we first introduce the *vacuum* vector of the system ${\phi}_{{\overrightarrow{0}}_{a},{\overrightarrow{0}}_{b}}$, where ${\overrightarrow{0}}_{a}={\overrightarrow{0}}_{b}=(0,0,\cdots ,0)$ are two *L*^{2}-dimensional vectors. The *vacuum* is annihilated by all the operators *a*_{α}, *b*_{α}, i.e.,

$$\begin{array}{c}\hfill {a}_{\alpha}{\phi}_{{\overrightarrow{0}}_{a},{\overrightarrow{0}}_{b}}={b}_{\alpha}{\phi}_{{\overrightarrow{0}}_{a},{\overrightarrow{0}}_{b}}=0,\phantom{\rule{1.em}{0ex}}\forall \alpha =1,..,{L}^{2}.\end{array}$$

We then construct the states of the basis of ℋ by acting with the operators ${a}_{\alpha}^{\u2020},{b}_{\alpha}^{\u2020}$ over ${\phi}_{{\overrightarrow{0}}_{a},{\overrightarrow{0}}_{b}}$,

$$\begin{array}{c}\hfill {\phi}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}={\left({a}_{1}^{\u2020}\right)}^{{m}_{a}\left(1\right)}\cdots {\left({a}_{{L}^{2}}^{\u2020}\right)}^{{m}_{a}\left({L}^{2}\right)}{\left({b}_{1}^{\u2020}\right)}^{{m}_{b}\left(1\right)}\cdots {\left({b}_{{L}^{2}}^{\u2020}\right)}^{{m}_{b}\left({L}^{2}\right)}{\phi}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}},\end{array}$$

(6)

where ${\overrightarrow{m}}_{a}=({m}_{a}\left(1\right),\cdots ,{m}_{a}\left({L}^{2}\right))$, ${\overrightarrow{m}}_{b}=({m}_{b}\left(1\right),\cdots ,{m}_{b}\left({L}^{2}\right))$ are all the possible *L*^{2}-dimensional vectors whose entries are only 0 or 1. In fact, ${{a}_{\alpha}^{\u2020}}^{2}={{b}_{\alpha}^{\u2020}}^{2}=0$, any other choice would destroy the state. Notice that we can construct at most 2^{L2} vectors ${\overrightarrow{m}}_{a}$ and 2^{L2} vectors ${\overrightarrow{m}}_{b}$, so that we can build 4^{L2} possible pairs $({\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b})$.

The set ℱ_{φ} of all the vectors obtained by this construction forms an orthonormal basis of ℋ. Moreover:

$$\begin{array}{ccc}\hfill {\widehat{n}}_{\alpha}^{\left(a\right)}{\phi}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}& =& {m}_{a}\left(\alpha \right){\phi}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}},\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill {\widehat{n}}_{\alpha}^{\left(b\right)}{\phi}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}& =& {m}_{b}\left(\alpha \right){\phi}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}},\hfill \end{array}$$

(8)

for all *α* = 1, , *L*^{2}. For more details on CAR, we refer to Roman [23].

The vectors ${\phi}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}$ can now be interpreted similarly to _{j, k} in the previous section. The vacuum ${\phi}_{{\overrightarrow{0}}_{a},{\overrightarrow{0}}_{b}}$, for instance, describes a situation where all the lattice cells hold few of both humans and resources. The vector ${\phi}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}$ with ${\overrightarrow{m}}_{a}=(1,0,0\cdots ,0)$ and ${\overrightarrow{m}}_{b}=(0,0,\cdots ,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 ℱ_{φ}:

$$\begin{array}{c}\hfill \Psi \left(0\right)=\sum _{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}{c}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}{\phi}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}},\end{array}$$

(9)

where ${c}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}$ are complex scalars such that $\sum _{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}{\left|{c}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}\right|}^{2}=1$.

From Eqs (7) and (8), the ${\phi}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}$ 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*:

$$\begin{array}{ccc}\hfill {n}_{\alpha}^{\left(a\right)}\left(t\right)& =& \u27e8\Psi \left(t\right),{\widehat{n}}_{\alpha}^{\left(a\right)}\Psi \left(t\right)\u27e9={\parallel {a}_{\alpha}\Psi \left(t\right)\parallel}^{2},\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill {n}_{\alpha}^{\left(b\right)}\left(t\right)& =& \u27e8\Psi \left(t\right),{\widehat{n}}_{\alpha}^{\left(b\right)}\Psi \left(t\right)\u27e9={\parallel {b}_{\alpha}\Psi \left(t\right)\parallel}^{2},\hfill \end{array}$$

(11)

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

To determine the densities Eqs (10) and (11), we need first to compute the time evolution of the state of the system Ψ(*t*). If $\Psi \left(0\right)={\Psi}_{0}={\displaystyle \sum _{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}{c}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}{\phi}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}}$, then Ψ(*t*) can be written as $\Psi \left(t\right)={\displaystyle \sum _{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}{c}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}\left(t\right){\phi}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}}$.

The time-depending functions ${c}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}\left(t\right)$ are obtained through the Schrödinger equation $i\frac{\partial \Psi \left(t\right)}{\partial t}=H\phantom{\rule{0.166667em}{0ex}}\Psi \left(t\right)$, 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 ${\phi}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}$, we obtain the following ordinary differential equation (ODE) system for the coefficients ${c}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}\left(t\right)$

$$\begin{array}{c}\hfill i\frac{\partial {c}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}\left(t\right)}{\partial t}=\u27e8{\phi}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}},H\Psi \left(t\right)\u27e9,\end{array}$$

(12)

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

$$\begin{array}{l}{n}_{\alpha}^{\left(a\right)}\left(t\right)={\displaystyle \sum _{\begin{array}{c}{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}\\ {m}_{a}\left(\alpha \right)=1\end{array}}{\left|{c}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}\left(t\right)\right|}^{2}},\hfill \\ {n}_{\alpha}^{\left(b\right)}\left(t\right)={\displaystyle \sum _{\begin{array}{c}{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}\\ {m}_{b}\left(\alpha \right)=1\end{array}}{\left|{c}_{{\overrightarrow{m}}_{a},{\overrightarrow{m}}_{b}}\left(t\right)\right|}^{2}}.\hfill \end{array}$$

(13)

Note that the above general construction implies that both ${n}_{\alpha}^{\left(a\right)}\left(0\right)\le {\parallel \Psi \left(0\right)\parallel}^{2}=1$ and ${n}_{\alpha}^{\left(b\right)}\left(0\right)\le {\parallel \Psi \left(0\right)\parallel}^{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 *H*_{M}, so that the full Hamiltonian is given by: *H* = (∑_{α}
*H*_{α})+*H*_{M}, where the last term is clearly zero in the first scenario, i.e., in absence of migration. A natural choice for *H*_{α} and *H*_{M} is:

$$\begin{array}{ccc}\hfill {H}_{\alpha}& =& {\omega}_{\alpha}^{a}\left(t\right){a}_{\alpha}^{\u2020}{a}_{\alpha}+{\omega}_{\alpha}^{b}\left(t\right){b}_{\alpha}^{\u2020}{b}_{\alpha}+{\lambda}_{\alpha}\left(t\right)\left({a}_{\alpha}^{\u2020}{b}_{\alpha}+{b}_{\alpha}^{\u2020}{a}_{\alpha}\right),\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill {H}_{M}& =& \sum _{\alpha \ne \beta}{p}_{\alpha ,\beta}{\gamma}_{\alpha ,\beta}\left(t\right){a}_{\alpha}{a}_{\beta}^{\u2020},\hfill \end{array}$$

(15)

where *p*_{α, β} ≠ 0 if *α* and *β* are neighboring cells, and *p*_{α, β} = 0 otherwise.

The meaning of the terms *H*_{α}, *H*_{M} has been explained in details and discussed in [21] and [20]. Here we simply recall that the terms in ${\omega}_{\alpha}^{a}\left(t\right){a}_{\alpha}^{\u2020}{a}_{\alpha}$ and ${\omega}_{\alpha}^{b}\left(t\right){b}_{\alpha}^{\u2020}{b}_{\alpha}$ are related to a sort of (time-dependent) *inertia* of the populations within the system, because increasing the values ${\omega}_{\alpha}^{a,b}\left(t\right)$ leads to a more static behavior of the associated populations keeping the densities close to their initial values. The terms ${\lambda}_{\alpha}\left(t\right)({a}_{\alpha}^{\u2020}{b}_{\alpha}+{b}_{\alpha}^{\u2020}{a}_{\alpha})$ 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}_{\alpha}^{\u2020}{b}_{\alpha}$ that increases the density of 𝒫_{a} whereas the density of 𝒫_{b} decreases. The hermitian conjugate term ${a}_{\alpha}^{\u2020}{b}_{\alpha}$ induces an opposite effect. The terms ${a}_{\alpha}{a}_{\beta}^{\u2020}$ in *H*_{M} rule the migration of humans between *α* and *β*, because ${a}_{\alpha}{a}_{\beta}^{\u2020}$ 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 ${\omega}_{\alpha}^{a,b}\left(t\right)$, *λ*_{α}(*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}_{\alpha}\left(t\right)={n}_{\alpha}^{\left(b\right)}\left(t\right)/{n}_{\alpha}^{\left(a\right)}\left(t\right)$ 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:

$$\begin{array}{ccc}\hfill {\omega}_{\alpha}^{a}\left(t\right)& =& {\sigma}_{\alpha}{\left({exp}^{-{({K}_{\alpha}{\left(t\right)}^{-1}/{\tau}_{\alpha})}^{2}}{K}_{\alpha}\left(t\right)\right)}^{1/2},\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill {\omega}_{\alpha}^{b}\left(t\right)& =& {\sigma}_{\alpha}{\left({exp}^{-{({K}_{\alpha}\left(t\right)/{\tau}_{\alpha})}^{2}}{K}_{\alpha}{\left(t\right)}^{-1}\right)}^{1/2},\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill {\lambda}_{\alpha}\left(t\right)& =& {\omega}_{\alpha}^{a}\left(t\right)+{\omega}_{\alpha}^{b}\left(t\right)+{\mu}_{\alpha},\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill {\gamma}_{\alpha ,\beta}\left(t\right)& =& {\omega}_{\alpha}^{b}\left(t\right)+{\omega}_{\beta}^{b}\left(t\right),\hfill \end{array}$$

(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 ${\omega}_{\alpha}^{a}\left(t\right)$ (resp. ${\omega}_{\alpha}^{b}\left(t\right)$). As a consequence, when *K*_{α}(*t*) is high the migratory effects are mainly given by the contribution ${\omega}_{\beta}^{b}\left(t\right)$ 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 ${\omega}_{\alpha}^{a}\left(t\right)$ or ${\omega}_{\alpha}^{b}\left(t\right)$, 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 ${\omega}_{\alpha}^{a}\left(t\right)$ an ${\omega}_{\alpha}^{b}\left(t\right)$ 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

$$\begin{array}{c}\hfill {n}^{\left(a\right)}\left(t\right)+{n}^{\left(b\right)}\left(t\right)=\sum _{\alpha}\left({n}_{\alpha}^{\left(a\right)}\left(t\right)+{n}_{\alpha}^{\left(b\right)}\left(t\right)\right)={n}^{\left(a\right)}\left(0\right)+{n}^{\left(b\right)}\left(0\right),\end{array}$$

(20)

so that, if the global population density increases, the global resources density decreases and vice-versa. The quantity $K=\frac{{n}^{\left(a\right)}\left(0\right)+{n}^{\left(b\right)}\left(0\right)}{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}_{\alpha}^{\left(a\right)}\left(0\right)$ in each cell randomly differs at most by 10% from *K*/*L*^{2}, 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}_{\alpha}^{\left(b\right)}\left(0\right)=(2K/{L}^{2}-{n}_{\alpha}^{\left(a\right)}\left(0\right))$.

Notice that ${\omega}_{\alpha}^{a}\left(t\right)$ or ${\omega}_{\alpha}^{b}\left(t\right)$ can theoretically blow-up when ${n}_{\alpha}^{\left(a\right)}=0$ or ${n}_{\alpha}^{\left(b\right)}=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.

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.

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 *L*^{2} 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.66*K*, 1.50*K*] interval, while in Fig 3b by Goldberg *et al*. they fell in the [0.70*K*, 1.57*K*] 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.

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 ${\omega}_{\alpha}^{a}\left(t\right)$ 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.

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 [3–6]. This “peaceful noble savage” view has already been challenged elsewhere [8–11, 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 [27–29], the lessons that can be derived from the modeling of prehistoric large-scale demographic dynamics could provide valuable insights also for today’s world.

This appendix presents an extended sensitivity analysis on the model.

(PDF)

Click here for additional data file.^{(957K, pdf)}

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

(PDF)

Click here for additional data file.^{(158K, pdf)}

(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)

Click here for additional data file.^{(644K, tif)}

(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)

Click here for additional data file.^{(705K, tif)}

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

(TIF)

Click here for additional data file.^{(261K, tif)}

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.

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

Data Availability

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

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**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |