Human-Plant Coevolution: A modelling framework for theory-building on the origins of agriculture

The domestication of plants and the origin of agricultural societies has been the focus of much theoretical discussion on why, how, when, and where these happened. The ‘when’ and ‘where’ have been substantially addressed by different branches of archaeology, thanks to advances in methodology and the broadening of the geographical and chronological scope of evidence. However, the ‘why’ and ‘how’ have lagged behind, holding on to relatively old models with limited explanatory power. Armed with the evidence now available, we can return to theory by revisiting the mechanisms allegedly involved, disentangling their connection to the diversity of trajectories, and identifying the weight and role of the parameters involved. We present the Human-Plant Coevolution (HPC) model, which represents the dynamics of coevolution between a human and a plant population. The model consists of an ecological positive feedback system (mutualism), which can be reinforced by positive evolutionary feedback (coevolution). The model formulation is the result of wiring together relatively simple simulation models of population ecology and evolution, through a computational implementation in R. The HPC model captures a variety of potential scenarios, though which conditions are linked to the degree and timing of population change and the intensity of selective pressures. Our results confirm that the possible trajectories leading to neolithisation are diverse and involve multiple factors. However, simulations also show how some of those factors are entangled, what are their effects on human and plant populations under different conditions, and what might be the main causes fostering agriculture and domestication.


