Resource Concentration and Clustering in Replicator Dynamics with Stochastic Reset Events

As a model for economic and ecological systems, replicator dynamics represent a basic form of agent competition for finite resources. Here, we investigate the effects of stochastic resetting in this kind of processes. Random reset events abruptly lead individual resources to a small value from which dynamics must start anew. Numerical results show that resource distribution over the population of competing agents develops highly nonuniform profiles, exhibiting clustering and fluctuations with anomalous dependence on the population size. This non-standard statistical behavior jeopardizes an analytical treatment based on mean-field assumptions. We propose alternative simplified analytical approaches which provide a stylized description of entropy evolution for the clustered distribution of resources and explain the unusually slow decrease of fluctuations.


Introduction
In theoretical biology, a replicator is an abstract unit capable of creating copies of itself through interaction with the environment [1,2]. This very generic concept-which provides a unified tool for studying evolutionary dynamics at several levels-encompasses such entities as nucleic-acid molecules (RNA and DNA), genes, cells, and, of course, living organisms. In the theory of cultural evolution, an analogous notion applies to memes, the units of cultural information, thus extending the same theoretical framework to social and economic phenomena [3]. The concept of replicator turned out to be especially fruitful within evolutionary game theory, as a model for biological evolution under natural selection. In this context, replicators represent strategies whose individual profit, measured by their relative reproduction success, depends on both their intrinsic fitness and their mutual interaction [4].
Replicator dynamics is a mathematical model, used in evolutionary game theory, that describes how the relative prevalence of different strategies changes in time [5,6]. If, in a large population, x i (t) is the fraction of players adopting strategy i at time t, replicator dynamics prescribe thatẋ (i = 1, 2, . . . , N), where f i (x) denotes the fitness of strategy i, and generally depends on all the components of x = (x 1 , x 2 , . . . , x N ). It can be seen that the N-dimensional simplex, given by ∑ i x i = 1 with x i ≥ 0 for all i, is invariant under Equation (1), and also acts as a global attractor for all non-negative initial conditions. From the perspective of population dynamics, Equation (1) can be interpreted as the time evolution of N interacting species with fitnesses f i (x), additionally subjected to a global mechanism of growth limitation, given by the second term in the brackets, which asymptotically constrains populations to the subspace where ∑ i x i = 1. In this work, we adopt a similar interpretation, where x i represents the resources (richness) of an economic agent i in a population of N interacting agents.
In the simplest version of replicator dynamics, all fitnesses are constant: f i (x) = λ i for all i [7]. In this situation, the first term in the right-hand side of Equation (1) induces an exponential growth of the resources x i , at rate λ i . The opposing effect of the second term, however, limits this growth. For sufficiently long times, in fact, the system approaches the N-dimensional simplex. The outcome of these contrary trends is that, asymptotically, the replicator with maximal fitness accumulates all the resources. Namely, for t → ∞, Thus, with constant fitnesses, the population always ends in a state where resources are trivially concentrated in just one agent. If two or more agents have identical maximal fitnesses, all the resources become shared between them in proportions depending on the initial values x i (0). Our aim in this paper is to study the effect of reset events on the replicator dynamics with constant fitnesses. Resetting is a stochastic mechanism by which a dynamical variablein the present case, x i (t)-is occasionally brought to a prefixed value, from which its dynamics start anew. This mechanism is able to severely modify the statistical behavior of a dynamical system [8]. In the present case, we expect it to inhibit the accumulation of resources by a single agent or a small group of agents, bringing about a nontrivial resource distribution over the replicator population. To gain insight into the overall behavior of our model, which we present in Section 2, Section 3 is devoted to the numerical and analytical study of the case of a single replicator. In Section 4, we show that the combined effect of replicator dynamics and resetting in a large population with identical fitnesses results in anomalous statistical properties, with an extremely slow decrease of fluctuations as the population size grows. This unusual feature is accompanied by clustering in the amount of individual resources, which, over time, sustains a highly heterogeneous resource distribution over the population. Analytical arguments based on a toy two-cluster model are proposed to explain these numerical observations. Finally, Section 5 is devoted to discussing our main results.

