Cooperation dynamics in the networked geometric Brownian motion

TL;DR

This paper explores how network structures influence cooperation among individuals sharing resources in fluctuating environments modeled by geometric Brownian motion, revealing complex effects on resource distribution and growth.

Contribution

It introduces a model for networked pooling and sharing in fluctuating environments, highlighting the impact of topology on cooperation dynamics and resource distribution.

Findings

Cooperation generally enhances individual growth rates.

Network topology causes significant disparities in resource endowments.

The interplay between cooperation and network structure affects evolutionary outcomes.

Abstract

Recent works suggest that pooling and sharing may constitute a fundamental mechanism for the evolution of cooperation in well-mixed fluctuating environments. The rationale is that, by reducing the amplitude of fluctuations, pooling and sharing increases the steady-state growth rate at which the individuals self-reproduce. However, in reality interactions are seldom realized in a well-mixed structure, and the underlying topology is in general described by a complex network. Motivated by this observation, we investigate the role of the network structure on the cooperative dynamics in fluctuating environments, by developing a model for networked pooling and sharing of resources undergoing environmental fluctuations, represented through geometric Brownian motion. The study reveals that, while in general cooperation increases the individual steady state growth rates (i.e. is evolutionary advantageous), the interplay with the network structure may yield large discrepancies in the observed individual resource endowments. We comment possible biological and social implications and discuss relations to econophysics.

Equations

$$\mathrm{d}{y}_i$$

$$g_i(y_i(t),t)$$

$$\underset{t\to \infty }{\lim }{g}_i(y_i(t),t)$$

$$\underset{N\to \infty }{\lim }{g}_i(⟨y⟩,t)$$

$$\mathrm{d}y$$

$$g_i(y_i(t),t)=\mu -\frac{{\sigma }^2}{2}\frac{1}{N}.$$

$$\mathrm{d}{y}_i$$

$$y_i(t+\Delta t)=\sum _j{A}_{ij}{y}_j(t)\left[1+\mu \Delta t+\sigma {\epsilon }_j(t)\sqrt{\Delta t}\right],$$

$${\bar{y}}_i(t+\Delta t)$$

$$⟨{\bar{y}}_i⟩=\frac{{\sum }_j{A}_{ij}{\bar{y}}_j}{{\sum }_j{A}_{ij}}.$$

$$g_i(y_i(t),t)$$

$$\mathrm{d}f(\mathbf{y},t)$$

$$g(⟨y{⟩}_{\mathcal{N}},t)$$

$$\underset{t\to \infty }{\lim }{\hat{y}}_i(t)$$

$$\underset{t\to \infty }{\lim }{g}_i(y_i(t),t)=\mu -\frac{{\sigma }^2}{2}\frac{⟨v^2⟩}{N}.$$

$$\underset{t\to \infty }{\lim }{g}_i(y_i(t),t)$$

$$=\underset{t\to \infty }{\lim }\frac{1}{t}\log \left(⟨y⟩\cdot z_i\right)$$

$$=\underset{t\to \infty }{\lim }\frac{1}{t}\log \left(⟨y⟩\right)+\underset{t\to \infty }{\lim }\frac{1}{t}\log \left(z_i\right)$$

$$=\underset{t\to \infty }{\lim }\frac{1}{t}\log \left(⟨y{⟩}_{\mathcal{N}}\right)$$

$$\doteq \underset{t\to \infty }{\lim }g(⟨y{⟩}_{\mathcal{N}},t).$$

$$\underset{t\to \infty }{\lim }{\hat{y}}_i(t+\Delta t)$$

$$\approx \underset{t\to \infty }{\lim }\sum _j{A}_{ij}{\hat{y}}_j(t).$$

Full text

††preprint: APS/123-QED

I Introduction

Cooperation has played a fundamental role in the evolution of systems consisting of individuals with different levels of complexity, ranging from simple cell to complex human behavior Axelrod (2006). However, natural selection imposes competition and thus the emergence of cooperation is predicated on the co-occurrence of a specific mechanism within the studied network of contacts Nowak (2006).

A standard approach for examining the effect of different mechanisms on the cooperation dynamics in complex networks is through evolutionary graph theory Lieberman et al. (2005). Under this setting, the individuals interacting in a network are given a set of strategies which they can choose from, and a set of payoffs (changes in the individual resource endowment) that result from interactions with other individuals and their chosen strategies. In the simplest situation, each individual can either be a cooperator or a defector. A cooperator is someone who sacrifices its own resources in order to achieve a better collective performance, whereas defectors are individuals who exploit this cooperative behavior.