Introduction
The domestication of plants and the origin of agriculture is a major change in human history, and it has been the focus of much theoretical discussion on why, how, when and where this change happened. Evidence from archaeobotany and plant genomics gathered during the last two decades expanded our knowledge on where this process happened and identified several centres of agricultural origin around the world [1][2][3]. Methodological advances in identification criteria [4] and the widespread recovery of plant remains from archaeological sites [5] better clarified the timing of this process in many areas. However, a better understanding of the why and how agriculture began seems to be still elusive [6][7][8].
Climate change [9][10][11], cognitive/symbolic change [12][13][14], or social competition and demography [15,16] have long been discussed as drivers for socio-ecological transformations called the Neolithic Revolution [17]. A major problem with these approaches is to bundle under the same explanation behavioural trajectories that do not necessarily share the same premises. Domestication and agriculture emerged from diverse historical contexts and the empirical record available is manifold, inherently biased and fragmentary due to preservation issues, and it can often also be contradictory in evidencing causality [18]. Furthermore, several models rely on ethnographic observations of contemporary traditional practices among indigenous peoples around the world [19][20][21][22][23]. While these practices make a useful basis for creating models of the past, they may greatly differ in context from those of the first communities engaging in agriculture within any given region, and therefore such "parallelisms" need to be used with care [24].
A current and lively discourse on how domestication (and eventually agriculture) came into being is that of protracted [25][26][27][28] versus expedite [14,29] domestication. Broad contextual analyses of the archaeobotanical record within macroevolutionary theory [18] and singlecrop approaches [30] started to bring new light on the process of domestication based on a fast-growing body of archaeological evidence. The analysis of this massive and relatively recent volume of data makes clear that it is now necessary to return to theory by revisiting the mechanisms allegedly involved in domestication, disentangling their connection to a diversity of trajectories [31,32], being those protracted or sudden, and identifying the weight of the social and ecological parameters. Approaches developed within human behavioural ecology [7,[33][34][35][36][37][38], such as niche construction or cultural niche construction theories, have gained momentum in this effort. These approaches emphasise "the capacity of organisms to modify natural selection in their environment and thereby act as co-directors of their own, and other species' evolution" [39]. However, such perspectives have been heavily criticised, especially as they are considered by some researchers indifferent to the role of human agency and intentionality [14,29,40,41]. The relevant, yet stale, century-long debate on human intentionality in plant domestication is a clear sign that the field still lacks a unifying theoretical framework.

Simulation approaches to human-plant coevolution
The study of the prehistoric human past is necessarily approached through the archaeological record, which does not always allow addressing historical processes and organizational dynamics. Information gaps as well as uncertainty in the record are behind the push for archaeology to participate in and take advantage of innovative methodological approaches, such as modelling and simulation. In a subject like domestication and the origins of agriculture, where the archaeological record is incomplete in both space and time, and real-world experiments are unrealistic, the use of modelling and simulation has become a useful alternative for testing hypotheses and building theory [42]. However, the most important contributions within this framework have focus on the representation of plant domestication in terms of genetic change [25,43] and the geographical spread of the Neolithic transition [44][45][46][47], mainly for testing hypotheses related to regional or species-wise case studies. Exceptionally, there have been key contributions from niche construction and optimal foraging theory as well as complex adaptative systems, but such contributions have been mostly centred on the human side of the process [31,38,[48][49][50]. Few simulation models have considered coevolution as the core mechanism producing changes in both plants and humans [51,52], while the first proposals in this line date back to almost forty years ago [53].
The current work explores hypotheses on plant domestication and the origin of agriculture by using a coevolutionary framework capable of accounting for both plant and human factors. Our model combines readily-available formal models for mutualism and evolution used in population ecology, sociology and economics. Despite sharing the term "coevolution", our approach is neither based on nor necessarily aligned with the gene-culture coevolution or dual inheritance theory. The latter concerns a coupled process of genetic and cultural change in the same population and species, typically humans and other primates, in which other populations and species, and their changes, are considered as factors rather than the subjects of coevolution [54]. Likewise, the model we propose can be distinguished from human behaviour ecology models in this field since these have been defined in terms of human behaviour only (e.g., focusing on decision-making criteria) while factoring other species primarily as resources [38,55]. We state our model assumptions explicitly and have worked intensively on documenting all implementation details to assure its reproducibility and facilitate re-use and future expansions. Our contribution is theoretical and explorative, thus it is not driven by the use of any specific dataset or case study. Furthermore, it does not carry the pretence -at least in its current form-of direct applicability to the many formats of empirical data.

The Human-Plant Coevolution (HPC) model
Human-plant interaction is a specific case of animal-plant interaction, which spans predatorprey, mutualistic and symbiotic relationships. All ecological relationships consistent in time are driven by coevolution, where each party exerts selective pressures on the other, eventually redefining their genetic (and cultural) construct [53,[56][57][58]. Under mutualistic coevolution, the interaction between two populations increases the total potential return or carrying capacity of the environment for each species. At the same time, it also modifies the selective pressures acting over the populations involved. In this light, plant domestication is similar to other mutualistic relationships, where coevolution made possible the emergence of certain traits, manifested at physiological, morphological and behavioural levels; e.g., insects and fungi [59] and ants and acacias [60].
The Human-Plant Coevolution (HPC) model is an ecological positive feedback system (mutualism), which can be reinforced by an evolutionary positive feedback (coevolution). The model is the result of wiring together relatively simple models of population ecology (Verhulst-Pearl model) and evolution (replicator dynamics), through a computational implementation using R programming language [61].
The HPC model embodies the dynamics of two interacting populations: one of humans and another of a given plant species. Here, we assume that population units are individual organisms. Because this model greatly simplifies the mechanisms involved in population dynamics, units could also be set to be groups of individuals or even population proxies (e.g. human working hours, plant-covered soil surface). However, the scale of population units is relevant when calibrating the parameters and interpreting results, and thus must be made explicit.
Each population unit may exploit the available resources in different ways, and may have a different utility for sustaining the other population. To represent this, we assume that each population can be divided into types ranging from the least (1) to the most (n) mutualistic, each corresponding to a value of baseline carrying capacity and utility per capita, which in turn range from population-specific minima and maxima. Each type can relate either to truly discrete units (e.g., presence/absence of trait), arbitrary degrees in a continuum (e.g., size of anatomy trait, frequency of behaviour), or a combination of both. In the case of human populations, types would consist majorly of different combinations of behaviours impacting the plant population, such as protection from predators, removal of competitors, enhancement of soil conditions, or transporting and storing propagules.
This simplification of population diversity gives the possibility to implement a relatively simple and straightforward mechanism of evolution, the replicator dynamics [62]. Under our specific version of this mechanism, the distribution of a population within types changes depending on three factors: undirected variation, inertia, and selection.
The HPC model was conceptualised as a highly symmetric structure (Fig 1). This model reduces the complexity of the human and plant populations to a point where these can be defined using the same terms (parameters and variables). The symmetry is only broken by the inclusion of a constraint specific to plants, the maximum number of plant units fitting  Tables 1 and 2). Populations are shown in yellow, their change in red, type-wise vector or array variables in blue, aggregate population variables in orange, and parameters in white.
https://doi.org/10.1371/journal.pone.0260904.g001 the area available (MaxArea or max_area), reflecting one of the main ecological differences between plants and animals: the latter are able to move and exploit multiple habitats within a lifetime.
The HPC model enables to reproduce a double positive feedback loop, where two populations increase their carrying capacity (mutualism) and empower this relationship by influencing each other's trait selection (coevolution). The consequence is that, given certain conditions, both human and plant populations shift to stronger mutualism types and increase their numbers, potentially moving far away from pre-coevolutionary levels (Fig 2).
All parameters and variables of the model are listed and defined in Tables 1 and 2, respectively. States of the system are evaluated and compared by a set of output variables, i.e. those not used to recalculate the state of the system (Table 3). Among the output variables, the coevolution coefficients are the most revealing. Each indicates if and how much the population type distribution has been modified by the coevolutionary process. Their values range between -1 (the entire population is of type 1) and 1 (the entire population is of type n).

Ecological relationships and population dynamics
The model can be expressed by a relatively simple system of two discrete-time difference equations Eq (1), based on the Verhulst-Pearl Logistic equation [63,64]. The change of both  Table 2) depends on an intrinsic growth rate (r H or intrinsic_grow-th_rate_humans, r P or intrinsic_growth_rate_plants), the population at a given time (H[t] or humans, P[t] or plants) and the respective carrying capacity of the  , which may also vary over time.
Human and plant populations engage in a mutualistic relationship, where one species is to some extent sustained by the other Eq Both populations are also sustained by an independent term, representing the baseline carrying capacity of the environment or the utility gain from other resources, which is timedependent (U bH [t] or utility_other_to_humans, U bP [t] or utility_other_-to_plants). While assuming that the growth of the human population has no predefined ceiling, the expansion of the plant population is considered limited by the area over which plants can grow contiguously (MaxArea or max_area), and represented as a compendium of both space and the maximum energy available in a discrete location Eq (2a).
Considering that mutualistic relationships involve a positive feedback loop, the population growth at time t improves the conditions for both humans and plants at time t + 1, sustaining their growth even further. See model assumptions in Table 4.
Population diversity. The HPC model contemplates a vector pop of length n for each population, containing the population fractions of each type (pop H [t] or type_propor-tions_humans, pop P [t] or type_proportions_plants). The lengths of these vectors or the numbers of types are population-specific and given as two parameters (n H or num-ber_types_humans, n P or number_types_plants). These vectors include all possible variations within a population so that they each amount to unity (i.e. P n i¼1 pop H i ¼ 1 and P n i¼1 pop P i ¼ 1).
To account for multiple types, we replace Eq (3) with Eq (4), where the utility of one population to the other at any given time (U HP [t] [t] or utility_other_to_plants) are calculated similarly, though using the utility that each type is able to gain from other resources (U bHi or type_utili-ty_other_to_humans, U bPi or type_utility_other_to_plants) Eq (5).

Domains Assumptions
On interacting populations A population of humans interacts with a population of plants.
On population growth Population growth is a self-catalysing process, where the population density in the present will contribute to its own increase in the future, depending on an intrinsic growth rate (r). Population growth is a self-limiting process, where the population density in the present will constraint its own increase in the future, depending on respective carrying capacity of the environment (K). The logistic growth model is acceptable as an approximation to the dynamics of populations, both human and plant, under constant conditions. The carrying capacity of the environment for a population depends on constant factors and on a time-varying factor (K[t]).

On positive ecological relationships
Positive ecological relationships exist, where an individual of one population increases by an amount the carrying capacity of the environment for another population.
Coupled positive ecological relationships (i.e., mutualism) exist, where two populations increase the carrying capacities for each other. There is variation in positive ecological relationships, so individuals of one population vary in terms of how much they increase the carrying capacity for the other population.
On human-plant mutualism A given plant species yield a positive utility for humans, e.g., as a source of food and raw materials. Humans return a positive utility for this plant species, e.g., by improving soil conditions. The utility given by one population adds value to the carrying capacity for the other, and vice versa. The carrying capacity for humans rely also on other resources, which are independent of the plant species (i.e., the baseline carrying capacity for humans). The carrying capacity for plants also relies on other conditions, which are independent of humans (i.e., the baseline carrying capacity for plants). The carrying capacity for plants is eventually constrained by the space available for it to grow contiguously as a population (i.e., maximum area). https://doi.org/10.1371/journal.pone.0260904.t004 The shares of population within types follow a one-tail distribution rather than a normal distribution, which would be more adequate but less straightforward to use in a theoretical model. Under this circumstance, the distribution of population within types will always be biased towards the intermediate types. ) and the expected proportion per type, assuming a uniform distribution among types (1/n H and 1/n P ) Eq (6).
By considering inertia as an evolutionary mechanism, we assume that the more frequent a type is, the more likely that it is transmitted. Selection is implemented by assigning a fitness score to each type (fitness Hi [t] or fitness_humans, fitness Pi [t] or fitness_plants), which in turn biases its transmission. Eq (7) [65].
The replicator dynamics described so far defines how a trait evolves in a single population. However, coevolution can also be represented when the selective pressure on one population is modified by the changing traits of another population. In order to link two populations, the fitness scores of one population are derived from the weight of the contribution or utility of the other population (U HP [t] or utility_humans_to_plants, U PH [t] or utility_-plants_to_humans) in relation to the base carrying capacity for the former (U bH [t] or utility_other_to_humans, U bP [t] or utility_other_to_plants) Eq (8).
As a consequence of the model design, types of both human and plant populations span from a non-mutualistic type (i = 1), which has the best fitness score when there is no positive interaction with the other population (e.g., type 1 plants when U HP [t] � 0), to a mutualistic type (i = n), which is the optimum when nearly the whole of the carrying capacity is due to such relationship (e.g., U HP [t] � K P [t]). See model assumptions in Table 5.

End-state condition
A simulation ends when both populations and their respective type distributions are stable; i.e. no further change occurs given current conditions. More specifically, we use the RelTol method to decide if the absolute difference between the populations between time t and t-1 is Table 5. Assumptions on population diversity and coevolution.

Domains Assumptions
On the evolution of traits A population can be divided into types according to one or more traits. The distribution of individuals among types can vary in time, due to factors affecting trait transmission.

On the factors affecting the evolution of traits
Change of the population distribution among types depends on the previous population distribution: the more frequent is a type, the more likely it will be imitated or transmitted to the next generation. Change of the population distribution among types depends on the relative fitness of types: the greater the fitness score associated to a type, the more likely it will be imitated or transmitted to the next generation.
Change of the population distribution among types depends on undirected variation.
On the coevolution of traits related to human-plant mutualism The utility given by an individual varies within types. The utility given by other resources to a population varies within its types. The fitness of human types is modified by the relative weight of plant utility in the carrying capacity for humans The fitness of plant types is modified to the relative weight of human utility in the carrying capacity for plants.
https://doi.org/10.1371/journal.pone.0260904.t005 very small, less than 10 −� where � = 6 in our default setting (see reltol_exponential in Table 1). End-states defined by unchanged variables are known as stationary points. Exceptionally, under certain parameter settings, the HPC model does not converge into a stationary point but enters an oscillatory state. To handle these rare cases and others producing extremely slow-paced dynamics, simulations are interrupted regardless of the conditions after time max iterations (max_iterations, in the implementation in R).

Output variables
The most important output variables are the coevolution coefficients (coevo H [t] or coevo-lution_coefficient_humans, coevo P [t] or coevolution_coefficient_plants), which measure the trend in the distribution of a population among its types Eq (9).
The dependency coefficients (depend H [t] or dependency_coefficient_humans, depend P [t] or dependency_coefficient_plants) express the direction and intensity of the selective pressure caused by the other population. It is calculated as the slope coefficient of a linear model of the fitness scores (fitness Hi [t] or fitness_humans, fitness Pi [t] or fit-ness_plants) using the type indexes (types H or type_indexes_humans, types P or type_indexes_plants) as an independent variable.
Positive values of both these coefficients reflect the tendency of a population towards the most mutualistic types (effective coevolution), while negative values indicate an inclination towards the non-mutualistic type due to a low selective pressure exerted by the mutualistic relationship.
We recorded the time step at the end of each simulation (time end or time_end), obtaining a measure of the overall duration of the process. Whenever applicable, we registered the duration of change towards stronger mutualism types in both populations (timing H or timin-g_humans, timing P or timing_plants). We consider change to be effective when at least half of a population is in the higher quarter of the type spectrum, with the respective coevolution coefficient being greater than 0.5 in a scale from -1 to 1 (coevo θ or coevolution_threshold).

Experimental design
Although relatively simple, the HPC model has a total of 17 parameters. We did not engage in fixing any of these parameters to fit a particular case study as a strategy to reduce the complexity of results. In turn, as our aim is to theoretically explore human-plant coevolution, we scrutinised the 'multiverse' of scenarios that potentially represent the relationship between any given human population and any given plant species. The complexity of the model was managed by exploring the parameter space progressively, observing the multiplicity of cases in single runs, two and four-parameter explorations, and an extensive exploration including 15 parameters (all, except ini H and ini P ). The latter modality of exploration was performed by simulating 10,000 parameter settings sampled with the Latin Hypercube Sampling (LHS) technique [66] and Strauss optimization [67]. All simulation runs were executed for a maximum of 5,000 time steps, but most reached the end condition much earlier.

Model implementation and additional materials
The source files associated with the HPC model are maintained in a dedicated online repository [68]: https://github.com/Andros-Spica/hpcModel. This repository contains several additional materials, including a web application to run simulations and the full report on the sensitivity analysis.
The Human-Plant Coevolution model can generate trajectories with or without the final occurrence of human-plant coevolution. Moreover, simulations revealed a broad spectrum of cases (Fig 3), including those where coevolution produces oscillatory or asymmetric change.
Throughout all conditions explored, the results show that a completely successful coevolutionary trajectory, where both populations effectively change, is relatively demanding and it can be deemed unlikely, considering the entirety of the parameter space explored. Furthermore, in light of these results, plant populations are systematically more sensitive to the selective pressure of mutualism than humans, arguing for the scarcity of cases of origins of agriculture in comparison to a relative abundance of effective domestication processes.

End-states
The wide variety of end-states produced by the HPC model can be classified in three general groups: • Coevolution does not occur. Simulation runs in which a stationary point is reached without successful coevolution, thus returning a stable state where humans and plants have a weak mutualistic relationship.
• Coevolution occurs. Both populations go through successful coevolution and become stable only once they have shifted towards stronger mutualism types.
• Coevolution occurs partially, encompassing two types of end-states: • Stationary suboptimal mutualism: One or both populations undergo a significant, but partial change, remaining relatively well distributed among types, or • Oscillatory coevolution: Both populations become trapped in an endless cycle alternating engagement (strong mutualism) and release (weak mutualism).
Coevolution does not occur. Under some conditions, equilibrium is reached without coevolution taking place and consequently both human and plant populations are kept at relatively low densities (Fig 4). Without coevolution, the plant population exists mainly in the non-anthropic niche (U bP � U HP , or utility_other_to_plants is much greater than utility_humans_to_plants) and in wild forms (pop P 1 � pop P n , or type_propor-tions_plants [1] is much greater than type_proportions_plants[number_-types_plants]), while the bulk of human subsistence comes from other resources and only marginally from gathering these plants (U bH � U PH , or utility_other_to_humans is much greater than utility_plants_to_humans), which most humans do opportunistically and with little impact (fitness P 1 � fitness P n , or fitness_plants [1] is much greater than fitness_plants[number_types_plants]). End-states of this type can still diverge significantly due to different parameter settings.
Coevolution occurs. As intended, the HPC model is able to generate trajectories where equilibrium is reached with coevolution, and mutualism between humans and plants is reinforced (Fig 5; Animation 2). The plant population relies more on the human contribution (U bP � U HP , or utility_other_to_plants much less than utility_humans_to_plants) and humans depend significantly on harvesting these plants (U bH � U PH , or utility_other_to_humans much less than utility_plants_to_humans).
As a general rule, the coevolved human and plant populations reach higher levels compared to their counterparts in non-coevolutionary end-states under similar conditions. The total contribution from one population to the other will increase when coevolution happens, In most cases where coevolution happens, the difference between the pseudo-stable and stable population levels before and after coevolution is fairly clear. These two levels are visible as the first and second plateaus in the double-sigmoid curve (see population plot in Fig 5, top  left). The steep slope that mediates between these two levels follows the change in the The coevolutionary trajectories can be divided into two phases: • Prior to coevolutionary shift: This is a period during which human and plant populations are effectively coevolving. During this phase, population levels approach their first plateau or pseudo-stable state value before coevolution takes effect while the distribution of types change-first slowly, then abruptly-towards the most mutualistic type. It ends when the change in the distribution of types can be considered completed in both populations; we define this moment to be the latest time step between timing H and timing P (in Fig 5, it is timing P , represented by the pink vertical dashed line).
• Following coevolutionary shift: This is a period characterized by the stabilization of the populations around the truly-stable state. During this phase, both populations can be considered "changed" or effectively coevolved, even though they have still not realised the full potential for population growth made possible by coevolution. Although, depending on the specific conditions set by parameters, this phase typically involves a 'boom' for one or both populations.
Under some conditions, coevolutionary trajectories can display a punctual decrease in carrying capacities towards the end of the first phase, during the change from the least to the most mutualistic types. These demographic "bumps" happen in a population when the stronger mutualism type is less capable of exploiting other resources than the least mutualistic type (e.g., if U bH1 > U bHn , then U bH [t] > U bH [t + 1] during coevolution), while the other population has still not grown enough to counterbalance the loss in carrying capacity. In the example given in Fig 5, the plant population is the one suffering this effect, starting at the vicinity of the shift of the human population (vertical dashed cyan lines). In this case, the most mutualistic plant type is far less capable of exploiting non-anthropic resources than the least mutualistic type (U bP1 or utility_other_to_type_1_plants = 100, U bPn or utility_ other_to_type_n_plants = 20) and the utility given by the human population at that point (U HP � 80) lies below the utility obtained from other resources when the least mutualistic types were the vast majority (U bP [t]�100, for t or time from 1 to 200).
Coevolution occurs partially. Simulation experiments revealed cases in which the coevolution towards stronger mutualism occurs only partially. These cases are relatively rare, considering the entirety of the parameter space explored. However, they illustrate the complexity of the interaction of some factors accounted for in the HPC model.
The two types of end-states that fall into this general category, stationary suboptimal mutualism and oscillatory coevolution, are produced under parameter configurations that generally contain strong asymmetries either between the population or between types within the same population. These asymmetries include, for instance, configurations where one population has the most mutualistic types contributing the same amount of utility per capita than the least mutualistic types. In this scenario, the positive feedback between population growth and change in the distribution of types is weakened, but only enough to impede the change in one population; this is the case of the settings shown in Fig 6 ( � U H 1 P or utility_per_capi-ta_type_1_humans_to_plants = 0.5 and � U H n P or utility_per_capita_ty-pe_n_humans_to_plants = 0.5).

Parameter explorations
The extensive exploration of parameters demonstrated that a multiplicity of factors are involved in plant domestication and the origins of agriculture. However, the results also shed light on the relative importance of each of the factors included in the model. We summarise the roles of the parameters of the model as 'facilitators', 'obstructors', and 'scalers' (Table 6). Under most conditions, increasing the values of any facilitator improves the chances of having a successful coevolution process, while greater values for obstructors will diminish it (respectively, positive and negative correlations with coevolution and dependency coefficients). Scalers vary the size of populations (H or humans, and P or plants) at the end-state and the duration of the processes (time end or time_end, timing H or timing_humans, and timing P or timing_plants). Some parameters fit in more than one of the above classes, depending on the setting of the other parameters. The initial populations (ini H or initial_population_humans, ini P or initial_population_plants) remain outside this classification, having virtually no effect on end-states. Within the range of values explored, all parameters but the initial populations and the intrinsic growth rates (r H or intrinsic_growth_humans, r P or intrinsic_-growth_plants) displayed tipping points, i.e. threshold values beyond which the endstates of simulations change drastically (non-linear effect). The exact location of a tipping point in one parameter depends on the values of all others parameters with tipping points, indicating a generally strong interaction between their effects, and hence no single-cause explanation for a given end-state can be accurate. For instance, in light of the trajectories in Figs 4 and 5, it would be correct to say that coevolution occurs in the latter because utili-ty_per_capita_type_1_plants is high enough (i.e. higher than a threshold value between 0 and 0.15). Yet, it is also correct to affirm that it is so because, simultaneously, uti-lity_other_to_type_n_humans is low enough, undirected_variation_plants is high enough, and so on.
Despite their shared explanatory role, parameters vary significantly in importance when predicting the values of the coevolution coefficients at the end-state. We were able to rank the explanatory power of each parameter by fitting Random Forest Regression models where parameters are inputted as predictors in respect to each coevolution coefficient separately (Fig 7).
The same procedure was applied for the dependency coefficients and timings; see section 5.2 in [68]. The assessment of parameter importance for the dependency coefficients displayed a similar pattern, only highlighting those parameters with a direct impact on the carrying capacity of the respective population (greens and blues). While the intrinsic growth rates have the highest impact on the timing of coevolution, all other parameters are scored similarly, having at least some importance for one or both populations. Parameter explorations revealed that timing indicators (timing H or timing_humans, timing P or timing_plants, and t end or time_end) are larger, the closer parameter values are to a tipping point. In those liminal cases, the coevolutionary process can take up to three times longer.
Number of types, undirected variation and intrinsic growth rate. The numbers of types in human and plant populations (n H or number_types_humans, n P or number_ty-pes_plants) facilitate change (i.e. facilitators). However, these two parameters stand out as the least important. Such a result is desirable given that the aspect regulated by these parameters-i.e. the discretionality of population variation-is a necessary artefact of the model and  can only translate to arbitrary classifications when regarding real populations. Ultimately, every individual in a real population could be a single instance of their own type. The overall low importance of these parameters warrants future explorations to treat these as constants, preferably setting them at values much greater than unity (n A � 1). The levels of undirected variation (v H or undirected_variation_humans, v P or undirected_variation_plants) are also facilitators. With higher variation, there are more individuals belonging to stronger mutualism types. Though unfit for the initial conditions, these are the pioneer individuals that may eventually build up the necessary selective pressure on the partner population and trigger coevolution. The positive relationship between undirected variation and the occurrence of coevolution agrees with Fisher's fundamental theorem of natural selection [72,73], according to which higher variance increases the rate of adaptation of a species; which, in this case, leads to stronger mutualism.
Intrinsic growth rates (r H or intrinsic_growth_rate_humans, r P or intrin-sic_growth_rate_plants) are scalers, conditioning how fast populations levels change. Generally, higher intrinsic growth rates return shorter periods of population growth and change of type distribution. However, because they also define how rapid is the feedback cycle regulating the mutualistic selective pressures, they show a mirrored pattern where the intrinsic growth rate of one population has its greatest impact on the timing of change of the other population.
Utility-related parameters. Overall, the most important parameters in the HPC model are those characterising the potential of the mutualistic interaction between humans and plants (Fig 7); i.e. the utility per capita of type n individuals to the other population ( � U P n H or utility_per_capita_type_n_plants_to_humans, � U H n P or utility_per_-capita_type_n_humans_to_plants). Although the correspondent values for type 1 individuals ( � U P 1 H or utility_per_capita_type_1_plants_to_humans, � U H 1 P or utility_per_capita_type_1_humans_to_plants) also play a significant role, coevolution is more often enabled by the utility given by the higher-end types in the mutualistic spectrum. The effect of these parameters is mirrored (greens in Fig 7): utility_per_-capita_type_n_plants_to_humans mostly affects change in the human population, and utility_per_capita_type_n_humans_to_plants does it in the plant population. However, utility_per_capita_type_n_plants_to_humans weights considerably on both humans and plants.
All four parameters related to the utility exchange between humans and plants set a range of utility per capita of each population type that amounts to population totals (e.g., U HP or utility_humans_to_plants). Whenever these totals overcome the totals given by the other resources (e.g., U bP or utility_other_to _plants), the fitness scores will favour stronger mutualism types and trajectories will shift towards a successful coevolution (Fig 8).
The parameters determining the utility given by other resources (U bH1 , U bHn , U bP1 , and U bPn ) are obstructors. Overall, the parameters corresponding to the human population (U bH1 , U bHn ) have a stronger effect than those related to plants (blues in Fig 7). The two parameters regulating the utility of other resources to plants (U bP1 , U bPn ) can also be facilitators depending on the conditions set by other parameters; however, their effect is the weakest of all eight parameters associated with utility (greens and blues in Fig 7).
The parameters associated with utility are also important scalers since they have a direct effect on carrying capacities. The parameters contributing to the carrying capacity for humans ( � U P 1 H , � U P n H , U bH1 , and U bHn ) are able to influence scale more freely because they are not capped by MaxArea. In particular, the utility of other resources to type 1 individuals (U bP1 , U bH1 ) can condition almost entirely the respective carrying capacity-and consequently the population levels-at the end state. These parameters alone can generate trajectories where the human population at the end-state varies from a few to thousands of individuals, without ever incurring in coevolution.
Trajectories with coevolution can be very different (compare Figs 3E to 5) mainly due to the amount of space available for plants (MaxArea) and the conditions regulating the mutual utility between humans and plants ( � U H 1 P , � U H n P , � U P 1 H , and � U P n H ). These are important facilitators, but also have the potential for producing end-states that differ dramatically in the sheer size of the human and plant populations (H, P). For instance, an overall low utility of plant types to humans ( � U P 1 H , � U P n H ) can still produce end-states with coevolution that are indistinguishable in terms of human population size from others without coevolution, where the overall utility of other resources to humans is sufficiently high.
Surprisingly, full-fledged coevolution can still happen when type n individuals contribute less than type 1 individuals (e.g., � U P n H < � U P 1 H ). For instance, when � U P n H ¼ 1:5 and � U P 1 H ¼ 3 in Fig 8. This happens whenever the population total (e.g., U PH ) overcomes the amount given by other resources (e.g., U bH ). This discovery indicates that, at least under the assumptions of this model, the adaptation to mutualism could cause the deterioration of the contribution of individual organisms while still increasing population numbers.

Discussion
Much of the groundwork that helped to understand the evolutionary dynamics of plant domestication comes from archaeology, and more specifically from archaeobotany. Harris [74] theorised the process of domestication as composed of three stages: 1) wild food procurement by hunting and gathering societies; 2) cultivation of wild plants; and the 3) domestication syndrome fixation that established true agriculture of domestic plants. The early plant datasets, mostly coming from the Fertile Crescent, were interpreted as suggesting a 'rapid transition' between these stages due to a strong and direct human selection favouring interesting characters, such as non-brittle spikelets in cereals [75] and suppression of seed dormancy in legumes [76]. However, the richer archaeological record of the last few decades suggests that such transitions could involve a period of pre-domestication cultivation lasting thousands of years [77,78], followed by fixation of the emerging domestic traits, a process that again can happen over thousands of years; see e.g. for cereals [79]. This mechanism, leading to the evolution of domesticated and commensal species, seems to have been a response to the emergence of human-modified environments appearing from the end of the last glaciation [80]. Both the domesticated plants and human populations benefited from this co-evolutionary process, leading to stronger mutualism [53].