Replicators with Resetting
Stochastic resetting was initially introduced as a mechanism of unbounded growth limitation in the context of demographic dynamics [9,10]. Remarkably, when combined with multiplicative (exponential) growth, it gives rise to long-time power-law distributions for the relevant variables [10,11]. It can therefore be used as a model for the emergence of such distributions in the broad class of phenomena where they are observed [12], ranging from biological taxon abundances [13] to economic resource sharing [14]. Since its introduction more than two decades ago, the statistical effects of stochastic resetting have been studied in a wide variety of dynamical processes, such as transport on networks [15], hydrologic phenomena [16], RNA kinetics [17], and active-particle motion [18], among many others [8].
As described in the Introduction, stochastic resetting acts on a variable x(t), whose evolution is otherwise governed by certain dynamical rules, instantaneously bringing its value to a prefixed level u. Reset events are distributed at random along time, and the evolution of x(t) begins de novo after each resetting. Such events emulate the effect of sudden crises or catastrophic occurrences, where the state of the system under study suffers an abrupt change in a short time [19]. This kind of phenomenon is not uncommon in social and economic contexts [11,20,21].
In replicator dynamics with constant fitnesses λ i , we introduce reset events by proposinġ (i = 1, 2, . . . , N; cf. Equation (1)). Here, represents a Poisson (or shot [22]) noise signal, δ(t) being the Dirac delta function. For each i, the reset times t i,k (k = 1, 2, . . . ) are randomly distributed with uniform frequency q i , so that the average lapse between t i,k and t i,k+1 is q −1 i for all k. The prefactor u i − x i in the last term of Equation (3) insures that each reset event brings x i (t) to the reset value u i . The Markovian stochastic Equation (3) can be dealt with by means of a series of standard methods, notably, the Chapman-Kolmogorov equation, which governs the joint probability distribution of the resources x i (t) [22]. It can also be treated numerically, by a rather intuitive implementation of the Poisson process along discretized time [23]. In the following sections, we use these techniques to study the collective dynamics of the replicator population with resetting.

Dynamics of a Single Replicator with Resetting
As a first step in the analysis of our model, it is instructive to study the case of a single replicator, N = 1. Equation (3) becomeṡ with P(t) = ∑ k δ(t − t k ). The random reset times t k have frequency q. The first term in the right-hand side of Equation (5) makes it clear that, for a single replicator, the deterministic contribution to the dynamics is equivalent to logistic growth [24]. Due to arbitrariness in the choice of time units, the system has two independent parameters only: the ratio q/λ, and the reset value u. Figure 1 shows a pair of realizations of x(t), for u = 0.01 and two values of q/λ, exhibiting qualitatively different behavior. For a relatively small resetting frequency, q/λ = 0.1 (upper panel), x(t) usually has enough time to reach the zone of logistic saturation, just below the level of maximal resources (x = 1). The evolution is only occasionally punctuated by reset events to x = u. On the other hand, when the resetting frequency is larger (q/λ = 2.5, lower panel), x(t) barely transits the zone of exponential growth before it is interrupted by a reset event. In this latter situation, the evolution is very similar to the case where the deterministic part of the dynamics is purely multiplicative, which we have analyzed in detail in a recent contribution [19].
Assuming that the stochastic process represented by Equation (5) reaches a stationary regime for long times, the stationary distribution for x, f st (x), can be obtained from the Chapman-Kolmogorov equation by fixing ∂ t f ≡ 0. In the left-hand side of this equation, the second term represents the probability drift induced by the deterministic logistic dynamics, with v(x) = λx(1 − x). The two terms in the right-hand side are gain and loss contributions originating in reset events. The positive gain term is different from zero only at the reset value x = u, while the negative term represents uniform probability loss at frequency q for all x. On the whole, of course, the two terms compensate each other. For x = u, the effect of the delta-like gain term can interpreted as a boundary condition which connects the solution in the intervals Using this boundary condition, the stationary solution reads for u ≤ x < 1, and f st (x) = 0 otherwise. This time-independent distribution behaves as a power law both for small and large values of x. For q/λ > 1, the exponent of 1 − x is positive, and the distribution has a maximum at x = u while it decays to zero as x → 1. For q/λ < 1, on the other hand, f st (x) exhibits a bimodal profile, with a local maximum at x = u and a divergence at x = 1. This case is illustrated in Figure 2, where we plot the distribution as a function of both x (left panel) and 1 − x (right panel) for u = 0.01 and q/λ = 0.1. The log-log axes emphasize the power-law dependence toward the two ends. Excellent agreement between analytical and numerical results supports the assumption of a well-defined long-time stationary regime for the stochastic process. The bimodal concentration of resources at the extreme values, with the ensuing depletion in the intermediate zone, is a direct consequence of the competing effect of logistic growth, which favors accumulation near the maximum, and of reset events, which populate the zone of lower resources.

