We use an individual-based model to simulate the dynamics of structured populations (with non-fixed size) consisting of cooperators with varying degrees of cooperation and non-cooperators (free riders). Cooperators produce benefits (group resources that are shared equally by the group members) at a cost to themselves. Individuals have different genetically determined levels of cooperative investment, which determine the group resources that they generate. Group resources are essential for individuals’ reproduction. In other words, more cooperators generate greater group resources, which leads to a higher average reproductive rate of the group members. In contrast, non-cooperators provide no benefit to their groups and bear no cost; they simply consume the group resources.

For a comparison with social populations, we also model non-social populations consisting exclusively of non-cooperators. For simplicity, we consider asexual populations with a mutation rate (m) equal to 0.001. At the beginning of each simulation, population size (N) is set to 300, and all cooperators and free-riders are randomly divided into 90 groups. We assume that the interaction among individuals happens within groups. For completeness, we also model the scenario of non-structured populations and find that cooperation cannot evolve in such a scenario (Supplementary Figure 1), as has been found in other studies (Zhang and Hui, 2011).

Within every simulation time step, cooperators generate group resources, all individuals consume resources and reproduce, and some of them die. After individuals produce offspring, the offspring disperse and randomly join a group. The total number of time steps in a simulation is 10,000, which means that the evolutionary process lasts 10,000 years (roughly several thousand generations), to ensure that the system settles into relatively stable dynamics (Supplementary Figure 2), as suggested in the literature on non-linear dynamics (Strogatz, 2001). We record the proportion of cooperators in populations, the average degree of cooperation that an individual exhibits, the population size, and the total and per capita reproductive output of each group throughout the process. These properties spontaneously emerge from individual-level interactions.

All variables and parameters are summarized in Supplementary Table 1.

The individuals will undergo the following process during each time step:

At the beginning, individuals equally share the group resources (R_{i, t}, where i denotes the i^th group and t denotes the t^th time step). Therefore, when there are N_{i, t} individuals in the i^th group, each individual’s resource consumption (s_{i, t}) is equal to R_i,t/N_i,t. Group resources are determined by environmental resource availability (R₀) and cooperative benefits (bK⁢∑j∅Ki,j). b_K denotes cooperation efficiency, and ∅Ki,j with eleven levels (0.0, 0.1, 0.2, ……, 1.0) denotes the j^th individual’s degree of cooperation in the i^th group. Group resources are a saturating function of cooperative benefits (bK⁢∑j∅Ki,j), which is analogous to the Monod equation (Monod, 1949):

where I is the maximum resource increment rate. This equation means that the more the group members invest in cooperation, the more resources they gain, but there is an upper limit for this beneficial effect. For the groups without cooperators generating benefits, R_{i, t} = R₀.

Next, individuals produce offspring that disperse and randomly join a group. This occurs because in order to simplify the model, we do not assume an explicit spatial structure of the population. However, if offspring are more likely to join their natal patch, this will generally increase the frequency of cooperators within the population (Brauchli et al., 1999; Nadell et al., 2010). The number of offspring the j^th individual can produce (reproductive rate F_{i, j, t}) is a saturating function of the amount of resources it consumed in the form of the Monod equation (Monod, 1949), which also depends on the cost of cooperartion. In other words, an individual’s resource consumption positively influences its reproductive rate with an upper limit α, and if it is a cooperator, it will invest in cooperation at the cost of its own reproduction:

where α denotes the maximum reproductive rate of an individual, β(0 < β ≤ 1) is defined as the percentage decrease in the reproductive rate caused by per unit cooperation degree (∅Ki,j), M is a constant and represents the metabolic consumption of an individual, and K_s is the “half-saturation constant,” which is the value of the individual energy for reproduction (s_{i, t} − M) at which the reproductive rate (F_{i, j, t}) is half of its maximum.

Finally, the system determines whether individuals survive. The survival rate (r_{i, j, t}) of the j^th individual in the i^th group decreases as it gets older:

where c is a constant between 0 and 1. For the offspring born at the t^th time step, r_{i, j, t} = c because age_{i, j, t} = 0.

Next, we introduce environmental fluctuation into the system. Environmental resource availability (R0t) fluctuates periodically, which is described by a sine function:

where A denotes the amplitude and P denotes the period. We record the time series of the proportion of cooperators, the average degree of cooperation, and the population size. To derive the trend in population dynamics, we average the time series from 500 replicate simulations. We further use time-lagged cross-correlation (TLCC) to quantify synchrony between environmental fluctuation and population dynamics at the relatively stable state, and evaluate the variation of the population size by converting the time series of the population size variation standardized by the mean to the frequency spectra using a fast Fourier transformation (Dillon et al., 2016).

We found that the evolution of cooperative behavior is determined jointly by the amount of available environmental resources and the benefits of intraspecific cooperation. Intraspecific cooperation, in terms of both the proportion of cooperators in the population and the average degree of cooperation that each individual performs, is more likely to evolve when environmental conditions are harsh (i.e., low resource availability) and when the benefits of cooperation are large (Figures 2, 3). This occurs because individuals are generally unable to produce offspring without the resources generated by cooperators in harsh environments (Figure 4A). Thus, free riders cannot persist without cooperators under harsh environmental conditions. In addition, although both cooperators and free riders share the group resources generated by cooperators (but only the cooperators have to pay personal costs), cooperators in groups with more group resources can still have more offspring than individuals in groups with fewer cooperators and group resources in harsh environments (Figure 4A). However, resources generated by cooperators play a smaller role in impacting reproduction in benign environments because available environmental resources are already abundant (Figure 4B). As a consequence, cooperation is maintained in harsh environments, particularly when the benefit of cooperating is high, a result consistent with previous models of environmental harshness and intraspecific cooperation (Epstein, 1998; Wang and Goldenfeld, 2011; Zhang and Hui, 2011; Smaldino et al., 2013).

By exploring how environmental quality and intraspecific cooperation jointly influence population dynamics, we found that population size in social organisms is affected by environmental quality both directly in terms of resource availability and indirectly by its effect on the number of cooperators and the degree of intraspecific cooperation within the population. When the benefit of cooperating is small, population size is largely determined by environmental quality, resulting in a population that is similar in size to one without cooperators (Figure 5A). However, as the benefit of cooperating becomes greater, population size is determined by both environmental quality and intraspecific cooperation (Figure 5B). When the benefit of cooperating becomes very large, population size increases abruptly with an increase in environmental quality (i.e., an increase in resources) and then stays relatively constant (Figure 5C). This result can be explained by the fact that the average degree of intraspecific cooperation is also modulated by environmental quality in such a way that individuals are less cooperative in benign than in harsh environments. Therefore, the positive effect of additional resources in benign environments is canceled out by the negative effect of additional free riders. Moreover, additional cooperators in harsh environments compensate for any negative effects of resource scarcity. Furthermore, we showed that social populations possess an advantage in harsh environments with low resource availability because the benefit of cooperating is great enough to outweigh any effects of resource limitation on population size. This result implies that social populations can have wider ecological niches (i.e., can occur in environments with a wider range of resource availability) that non-social populations due solely to the fact that individuals cooperate.

Finally, we compared the population dynamics of social (i.e., cooperative) and non-social (i.e., non-cooperative) populations in a fluctuating environment, finding that the dynamics of non-social populations tend to synchronize with environmental fluctuation (Figures 6A,C), whereas the dynamics of social populations do not (Figures 6B,D). Population size increases after environmental conditions become extremely harsh, but decreases as conditions become more benign. In addition, the peak values in the frequency spectra of population size variation in non-social populations are higher than in social populations (Figures 6E,F), indicating that social populations are more stable than non-social ones in a fluctuating environment because the pattern of intraspecific cooperation is also modulated by environmental conditions (i.e., they are more cooperative in harsher environments), which can buffer the effect of changing environmental conditions on population size (Supplementary Figure 3).