Multiple factors, multiple scenarios
The HPC model illustrates the multiplicity of the dynamics that, under its theoretical framework (Tables 4 and 5), can explain ecological and socio-economic shifts, such as the so-called neolithisation. The exploration of the model reinforces the premise that, to explain the domestication of plants and the adoption of agricultural practices, we must integrate the different degrees of complexity of the phenomenon itself, and accept that single-factor explanations do not fit this multiple and heterogeneous reality [1,5]. The great variety of scenarios regarding the characteristics of the crops and of the ecological milieus, as well as the different social, cultural and technological settings in human populations, highlights the complexity of the process and the inevitability of generating case-specific narratives when interpreting the evidence. However, the HPC model goes beyond the replication of multiple single-case idiosyncrasies and contains the formalisation of a general mechanism: the coevolution of humans and plants. This model is able to generate a wide diversity of simulated trajectories and end-states, expressed as aggregated quantitative variables, which we hope can be used in the future to produce explanatory frameworks for specific real-world cases. Therefore, the HPC model is not aimed at reproducing historical processes per se but different possible scenarios of humanplant coevolution, which can be searched in and contrasted with specific lines of evidence.
The model points to several aspects that can explain the emergence of agricultural systems. Some of these aspects, like the utility per capita to the other population, have been part of the archaeological and botanical discourse, albeit not as a formal model [75]. Furthermore, the model shows that a small increase or decrease around a threshold value can produce major changes in the system (tipping point) and that, for coevolution to occur, all parameters showing tipping points must be either beyond or below a particular threshold, which, in turn, depends on the values of all other parameters.
The HPC model also shows that certain differences between human and plant populations can have an important effect on the outcome of human-plant coevolution. The selective pressure of one versus the other may vary significantly among parameter settings, thus producing qualitatively different scenarios.
At one end of the mutualism spectrum, the model can generate scenarios where the subsistence relies heavily on the plant population and the selective pressure is sufficient to drive a substantial change on plant type frequency and population levels, thus leading to some form of agricultural system. At the other end, the model produces outcomes where there is low human-on-plant pressure and humans have many (and preferred) alternative food sources. In such instances, wild plant forms are maintained in the population and low densities are retained. Human subsistence in such cases relies mostly upon other resources, which might still allow for high population densities independently of the plant population; e.g., fishing and complex hunter-gatherers [81,82]. Between these extreme endstate scenarios, the model also simulates other "realities" in which only one population exerts enough selective pressure over the other for it to shift towards stronger mutualism types: societies cultivating plants that, though affected, remain not fully domesticated (cultivation without domestication), or those foraging plant populations that increase their productivity without humans investing more time in them (domestication without cultivation).