Fluctuations and Clustering in Large Homogeneous Populations
Turning now the attention to the case with N > 1, we consider homogeneous replicator populations, in which the parameters u i , λ i , and q i in Equation (3) are the same for all agents. In this situation, agents differ from each other in the individual realizations of the sequence of stochastic reset events only. This homogeneity implies that none of them has an a priori advantage based on fitness, or on the frequency and strength of resetting. Thus, any nontrivial emergent collective behavior should be ascribed to the randomness in the time distribution of reset events.
For a homogeneous population, Equation (3) readṡ with P i (t) given as in Equation (4) with the same resetting frequency q for all i. In turn, stands for the total resources over the population. Assuming that, as in the case of N = 1, the system attains a well-defined stationary state for long times, we expect that x T reaches a constant value if N is large enough. Of course, this requires that resource fluctuations are self-averaging over time and over the ensemble. If these conditions are fulfilled, the stationary distribution for individual resources satisfies Equation (6) The solution is for u ≤ x < 1 and 0 otherwise. The absence of a logistic nonlinearity in Equation (8) determines that f st (x) is now a pure power law; cf. (7). The value of x T in Equation (10) must be obtained self-consistently, requiring that it coincides with the total resources calculated from the distribution f st (x), namely The only positive solution to this self-consistency equation is For a given value of Nu, the total resources vary monotonically from In the first limit, when the resetting frequency is negligible, the population is driven by almost purely replicator dynamics, and one single agent typically concentrates all the resources. When, on the other hand, reset events are dominant, the N agents always have resources close to the minimal value u. The corresponding distributions are In the remaining of this paper, we fix the attention on the case q < λ. Indeed, much as in the case of N = 1 analyzed in Section 3, for q > λ-when reset events dominate over resource growth-the replicator dynamics hardly manifests itself, and evolution does not essentially differ from that of a system of non-interacting multiplicative elements with resetting ( [19], cf. Figure 1b). For brevity, numerical results are presented for just a few parameter sets, which we have found to be representative of more general situations.
Following the same numerical techniques used in the case of a single replicator, we have computed the stationary distribution of individual resources for populations of different sizes, with Nu = 0.01 and q/λ = 0.1. According to the analytical result of Equation (12), all these systems have the same total resources, x T ≈ 0.901. Symbols in Figure 3 show histograms of f st (x) for three values of N, analogous to those presented in Figure 2 for N = 1. Lines stand for the corresponding analytical prediction (10).
It is apparent that, although numerical and analytical results follow the same general trend in the distribution of resources, there are important systematic deviations along the whole interval of the variable x. The deviations decrease in magnitude as the population grows, but are still non-negligible for a large system of 10 5 replicators. For this size and large x, the slopes of the power-law tails in the numerical estimation and the analytical prediction are very similar but, as for the values of the distributions, the former are about one order of magnitude above the latter. The difference has the opposite sign at small x, as shown in the inset. We show in the following paragraphs that these discrepancies originate in the anomalous statistical behavior of the total resources x T (t). Its fluctuations along time, in fact, decay very slowly with the system size N. This indicates that our assumption that x T is constant, used to solve the stationary Chapman-Kolmogorov equation, may only hold for extremely large populations, drastically limiting the usefulness of the analytical approach in this kind of systems.  In all cases, f st (x T ) is sharply peaked around a large value x T ≈ 0.93, and exhibits a broad shoulder for smaller x T . Overall, this behavior is compatible with the analytically predicted value, x T ≈ 0.901, obtained from Equation (12). Note however the rather slow change of the shoulder at small x T as N grows: a variation by a factor of 10 3 in the size of the population leads to a decrease of just above one order of magnitude in the height of the distribution in that zone.

Anomalous Fluctuations of Total Resources
x T This weak dependence on N is remarkably apparent in the coefficient of variation of x T , defined as where is the time average of x T (t), and T is a sufficiently long averaging interval. The coefficient C V encompasses overall statistical properties of f st (x T ) in a single quantity, as a measure of the fluctuations of x T (t) relative to its average. Figure 4b is a log-log plot of C V as a function of N. Across the five orders of magnitude covered by the system sizes, the coefficient of variation only decreases by a factor of 3, and there is no clear indication that it might approach zero as N → ∞. In fact, within this rather wide interval of N, it lacks the typical power-law trend that characterizes the system-size dependence of fluctuations in self-averaging statistical systems (usually, N −z with 0 < z < 1) [25]. This hints at a strongly heterogeneous behavior within the population, and calls for a closer look at the time evolution of individual replicators.