Since the pioneering works of Axelrod Axelrod (2006), and later Nowak et al. Nowak and May (1993); Nowak and Roch (2007); Nowak (2006); Allen et al. (2017), on matrix games, i.e. pairwise interactions between individuals, a lot of effort has been put into uncovering the mechanisms required for cooperators to survive the invasion of defectors in networked societies. In particular, more general forms of interaction structures which capture group interactions have been discussed in Perc et al. (2017, 2013); Santos et al. (2008). In this context, it has been found that the introduction of spatial randomness represented by heterogeneous resource endowments between individuals may unconditionally facilitate the evolution of cooperation Kun and Dieckmann (2013); McAvoy and Hauert (2015).

Despite the abundance of studies which capture such spatial stochasticity, a ubiquitous, yet largely unexplored scenario remains the one of cooperative interactions on complex networks in fluctuating environments – where the temporal evolution of resource endowments is strongly affected by their relative growth. In such situations, fluctuations have a net-negative effect on the time-averages, although having no effect on the ensemble (spatial) properties Peters and Klein (2013). This observation, which is a result of the non-ergodicity of the fluctuation-generating process Peters and Klein (2013); Peters and Gell-Mann (2016), yields evolutionary behavior which essentially differs from the one observed in standard models Radicchi and Castellano (2018); Stollmeier and Nagler (2018).

On this basis, it has been hypothesized that repeated pooling and sharing of resources which previously exhibit a fluctuating growth may constitute a fundamental mechanism for the evolution of cooperation in a well-mixed population. The rationale is that, by reducing the amplitude of fluctuations, pooling and sharing increases the steady state growth rate at which the individual cooperating entities self-reproduce Yaari and Solomon (2010); Liebmann et al. (2017); Peters and Adamou (2015). A crucial real-life observation is, however, that interactions between individuals are seldom realized in a well-mixed structure, and they are instead driven by a complex network of contacts Allen et al. (2017).

Motivated by this observation, here we investigate the impact the complex network topology on the cooperative dynamics in fluctuating environments, with networked individuals performing pooling and sharing of resources undergoing a geometric Brownian Motion (GBM). The noisy resource growth produced by GBM is a common model for fluctuations Stollmeier and Nagler (2018); Zheng et al. (2018); Cvijović et al. (2015). The interactions are modeled by considering each individual to also be a pool through which its (direct) neighbors share resources. The model is evaluated analytically and numerically on four types of random graphs: random d-regular graph (RR) Bollobás (2013), Erdos-Renyi Poisson graph (ER) Erdos and Rényi (1960), Watts-Strogatz small-world network (WS) Watts and Strogatz (1998) and Barabasi-Albert scale-free network (BA) Barabási and Albert (1999). Our findings suggest that, while there remains the general trend that cooperation increases the steady state growth rate of each individual (i.e. is evolutionary advantageous), the unique interplay between the non-ergodic fluctuation-generating process and the network topology may generate large discrepancies in the resource endowments. When present, this inequality has a negative effect on the growth rates of the individual entities, hampering their evolutionary performance. Parallels can be made to current societal discussions on wealth inequality Stiglitz et al. (2018).

The remaining of the paper is structured as follows. In Section II, we describe the system model by providing details about the pooling and sharing mechanism, the networked interactions, and the properties of GBM. In Section III we provide analytical results for the growth rate and the steady state behavior of the individual resource endowments. In Section IV we perform numerical experiments and comparison with the analytical results derived in the previous section. Finally, in Section V we discuss our findings and give directions for future work. Some additional technical details are provided in the Appendix.

II Model

II.1 Preliminaries

Formally, we assume that there is a population of non-cooperative individuals, where the dynamics of resources $y_{i}(t)$ of each individual $i$ at time $t$ follow a geometric Brownian motion (GBM),

[TABLE]

with $\mu$ being the drift term, $\sigma$ the noise amplitude, and $\mathrm{d}W_{i}$ is an independent Wiener increment, $W_{i}(t)=\int_{0}^{t}\mathrm{d}W_{i}$. Without noise ($\sigma=0$), the model is simply exponential growth at rate $\mu$. With $\sigma\neq 0$ it can be interpreted as exponential growth with a fluctuating growth rate.