Intensification and the coevolutionary dynamics of prehistoric plant management
In most early cases, the adoption of agriculture seems to be the culmination of a long process with deep roots in hunter-gatherer societies [83]. Archaeological literature traditionally considers this process to be fuelled by a series of changes related to food resource diversification [84,85] and, particularly for plants, intensification [86][87][88][89]. Within this context of change, intensive gathering and cultivation have been considered economic practices within a continuum, where some plant species are gathered opportunistically and others systematically exploited. At the beginning of every transition to agriculture, predatory strategies (fishing, hunting, and gathering) were central to human subsistence, while mutualism (plant tending and animal husbandry), if any, were complementary [32].
The theoretical continuum between resource management, domestication, and agriculture assumes that the existence of each forgoer component is paramount for the development of the next "step". However, any one of these phenomena does not inevitably lead to the next [4]. Assuming that in some cases there is an effective transition to agriculture, that means that the focus shifts from a wide range of prey-like resource use to a relatively small number of very successful mutualism partners, among which domesticated plants eventually become the basic source of staple food. In this framework, the coevolution between humans and plants can be defined as a process mediating between weaker and stronger mutualism that can involve many stages, each with a qualitative change in the distribution of types and consecutive boom and stabilisation of both populations.
The HPC model allows identifying various regimes of mutualism between humans and plants. The model, in fact, can give rise to a wide range of scenarios that, from the human point of view, consist of different combinations of wild/domesticated plant food resources and modes of exploitation of such resources, with variable commitments in terms of diet and investment. These strategies can be interpreted as mixed economies, which have been shown to be possible, viable and even resilient socio-economic choices. Within the specialized literature, mixed economies are usually understood as minor or marginal socio-economic systems, defined either as the combination of different strategies of low-level food production [90] or as by-products of a transitory, and thus not stable, stage [91].
These strategies are not necessarily implemented as static combinations, but also as seasonal or periodical activities, shifting from one strategy to another [92,93]. In addition, the pursuing of such strategies might not be a clear and rational decision adopted by specific social agents or groups in charge of the economic activities, but a scenario arising by the aggregation of multiple decision-making processes at the community level, throughout generations.
There is a strong relationship between richness of viable economical options and the specialisation and diversification in subsistence strategies [94][95][96][97]. Specialisation and diversification are hypothesised to have first occurred during the Mesolithic as a mean to intensify the acquisition of resources and they are considered a preamble for the implementation of agricultural practices [98]. Although the concept of intensification could support the continuum concept, there is a strong debate about the reasons and conditions under which intensification takes place in hunter-gatherer societies [99].
With the current work, we aim at showing how the succession of mixed economies are an intrinsic part of the coevolutionary dynamics between human and plants, and shed some light on why can these culminate, in many cases, in the emergence of agriculture.