Heterogeneity and Clustering in the Evolution of Resources
The darkest curve in Figure 5a shows the evolution of total resources x T (t) in a population of N = 10 4 replicators, with Nu = 0.01 and q/λ = 0.1. At the initial time, all the replicators have identical resources, x(0) = u. We see that, most of the time, x T (t) fluctuates close to its maximum value. Intermittently, however, total resources exhibit sharp collapses where x T (t) suddenly drops to a small value, followed by a rapid recovery.
Other curves in Figure 5a show x i (t) for the three agents with highest resources at each time. These curves demonstrate the typically heterogeneous resource distribution over the population: most of the time, these three replicators accumulate a large fraction of the total resources. Comparison with x T (t), moreover, illustrates how collapses in total resources usually coincide with a reset event of the richest replicator. As a more compact characterization of heterogeneity in the distribution of resources over the population, we have computed the entropy of the individual shares x i /x T as a function of time: This quantity is depicted in Figure 5b for the same realization as in the upper panel. It shows that, in the intervals between collapses of x T (t), resources progressively accumulate in less and less replicators. Resetting of one of the replicators with high resources, in turn, entails a sudden growth of H(t), with an ensuing decrease as resources become increasingly concentrated. During the intervals between collapses, we expect the population to be divided into at least two groups with different resource distributions inside each group. Those replicators that have undergone a reset event since the latest collapse should have low resources, close to the resetting level u. On the other hand, replicators that have evolved without resource resetting in the same period should possess, on the average, relatively higher resources, with a distribution closer to the equilibrium profile of Equation (10). In a succession of several consecutive collapses, the same mechanism may generate more than two groups, leading to a clustered, markedly heterogeneous resource distribution.
Clustering in the resource distribution is well illustrated by a Zipf plot, in which individual resources are represented against the rank of each replicator in a list sorted by decreasing values of x i . Figure 6 shows snapshots of this kind of plot at four times, in a system of N = 5000 replicators. Other parameters are as in Figure 5. For λt = 89, the first collapse has not taken place yet. In this situation, except for the first-rank replicator which already monopolizes practically all resources, the distribution over the population closely follows the equilibrium profile, whose slope is shown by the dashed line. As time elapses, the occurrence of collapses creates clusters, which in the Zipf plots appear as more or less flat plateaus separated by much sharper steps. In the Supplementary Video S1, which shows an animation of the Zipf plots for the same realization along time, the appearance, evolution, and fading of these plateaus is apparent. Intermittent collapses of total resources and the consequent clustering of resource distribution, leading to an overall highly non-uniform behavior inside the population, are likely determinants of the differences observed between analytical and numerical results, as illustrated by Figure 3, and the slow decay of fluctuations of Figure 4b. In the following, under a few simplifying assumptions, we provide a stylized description for the behavior of the entropy H(t) and a prediction for the typical time between collapses, as well as an argument which explains the extremely slow decay of fluctuations in total resources as the system size grows.