The advantage of modelling through GBM lies in its universality, as it represents an attractor of more complex models that exhibit multiplicative growth Aitchison and Brown (1957); Redner (1990). Its non-ergodicity manifests as the difference between the growth rate observed in an individual trajectory and the ensemble average growth Peters and Klein (2013); Peters and Gell-Mann (2016). In particular, the estimator for the growth rate, $g_{i}(y_{i}(t),t)$, of a single GBM trajectory is defined as

[TABLE]

where $y_{i}(0)$ is the initial condition. For simplicity, we assume $y_{i}(0)=1$.

The time-averaged growth rate is found by letting time remove the stochasticity in the process, i.e., taking the limit as $t\to\infty$, which results in

[TABLE]

The ensemble growth rate, on the other hand, is found by substituting $y_{i}(t)$ with the average $\langle y\rangle$ of an infinite ensemble, where $\langle\cdot\rangle$ is the averaging operation. In other words, one lets the spatial dimension remove the stochasticity by averaging across all possible realizations. Mathematically, the solution is

[TABLE]

where $N$ is the ensemble size.

If only a single system is to be modeled, in steady state only the time-averaged growth rate, Eq. (3), is observed. As discussed in Peters and Klein (2013), the ensemble average growth rate (4) is fictive, as it assumes averaging over “imagined parallel universes”. Hence, in reality, it is the time-averaged growth rate that determines the evolutionary performance of an individual GBM trajectory. Simultaneously, it provides parallels to real-life phenomena. For instance, in evolutionary games the time-averaged growth rate is the geometric mean fitness for the accumulated payoff (resources) of a particular phenotype Sæther and Engen (2015). In economic decision theory, where wealth (resources) dynamics follows a multiplicative process, the same growth observable arises naturally as the unique utility measure Peters and Gell-Mann (2016).

II.2 Pooling and sharing of resources

From an evolutionary perspective, individuals with lower noise amplitude should exhibit higher steady state growth rates and should thus be favored. In this regard, pooling and sharing may constitute a fundamental mechanism for the evolution of cooperation in well-mixed fluctuating environments since it has been found that it reduces the uncertainties in future growth and, hence, brings closer the observed growth rate to the ensemble value Yaari and Solomon (2010); Liebmann et al. (2017); Peters and Adamou (2015). For GBM dynamics, this has been nicely evidenced in Peters and Adamou (2015).

Concretely, the pooling and sharing mechanism can be described as follows. A mutation introduces cooperative dynamics in a population of $N$ individuals whose resource growth is given by a GBM trajectory. In a discretized version of (1), after a period of growth, the individuals pool their resources and subsequently share them equally, resulting in the following dynamics for the resources

[TABLE]

In (5) the subscript $i$ has been dropped due to the equal sharing and $\mathrm{d}W=\frac{1}{\sqrt{N}}\sum_{i}\mathrm{d}W_{i}$ represents the pooled Wiener increment. Evidently, equation (5) is a GBM with an amplitude of $\sigma/\sqrt{N}$, thus yielding a time-averaged growth of

[TABLE]

Notice that as the number of cooperating individuals increases, the time-averaged growth rate converges to the ensemble average growth. This implies that in finite populations, the introduction of new individuals always produces a net performance gain. As a result, one may conjecture that the evolution of group formation and simple multicellularity, where a class of non-cooperating unicellular species mutates to a new trait capable of forming multicellular organisms, could be a consequence of the fact that larger number of cooperators in a fluctuating environment effectively enhances the growth rate (or reduces the drift) Short et al. (2006); Roper et al. (2013). Similar analogy may hold at higher levels of intelligence. As an illustration, consider situations where individuals join a community-supported agriculture to exchange their produced goods for a fixed basket of products, thereby reducing the risks in farming Adam (2006). Another example are nations joining unions to assure sustainable economic growth through common goals Sapir et al. (2004). However, being a model of unconstrained multiplicative growth, GBM has limitations when modeling additive environments or circumstances where growth opportunities are limited due to resource or spatial constraints.

II.3 Networked GBM