Insights on the Neolithic Demographic Transition
In archaeological theory, the origins of agriculture is often defined as the birth of a new socioeconomic paradigm involving key changes in human demography and social organization, such as increased hierarchy and division of labour. Among these changes, the most striking is the unprecedented population growth that usually followed the adoption of agriculture, i.e. the Neolithic Demographic Transition [100][101][102].
The HPC model considers the relationship between plant utility and human needs (population pressure) but also the positive effects humans can have on plant growth. The latter involves a delayed improvement of plant utility to humans, through the evolution of traits and sheer population growth, and an increasing human population growth, putting pressure on old and new food resources. Low population pressure, given by either low population density or abundance of food resources, has been argued as a precondition for increasing growth rates in human populations [93]. The demographic increase by the end of the Upper Palaeolithic, as shown by the archaeological record, has been considered a possible cause for a series of intensification processes (such as the intensification of plant gathering or the expansion in coastal populations and an increase in the consumption of coast and marine resources). At the same time, either the intensification of resource exploitation and/or the adoption of agricultural practices (both increasing the productivity per area but also involving labour-intensive, timesensitive activities) might have fostered the abandonment of a series of measures controlling fertility, resulting in a population increase.
A few studies have recently focused on the various demographic booms and busts identified during the Early Neolithic in Europe [16] and which may be interpreted as the possible diverse outcomes of the neolithisation process. While neolithisation intuitively implies a population boom due to the overall increase in food availability, not all instances of shifting to an agricultural economy appear to have been demographically successful. The HPC model suggests a possible explanation for population busts within its formal framework: a momentary decrease in the adaptive fitness of the population and, thus, of the carrying capacity of the environment.
The growth of the human population can have a series of implications. First, a higher demand of the available resources that become manifest in the selective pressure on the plant population or other resources (mixed economy). This may have positively affected the domestication process, by increasing plant bulk productivity, but also produced a series of changes fostering the hunter-gatherer strategy to be less effective when combined with a more invested plant cultivation. When cultivation becomes a priority, there is an expectation for societies or groups within societies to become more sedentary, at least seasonally, so that crops are properly monitored during growth. As a consequence, there would be a reduction in the fitness of the hunter-gatherer strategies. Firstly, because some expertise may be lost, even within a generation, as a considerable part of the labour and efforts for cultural transmission would be focused on cultivation. Secondly, with sedentism (or partial sedentism), the catchment area available for foraging would shrink and quickly be impoverished, having less time to recover and at the same time suffering the effects of expanding cultivation practices. Thirdly, the human population will be pressured to adapt to the needs and schedule of the cultivated plant species and the associated labour bottlenecks, which might be incompatible with the dynamics required for gathering or hunting specific wild resources.