Our results show an unexpected relationship between environmental quality and population size in social organisms. Despite resource scarcity, carrying capacity and the size of a social population can be larger in harsh environments than in benign ones if the benefit of cooperating is high enough. This pattern contradicts the Malthusian view of resource-constrained population dynamics, which does not consider the impact of intraspecific cooperation on population size. It also differs from Boserup’s idea that high population density drives intraspecific cooperation to facilitate population growth, since the degree of intraspecific cooperation is also determined by environmental conditions. In a benign environment, free-riders account for the majority of individuals in a population, and thus there is little benefit of cooperating for further population growth. Therefore, our model synthesizes Malthus’ view of resource-constrained population dynamics (Malthus, 1798) with Boserup’s idea that intraspecific cooperation drives population growth (e.g., via agricultural innovation in humans) (Boserup, 1965). We show that environmental quality influences population dynamics both directly (via resource availability) and indirectly (via the degree intraspecific cooperation within the population).

The magnitude of the benefit of cooperation plays an important role in shaping the relationship between environmental quality and population size. When the benefit of intraspecific cooperation is low, the impact of cooperation on population size is weak, and thus the size of a social population—similar to that of non-social one—is mainly determined by environmental conditions. However, when the benefit of cooperation is high, intraspecific cooperation can strongly influence population size. In addition, the degree to which cooperators invest in creating group resources depends on the environmental conditions, such that they contribute more in harsh environments and less in benign environments. Consequently, the direct relationship between environmental quality and population size is less clear than has been previously assumed (Rodenhouse et al., 1997; Rose, 2000).

Two important implications can be drawn from this observation. First, ecologists often assume that better environmental conditions lead to larger populations (Gotelli, 2008; Molles, 2016). Yet, our finding highlights the need to empirically test such assumptions in social species by quantifying lateral forces—the degree of intraspecific cooperation—to understand their impact on population size. Many theoretical and empirical studies have found that social organisms are often more cooperative in harsh environments (Andras et al., 2007; Connelly et al., 2017), yet it remains unknown whether cooperative behavior allows these species to have higher fitness and generate large populations in harsh vs. benign environments. A rare empirical example comes from a study of cooperatively breeding birds, which found that individuals engaged in less social conflict in harsh environments, resulting in higher group productivity (Shen et al., 2012), thus providing some preliminary support for the role of lateral forces on population fitness. Second, the booming human population sizes of the past century are often considered as evidence that human populations are not limited by resources (Steinmann et al., 1998; Kögel and Prskawetz, 2001). However, we caution against such a view (Boserup, 1965), since our model suggests that the effect of resource availability on the degree of cooperation within a population can also influence population dynamics and can constrain population sizes.

Empirically testing the direct and indirect relationships among environmental quality, cooperation, and population dynamics is urgently needed in social species, including our own. Studies examining the Allee effect and how population density of social species influences individual, group, and population fitness have provided some empirical tests of these relationships (Stephens and Sutherland, 1999; Stephens et al., 1999; Angulo et al., 2018). However, the mechanisms (e.g., cooperative and competitive strategies) by which population density affects population dynamics of social species are often unclear (Angulo et al., 2018). Since our model explores the effects of cooperative behavior on population dynamics under different environmental conditions, it can potentially be used to understand the population dynamics of social species that show Allee effects (Angulo et al., 2018), including meerkats (Paniw et al., 2019), wild dogs (Angulo et al., 2013), and Arabian babblers (Keynan and Ridley, 2016). Similarly, the same rule can be applied to interspecific relationships. For example, the stress gradient hypothesis argues that a harsh environment promotes mutualistic relationships between species (Bertness and Callaway, 1994; Callaway et al., 2002). However, how such a mutualistic and competitive relationship over an environmental gradient affects population size is still unclear and deserves further study (Holland et al., 2002; He et al., 2013).