Real-life interactions between individuals are, however, seldom realized in a well-mixed structure, and are instead driven by a complex network of contacts Allen et al. (2017). To model this situation, we characterize each individual $i$ with participation in $d_{i}$ pools. In a discretized version of the model, each round $t$ begins with a growth phase where the resources $y_{i}(t)$ of $i$ grow to $\bar{y}_{i}(t+\mathrm{d}t)$. The growth phase is followed by a cooperation phase where each individual pools an equal fraction of its resources in each of the pools it belongs to. Afterwards, each pool returns an equal fraction of the pooled resources to each individual. The resulting mechanism is illustrated in Fig. 1.

The interaction structure is modeled by a connected bipartite random graph $\mathbf{B}$ between finite sets $\mathcal{N}$ of $N$ individuals and $\mathcal{M}$ of $M$ pools, with binary edge variables $B_{im}\in\left\{0,1\right\}$ between pairs of individuals $i\in\mathcal{N}$ and pools $m\in\mathcal{M}$ ($\mathrm{B}_{im}=1$, indicating participation of $i$ in pool $m$). The bipartite representation offers a principled way of capturing wider information regarding the group composition and network interactions Perc et al. (2013). In this regard, the model can be related to games of public goods played on networks, with the main difference that in our model the growth of resources of each individual precedes the pooling phase Perc et al. (2013); Santos et al. (2008); Stojkoski et al. (2018a)111It can be argued that this is more realistic for the examples of cell mutation and agricultural societies given above. In particular, cells first gather nutrients (grow), then share them. Similarly, members of community-supported agricultural first produce their goods then share them in a common pool..

By setting $\mathrm{d}t\to 0$, the dynamics can be explained as

[TABLE]

where $\mathbf{A}$ represents a transition matrix of the network with entries $\mathrm{A}_{ij}=\sum_{m}^{M}\frac{B_{im}}{d_{m}}\frac{B_{jm}}{d_{j}}$ determining the total allocated resources from individual $j$ to individual $i$. Equation (7) resembles the Bouchaud–Mezard wealth reallocation model Bouchaud and Mézard (2000); Garlaschelli and Loffredo (2008); Berman et al. (2017); Ichinomiya (2012), with the note that now the reallocation happens after the growth phase.

III Analytical Results

III.1 Time-averaged growth rate

For tractability, we proceed by examining a discrete version of equation (7),

[TABLE]

where $\varepsilon_{j}(t)$ is a random variable following the standard Gaussian distribution, and utilize a mean-field approach. For this purpose, we define two variables. First, the grown resources of each individual $i$ are given as

[TABLE]

For large $t$ the time-averaged growth rate of this variable should be the same as $g_{i}(y_{i}(t),t)$ as its value will be dominated by $y_{i}(t)$. Second, we define the mean-field around individual $i$ as the average grown resources of each of its neighbors weighted by their contributions to $i$, i.e.,

[TABLE]

By combining the last two equations and adapting the time scale such that $\Delta t=1$, the growth of $i$ can be approximated as

[TABLE]

Two implications arise from equation (9). First, in the transient regime there is an additive term in the growth rate which is solely dependent on the network structure. Hence, during this regime, individuals which are better connected in terms of $\sum_{j}A_{ij}$ should have faster growth rates. The second observation is that the second term on the right-hand side (RHS) of equation (9) eventually converges to the same value for each individual. This is because we study a connected graph where participation in a pool implies that there is a path between any pair of individuals. Due to this interconnectedness, we expect that the steady state time-averaged growth of each $\langle\bar{y}_{i}\rangle$ will be dominated by the growth of the wealthiest individual in the network.

The convergence of the growth rates between individuals provides a direct equivalence with the time-averaged growth rate $g(\langle y\rangle_{\mathcal{N}},t)=\frac{\mathrm{d}\log(\langle y\rangle_{\mathcal{N}})}{\mathrm{d}t}$, which is derived from the partial ensemble average $\langle y\rangle_{\mathcal{N}}$. This object is constructed from all individuals present in the network. As a consequence, one can use Itô’s lemma to directly calculate the time-averaged growth rate in the network. Formally, the lemma states that the differential of an arbitrary one-dimensional function $f(\mathbf{y},t)$ governed by an Itô drift-diffusion process (such as equation (7)), is given by

[TABLE]