Conclusions
Considering the potential of the modeling results, we would like to underline the bullet and conservative nature of the HPC model. All the diversity observed in terms of both attractors and trajectories was generated by the combination of only two submodels, the Verhulst-Pearl Logistic equation and the Replicator Dynamics, which are straightforward benchmark models in theoretical biology. The sole fact that a relatively bullet model can greatly help to understand complex phenomena, such as the origin of agriculture, argues for the use of formal models, and specifically for simulation approaches, in archaeology.
As other examples in the past [53,55], the HPC model demonstrates that population-level (top-down) theory can still produce useful insights. Strong explanatory frameworks can be achieved without the fine insights of case-wise detail; an approach often resisted by archaeologists, but which is at the same time accepted whenever data is interpreted. In this sense, we consider that formal models are fundamental tools to present, demonstrate and explore any theoretical proposal. The HPC model also offers a solid basis for the design and further development of generative (bottom-up) models [51,52,[103][104][105], and is complementary to approaches focusing on plant domestication syndrome through phenotypic and genetic characterisation [106,107].
According to the HPC model, there are several factors involved in the facilitation or obstruction of the emergence of agricultural systems. Although the model confirms the expectation of attributing several causes to the origin of agriculture, it also further explains how multiple factors could be compatible with asserting causation in a historical sense (i.e., concatenation of events).
In the HPC model, the state of the system connecting humans and the plant species is sensitive to almost the totality of the thirteen parameters. More precisely, this sensitivity is expressed as a rather abrupt shift (tipping point) from a weak to a strong mutualistic state, or vice-versa, depending on the threshold values for each parameter, which are in turn dependent on the current values of every other parameter. Then, according to our model, the emergence of agriculture could be explained by the confluence of all these conditions at specific times and places. However, it seems unlikely that, for the same case of emergence, all these conditions change and cross multiple thresholds simultaneously. Conversely, still within the HPC model, we may envisage scenarios in specific regions at specific moments (i.e. under a specific set of other conditions) where the change in few or even one condition triggered the emergence of agriculture. In this case, certain factors may be considered the cause of the phenomenon in a more deterministic sense.
Beyond the identification of factors that play a role in the human-plant coevolutionary dynamics, the HPC model allows assessing the differences in scale and timing between case trajectories. This capability seems to be especially relevant to understand the many cases of non-industrial agricultural systems documented by archaeology and ethnography. By controlling parameter on a case-by-case basis, further work with the HPC model would yield insight on the reliability of particular hypotheses of how agricultural systems emerged in the past, and help explaining why some origins are more observable than others.