Two-Cluster Model and the Decay of Fluctuations
As a simplified analytical approach to heterogeneity in the replicator population, we propose a toy model in which, at all times between collapses, total resources have the value x T given by Equation (12), and the ensemble is divided into just two clusters. The first cluster contains the N r (t) replicators whose resources have been reset after the latest collapse, occurred at time t c . The second cluster comprises the N − N r (t) remaining replicators. Moreover, we assume that the individual resources in the first cluster are all equal to the reset level u, while the remaining resources are homogeneously distributed over the second cluster. This implies that the total resources in each cluster are N r (t)u and x T − N r (t)u, respectively. With these assumptions, Equation (16) yields where the approximation of the rightmost side holds for u x T . As successive reset events occur, replicators from the cluster of high resources are transferred to the other cluster at rate q so that, on the average, the number of replicators in the former satisfies the equation with N − N r (t c ) = N at the time of the latest collapse. Namely, Replacing into the approximation for the entropy in Equation (17), we find which predicts an approximate linear decay between collapses. The slanted dashed segment in Figure 5b has the slope predicted by this result, displaying very good agreement with the behavior of the numerically obtained signal for H(t).
Our approximation for the entropy H(t) makes it also possible to estimate the typical time between collapses, τ. In fact, in the two-cluster model a collapse will occur when just a single replicator remains in the high-resource cluster, N − N r (t) = 1, accumulating essentially all the resources. In this case, H = 0 which, according to Equation (17), is the entropy attained at time t = t c + q −1 ln N. On the average, the last replicator will be reset after an additional time q −1 . Thus, we have In our simplified picture, τ is nothing but the period of the successive decays of H(t) between its maximum and its minimum. Figure 7a shows the power spectrum P(ν) of an actual numerical calculation of H(t) in a system with N = 1000, Nu = 0.01, and q/λ = 0.1. Its broad profile exposes the stochastic nature of the mechanisms at play in the variation of the entropy, but shows a clear peak at a well-defined frequency, which reveals an underlying time-periodic pattern. The vertical dashed line demonstrates that this frequency coincides quite sharply with the prediction of Equation (21), ν = τ −1 = q/(1 + ln N). We have performed this same comparison for different values of N, evaluating the main period of of numerical signals for the entropy from the position of the highest peak in their power spectra. In Figure 7b, results are compared with Equation (21), represented by the dotted line, with very good agreement. Finally, along the same lines of approximation, we are able to give an explanation for the extremely slow decay of fluctuations in the total resources x T as the system size N grows, revealed by the weak dependence on N of the stationary resource distributions f st (x) and f st (x T ) (Figures 3 and 4a) and explicitly illustrated in Figure 4b. The time signal of x T (t) shown in Figure 5a suggests that fluctuations in total resources are mainly dominated by the collapses associated with resetting of the replicators that accumulate most of the resources. In a highly stylized model for the signal x T (t), we can assume that the statistical distribution of total resources is given by a dichotomic process, where-in the interval between collapses-x T stays at its minimum value Nu during a "recovery time" t R , and at its (approximate) equilibrium value 1 − q/λ during the (average) remaining time τ − t R . Namely, From this Ansatz, the calculation of the mean value and the standard deviation of x T is straightforward. In the limit Nu 1, we find which yields a coefficient of variation If t R is interpreted as the time needed by x T (t) to recover from its small value just after a collapse up to its equilibrium value, we do not expect t R to depend on N, at least for sufficiently large systems. Indeed, according to Equation (8), total resources should approximately obeyẋ T = λx T (1 − x T ) − qx T , which is independent of N. If this is the case, Equations (21) and (24) imply that the coefficient of variation of x T decays as for N → ∞. Symbols in Figure 8 correspond to results for C V as a function of ln N for three different values of q/λ, obtained from numerical solutions of Equation (8) analogous to those of Figure 4b. Dashed lines stand for the asymptotic behavior predicted by Equation (25). Numerical results closely follow the prediction, even for relatively small values of N. On the one hand, Equation (25) shows that C V converges to zero as N grows, which validates the Chapman-Kolmogorov formulation for sufficiently large systems. On the other, the same result proves the extremely slow decay of fluctuations with the population size. Just as an illustration, suppose that one wants to diminish fluctuations in x T by a factor of 10, starting from results for a system of 10 4 replicators. The new system should have nothing less than 10 400 replicators (!), a size clearly beyond the reach of any presently available computational means.

Conclusions
Replicator dynamics with constant fitnesses is a basic model of agent competition, where one or a few agents eventually accumulate all the available resources. In this paper, we have investigated whether this concentration can be mitigated by stochastic resetting in the case of a homogeneous population. Reset events are randomly distributed in time, and force the dynamics of randomly drawn agents to start anew from a small value. Analytical results based on the Chapman-Kolmogorov equation show that, in fact, the long-time distribution of individual resources approaches a smooth profile, with a power-law decay of probability as the amount of resources grows.
However, numerical evidence reveals that-even for long times and large populationsthe analytical prediction is, at most, an approximation to the actually observed resource distributions. A closer inspection of the dynamics of individual agents shows that the overall behavior is still governed by a few agents, which occasionally accumulate most of the total resources. When the resources of one of these wealthier agents are reset, total resources "collapse", and the resource distribution suddenly becomes much more even. Subsequent collapses of this kind lead the distribution to develop clustering, separating the population into groups of agents with similar individual resources. This heterogeneity is responsible for the sustained differences between numerical and analytical results. These collapse-driven dynamics are also responsible for the extremely slow decay of fluctuations with the system size, which jeopardizes the use of the mean-field approach implicit in the Chapman-Kolmogorov Equation (6) for any practically attainable number of agents. Such anomalous statistical behavior is reminiscent of extreme-value statistics, whose relevance to economic processes has been emphasized in various contexts [20,26,27].
The present study complements recent work on cooperative agents subject to stochastic resetting [19], where we have shown that cooperation leads to resource redistribution, distorting the power-law distributions derived from the sole effect of reset events. These contributions represent a first attempt to characterize the collective behavior of interacting agents under the action of resetting, thus combining deterministic dynamics with stochastic ingredients. Other interactions of economic and ecological interest (e.g., parasitism, predator-prey, etc.) are worth considering in future work on the subject.