In the case of $g(\langle y\rangle_{\mathcal{N}},t)$, we have that $f(t,\mathbf{y})=\log(\langle y\rangle_{\mathcal{N}})$. Then, the first and second derivative of $f$ with respect to $y_{i}$ and $y_{j}$ are easily calculated as $\frac{\partial f}{\partial y_{i}}=\frac{1}{N}\frac{1}{\langle y\rangle_{\mathcal{N}}}$ and $\frac{\partial^{2}f}{\partial y_{i}\partial y_{j}}=-\frac{1}{N^{2}}\frac{1}{\langle y\rangle_{\mathcal{N}}^{2}},$ Peters and Adamou (2018). Moreover, this transformation makes the differential $\mathrm{d}f(\mathbf{y},t)$ ergodic, and since we are looking at steady state averages, $\mathrm{d}y_{i}$ and $\mathrm{d}y_{i}\mathrm{d}y_{j}$ can be substituted with their expected values $\langle\mathrm{d}y_{i}\rangle$ and $\langle\mathrm{d}y_{i}\mathrm{d}y_{j}\rangle$. To estimate these expectations we utilize the independent Wiener increment property $\langle\mathrm{d}W_{i}^{2}\rangle=\mathrm{d}t$, and make use of the fact that $\sum_{k}A_{kj}=1$. Further, we omit terms of order $\mathrm{d}t^{2}$ as they are negligible. As a result, we obtain that $\langle\mathrm{d}y_{i}\rangle=\left[(1+\mu)\sum_{j}A_{ij}y_{j}-y_{i}\right]\mathrm{d}t$ and $\langle\mathrm{d}y_{i}\mathrm{d}y_{j}\rangle=\sigma^{2}\mathrm{d}t\sum_{k}A_{ik}A_{jk}y_{k}^{2}$. By inserting the estimates in equation (10) we can approximate the time-averaged growth rate as

[TABLE]

where $\hat{y}_{i}=y_{i}/\langle y\rangle_{\mathcal{N}}$ are the rescaled resources of individual $i$. This is a dimensionless quantity which compares the endowment of resources of an individual with the population average and as such has been particularly useful in analyses related to wealth inequality Bouchaud and Mézard (2000). In fact, equation (11) indicates that the variance of the rescaled resources $\langle\hat{y}^{2}\rangle_{\mathcal{N}}$ dictates the time-averaged growth rate. Under this model, networks with larger resource inequality, i.e. higher $\langle\hat{y}^{2}\rangle_{\mathcal{N}}$, are expected to have lower steady state growth rates than those where the resources are distributed more equally.

Additional technical details which suggest the usage of the growth rate of the partial ensemble average $g(\langle y\rangle_{\mathcal{N}},t)$ as the growth rate of each individual is provided in the Appendix.

III.2 Steady-state behavior

When deriving the individual growth rate we utilized a steady state property of the system. Such properties are key to understanding the role of complex networks within the pooling and sharing mechanism. In particular, notice that in the limit we can substitute the product of $y_{j}(t)$ and the exponential of (11) for each $\bar{y}_{j}(t+\Delta t)$, divide both sides of the equation by the population average resources and conclude that the steady state rescaled resources of individual $i$ are

[TABLE]

where $v_{i}$ is the $i$-th element of the right-eigenvector of $\mathbf{A}$ associated with the largest eigenvalue normalized in a way such that $\sum_{i}v_{i}=N$. A direct corollary is the equilibrium individual growth rate

[TABLE]

We emphasize that the quantity on the RHS of equation (13) is always greater than $\mu-\sigma^{2}/2$. This can be concluded by examining the optimization problem of maximizing $\langle v^{2}\rangle$ constrained on $\sum_{i}v_{i}=N$, and noting that the global maximum is always less than $N$. Therefore a network of pooling and sharing individuals on the long run will always outperform non-cooperating GBM trajectories. While this indicates that cooperation is a dominant trait in the population, it also asserts that, depending on the distribution of $v$, pooling and sharing may produce societies where the distribution of resources differs to a great extent from the one observed in individual trajectories 222The distribution of resources in non-cooperating GBM trajectories is log-normal.

IV Numerical results

IV.1 Settings

In the numerical analysis we compare the simulated dynamics of the discrete version of the networked GBM, as described with Eq. (8), with the analytical results presented in the previous section. Due to the fact that we can only simulate for finite amount of time and as a consequence may fail to completely remove the stochasticity, we construct partial ensemble averages by averaging the results across 100 realizations of pooling and sharing.