If we consider resource availability as a dimension of niche space, we can deduce that social species should have a wider niche breadth than non-social species. If the benefit of cooperating is large, intraspecific cooperation helps social species to maintain positive population growth even when the environments are harsh and have scarce resources. This result is consistent with the social conquest hypothesis (Wilson, 1987, 2012), which argues that ants and humans are the two most dominant species on earth because of their ability to cooperate and form groups (Wilson, 1987, 2012). Despite few formal analyses of this idea, a study in Asian burying beetles showed that social groups are more cooperative in harsher, hotter environments with more intense interspecific competition than in benign, cooler environments where interspecific competitors are absent but intraspecific competition is intense (Sun et al., 2014). As a consequence, social groups are better able than the solitary individuals to expand their thermal niche into harsher environments. Our model supports this idea by predicting that social species should have a wider niche breadth than non-social ones.

We found that in fluctuating environments, the size of social populations exhibiting intraspecific cooperation is more stable than that of non-social populations that do not cooperate because environment-associated cooperation buffers the impact of environmental fluctuation on social populations. In other words, the key mechanism leading to population resilience of social species is that they are more cooperative in harsh than benign environments. A study of social microbes also found that social populations are more resilient to environmental disturbance (i.e., experimentally lower population density) than non-social ones because more cooperators generate greater group resources (public goods) (Sanchez and Gore, 2013). However, this study assumed that environmental resource availability was stable and the only change in the experiment was lowering population density. Thus, population resilience in this study meant that populations could survive through periods of low density but not environmental harshness per se, and that the population dynamics were driven by density-dependent processes such that more cooperators in the population generated greater public goods. Several theoretical studies have investigated this density-dependent process (Epstein, 1998; Hauert et al., 2006; Zhang and Hui, 2011), showing similar patterns of dynamics: abundant public goods generated by cooperators favor the rise of free-riders who do not invest in producing public goods, eventually resulting in an overall reduction of public goods. As public goods become scarce, cooperators are favored by selection again. However, these models exploring the effects of cooperation on population dynamics often assume that environmental conditions are fixed, so that the payoff to cooperators and free-riders is largely determined by the strategies of other individuals. Therefore, such models often find oscillating patterns of evolutionary dynamics of cooperators and free-riders. That is, when there are more cooperators, the free-riders are more likely to encounter the cooperator and get a higher payoff, and thus selection will favor the emergence of free-riders. In reality, however, the benefits gained by the free-riders encountering cooperators should still be influenced by environmental resources. Our results suggest that if reproduction is influenced by resource status, then the number of cooperators and free riders may not show oscillating patterns in the population. Since environmental conditions such as resource availability are assumed constant, these studies cannot determine how bottom-up forces influence the interaction between population and evolutionary dynamics. In contrast, our model suggests that environmental quality can influence both the evolution of intraspecific cooperation and population dynamics, a result that should be incorporated in future empirical studies.

As the earth continues to warm, its climate is becoming increasingly unpredictable (Easterling et al., 2000; Rohr and Raffel, 2010). Some studies have argued that climate change-driven resource scarcity will lead to increased armed conflict in human societies as a result of resource scarcity, a truly neo-Malthusian view (Nordås and Gleditsch, 2007; Scheffran and Battaglini, 2011; Raleigh and Kniveton, 2012). Yet, other studies have argued that the environmental problems caused by climate change will not exacerbate violent conflict (Benjaminsen et al., 2012; Gleditsch, 2012), and may even promote peace and greater cooperation (Slettebak, 2012). Based on our models exploring environmental quality, social interactions, and population dynamics, we predict that harsh environments—those with low resource availability—will also promote more cooperation in human societies depending on the types of cooperative benefits that can be generated under different environmental conditions. Empirical studies testing our model predictions by comparing patterns of intraspecific cooperation under different environmental scenarios, as well as those examining population fluctuation and stability between social and non-social species, will be of great importance for understanding the future dynamics of social species—including our own—in a period of continued global change.