To make the analysis simpler, we formulate the interactions by considering each individual to also be a pool through which its (direct) neighbors share resources. This results in a bipartite graph where the average degree $\langle d\rangle_{\mathcal{N}}$ between individuals is equal to the average degree between pools $\langle d\rangle_{\mathcal{M}}$, i.e. $\langle d\rangle_{\mathcal{N}}=\langle d\rangle_{\mathcal{M}}=\langle d\rangle$. Fig. 2 depicts the process of mapping the original undirected random graph to a directed replacement graph, via a bipartite graph representation which models the pooling and sharing mechanism. The edges in the replacement graph capture the elements $A_{ij}$ in the transition matrix $\mathbf{A}$.

The evaluation of the model properties is done on four types of random graphs: random d-regular graph (RR) Bollobás (2013), Erdos-Renyi Poisson graph (ER) Erdos and Rényi (1960), Watts-Strogatz small-world network (WS) Watts and Strogatz (1998) and Barabasi-Albert scale-free network (BA) Barabási and Albert (1999). To capture the representative graph of each random graph that we study, for each random graph type we average the results across 100 realizations. Moreover, to distinguish the performance of the model in graphs of different size we analyze the model on both small graphs ($N=10$) and large graphs ($N=100$).

IV.2 Experiments

To evaluate the performance of the model under different graph settings we conduct three experiments.

Experiment 1. In the first experiment we examine the transient behavior and the convergence properties of the derived growth rate as described with Eq. (9). The results for small and large networks are respectively given in Fig. 3a and Fig. 3b. Even though we observe that there are some discrepancies at the beginning of the simulation for each graph type an size, eventually the analytical and the numerical growth rate converge to the same value. This result holds for both small and large networks and for each random graph type, thus suggesting the plausibility of our conjecture for the convergence of the growth rate.

Experiment 2. The second experiment compares the distribution of the rescaled resources in steady state, $P_{\mathrm{\hat{y}}}(\hat{y})$, among the graphs. Samples of the corresponding probability density functions (PDFs) are depicted in Fig. 4. Fig. 4a shows the results for small graphs while in Fig. 4b the corresponding results for large random graphs are provided. For both graph sizes, we notice the agreement between the analytical solution in (12) (the value of $v_{i}$) and the simulated rescaled resources, $\hat{y}_{i}$. In addition, independently of the graph size, we observe that the RR graph exhibits no inequality across the resources (point mass PDF), whereas the distributions of the rescaled resources in ER and WS graphs have exponential tails. Finally, the distribution of the rescaled resources in the BA graph resembles a fat tail, i.e. the resources of the individuals exhibit larger variances. As a consequence, the BA graph has the lowest steady state growth rate, followed by ER and WS, as depicted in the inset plots in Fig. 4. This acts as a confirmation for our second analytical finding that steady state growth rate of Eqs. (7) and (8) is uniquely determined by the variance of the right eigenvector associated with the largest eigenvalue of the network transition matrix $\mathbf{A}$.

Experiment 3. The last experiment investigates the role of network sparsity (measured through the average degree $\langle d\rangle$), on the resource distribution. In this respect, it relates the analytical predictions described by Eq. (11) with the numerical solutions of Eqs. (7) and (8). Fig. 5 depicts the variance of rescaled resources $\langle\hat{y}^{2}\rangle$ as a function of $\langle d\rangle$, for small (Fig. 5a) and large graphs (Fig. 5b). Moreover, the inset plots give the ratio of the individual growth rate and the drift parameter, as a function of the same variable. For both graph sizes we observe that denser ER, WS and BA graphs yield more equal resource distribution compared to their respectively sparser counterparts, whereas in the RR graph the resource distribution is invariant to the average degree. As illustrated, there is an alignment between the numerical and the theoretical results for the variance of the rescaled resources both across and within graph types. We note the slight difference between the observed (numerically obtained) growth rate and the analytical solution which, we argue, is due to the fact that the simulations run for a finite amount of steps.

V Discussion

Our findings suggest that interactions on complex networks in a fluctuating environment play a critical role in the observed time-averaged growth rates and resource distribution, both in transient regime and in steady state. The cooperation dynamics is dictated by the properties of the underlying bipartite graph which models the network interactions in the pooling and sharing mechanism.

A startling example is the dynamics taking place on a BA scale-free graph, where largest discrepancies between the individual growth-rates are observed in the transient regime, as compared to ER, WS and RR graphs. Furthermore, the BA graph has the smallest time-averaged growth once the equilibrium is reached, and the most unequal resource distribution. From an evolutionary perspective, a network structure which presents with lower time-averaged growth may be interpreted as being less supportive to cooperation. It is intriguing whether there is any relationship between the apparent lower propensity to cooperation of BA scale-free networks (under the here considered interaction model) and the recent empirical evidence regarding the low-prevalence (i.e. rarity) of scale-free networks in nature Broido and Clauset (2018); Clauset et al. (2009).

As a takeaway, we conclude that inequality may arise as a result of the interwoven relationship between complex networks and cooperative dynamics in fluctuating environments. While it is known that certain network topologies promote inequality Barabási and Albert (1999); Salganik et al. (2006), the effect of cooperative behavior in structured populations is still to be determined Nishi et al. (2015); Chiang (2015); Tsvetkova et al. (2018). As such, our investigations aim at providing deeper understanding on the nature of the relationship between these two occurrences.

Besides providing a basic model of self-reproducing living entities with temporal fluctuations, multiplicative processes are also excessively used to model self-financing investments Peters (2011), gambles Peters and Gell-Mann (2016) and wealth allocation Berman et al. (2017); Bouchaud and Mézard (2000). In this respect, our findings may provide insights to economic utility theory with applications to finance, portfolio management, risk-evaluation and decision-making. In addition, they contribute to the ongoing discussions in economics and econophysics regarding the potential negative effects of wealth inequality on economic growth and development Bouchaud and Mézard (2000); Herzer and Vollmer (2012) and on the individual well-being in general.

A straightforward direction for future work is a scenario where individuals exhibit heterogeneous drifts and volatilities. There, cooperation is evolutionary advantageous only in certain parameter regimes, and thus one should investigate the dynamics under a more general model where individuals are allowed to contribute only a fraction of their resources to the pool. In this context, relations to simplistic behavioral rules that model partial cooperation, e.g. Utkovski et al. (2017); Stojkoski et al. (2018b), may be of particular relevance.

Acknowledgement

This research was supported in part by DFG through grant “Random search processes, Lévy flights, and random walks on complex networks”.

VI Appendix

Here we provide further mathematical logic behind our intuition to use the growth rate of the partial ensemble average $g(\langle y\rangle_{\mathcal{N}},t)$ as the growth rate of each individual. In particular we derive two propositions which describe the dynamics of the system. The first proposition tells us that if the rescaled resources of each individual converge to a certain value, then the growth rate of every individual will also converge to the growth rate of the partial ensemble average. The second one, on the other hand, shows that if the growth rate of each individual converges to the same value, then the rescaled resources converge to $v_{i}$.

Proposition 1: If the rescaled wealth of every individual converges to a certain real value $z_{i}$, i.e. if $\lim_{t\to\infty}\hat{y}_{i}(t)=z_{i}$, with $0<z_{i}<N$, then the steady state growth rate of each individual converges to the growth rate of the partial ensemble average $g(\langle y\rangle_{\mathcal{N}},t)$.

Proof: Suppose that $\lim_{t\to\infty}\hat{y}_{i}(t)=z_{i}$ and that the initial resources $y_{i}(0)=1$ for all $i$, then

[TABLE]

Proposition 2: If the steady state growth rate of each individual converges to the same value, i.e. if $\lim_{t\to\infty}g_{i}(y_{i}(t),t)=g$ for all $i$, then the steady state rescaled resources of individual $i$, $\hat{y}_{i}$ is given by $v_{i}$, where $v_{i}$ is the $i$-th element of the right-eigenvector of $\mathbf{A}$ associated with the largest eigenvalue and normalized in a way such that $\sum_{i}v_{i}=N$.

Proof: Suppose that $\lim_{t\to\infty}g_{i}(y_{i}(t),t)=g$. Then by dividing equation (8) with the average resources $\langle y(t+\Delta t)$ in period $t+\Delta t$, for the discrete version of the model it follows that:

[TABLE]

The last expression gives a Markov chain for which it is widely known that the stationary distribution is given by the right-eigenvector of $\mathbf{A}$ associated with the largest eigenvalue Franzke and Kosko (2011), where its entries $v_{i}$ are normalized in a way such that $\sum_{i}v_{i}=N$.