Universality in survivor distributions:
Characterising the winners of competitive dynamics
Abstract
We investigate the survivor distributions of a spatially extended model of competitive dynamics in different geometries. The model consists of a deterministic dynamical system of individual agents at specified nodes, which might or might not survive the predatory dynamics: all stochasticity is brought in by the initial state. Every such initial state leads to a unique and extended pattern of survivors and nonsurvivors, which is known as an attractor of the dynamics. We show that the number of such attractors grows exponentially with system size, so that their exact characterisation is limited to only very small systems. Given this, we construct an analytical approach based on inhomogeneous meanfield theory to calculate survival probabilities for arbitrary networks. This powerful (albeit approximate) approach shows how universality arises in survivor distributions via a key concept – the dynamical fugacity. Remarkably, in the largemass limit, the survival probability of a node becomes independent of network geometry, and assumes a simple form which depends only on its mass and degree.
I Introduction
The fate of agents in any situation where death is a possibility attracts enormous interest: here, death does not have to be literal, but could also refer to bankruptcies in a financial context or oblivion in the context of ideas. This panorama of situations is usually captured by agentbased models including predatorprey models, the best known being the LotkaVolterra model Murray (1989); Berryman (1992); Szabó and Fáth (2007) or indeed, winnertakesall models, which embody the survival of the largest, or the fittest. The model we study here belongs to the second category; the most massive agents grow at the cost of their smaller neighbors, which eventually disappear. Motivated by the physics of interacting black holes in braneworld cosmology Majumdar (2003); Majumdar et al. (2005), its behavior in meanfield and on different lattices was investigated at length in Luck and Mehta (2005), followed by simulations on more complex geometries Thyagu and Mehta (2011); Aich and Mehta (2014). The essence of the model in its original context Majumdar (2003); Majumdar et al. (2005) involved the competition between black holes of different masses, in the presence of a universal dissipative ‘fluid’. The model also turned out to have a ‘richgetricher’ interpretation in the context of economics, where it was related to the survival dynamics of competing traders in a marketplace in the presence of taxation (dissipation) Mehta (2013); Mehta et al. (2005).
One of the most important questions to be asked of such models concerns the distribution of survivors, i.e., those agents who survive the predatory dynamics. In spatially extended models, the presence of multiple interactions makes this a difficult question to answer in full generality. The seemingly universal distributions of survivor patterns put forward in Thyagu and Mehta (2011) motivated us to ask: can one provide a theoretical framework for the appearance of some universal features in survivor distributions? The answer is yes, as we will demonstrate in this paper.
We first introduce a much simpler version of the model investigated earlier Luck and Mehta (2005); Thyagu and Mehta (2011); Mehta (2013); Aich and Mehta (2014), which, however, preserves features such as the exponential multiplicity of attractors (Section II) essential to its complexity. We characterize exactly the attractors reached by the dynamics on chains and rings of increasing sizes (Section III), an exercise which illustrates how a very intrinsic complexity makes such computations rapidly impossible. This leads us to formulate an approximate analytical treatment of the problem on random graphs and networks (Section IV), based on inhomogeneous meanfield theory, which yields rather accurate predictions for the survival probability of a node, given its degree and/or its initial mass.
Ii The model
The degrees of freedom of the present model consist of a timedependent positive mass at each node of a graph. These masses are subject to the following firstorder dynamics:
(1) 
Each node is symmetrically (i.e., nondirectionally) coupled to its neighbors, which are all nodes connected to by a bond; is a positive coupling constant.
The dynamics are deterministic, so that all stochasticity comes from the distribution of initial masses . This model is a much simpler version of the one investigated in Luck and Mehta (2005) which was inspired by blackhole physics Majumdar (2003); Majumdar et al. (2005). In particular the explicit time dependence and initial bigbang singularity of that model are here dispensed with; only a simple nonlinearity, quadratic in the masses, is retained. Despite these simplifications, the present model keeps the most interesting features of its predecessor, such as those to do with its multiplicity of attractors.
The coupling constant can be scaled out by means of a linear rescaling of the masses:
(2) 
The new dynamical variables indeed obey
(3) 
As a matter of fact, the model has a deeper dynamical symmetry. Setting
(4) 
where is a global proper time, so that , i.e.,
(5) 
the dynamical equations (1) can be recast as
(6) 
Remarkably, the dynamics so defined are entirely parameterfree. Quadratic differential systems such as the above have attracted much attention in the mathematical literature, such as in discussions of Hilbert’s 16th problem (see e.g. Jenks (1969); Li (2003)).
Equations (6) can be formally integrated as
(7) 
The amplitudes are therefore decreasing functions of . For each node , either of two things might happen:

Node survives asymptotically. This occurs when the integral in (7) converges in the limit. The amplitude reaches a nonzero limit , so that the mass grows exponentially as
(8) 
Node does not survive asymptotically. This occurs when the integral in (7) diverges in the limit. This divergence is generically linear, so that the amplitude falls off to zero exponentially fast in , while the mass falls off as a double exponential in time .
The dynamics therefore drive the system to a nontrivial attractor, i.e., an extended pattern of survivors and nonsurvivors. This attractor depends on the whole initial mass profile (although it is independent of the overall mass scale). The formula (7) generically implies the following local constraints:

each survivor is isolated (all its neighbors are nonsurvivors),

each nonsurvivor has at least one survivor among its neighbors.
Conversely, every pattern obeying the above constraints is realized as an attractor of the dynamics, for some domain of initial data. This situation is therefore similar to that met in a variety of statisticalmechanical models ranging from glasses to systems with kinetic constraints. Attractors play the role of metastable states, which have been given various names, such as valleys, pure states, quasistates or inherent structures Thouless et al. (1977); Kirkpatrick and Sherrington (1978); Stillinger and Weber (1982); Kirkpatrick and Wolynes (1987); Franz and Virasoro (2000). In all these situations the number of metastable states grows exponentially with system size as
(9) 
where is the configurational entropy or complexity. This quantity is not known exactly in general, except in the onedimensional case where it can be determined by means of a transfermatrix approach (see Appendix A). It is relevant to mention the Edwards ensemble here, which is constructed by assigning a thermodynamical significance to the configurational entropy Edwards (1994). According to the Edwards hypothesis, all the attractors of a given ensemble (e.g. at fixed survivor density) are equally probable. This hypothesis holds generically for meanfield models, while it is weakly violated for finitedimensional systems Franz and Virasoro (2000); Dean and Lefèvre (2001); Berg and Mehta (2001); Prados and Brey (2001); de Smedt et al. (2002); Godrèche and Luck (2005).
Iii Exact results for small systems
In this section we consider the model on small onedimensional graphs, i.e., closed rings and open chains of nodes. In one dimension, (1) and (6) read
(10)  
(11) 
for , with appropriate boundary conditions: periodic (, ) for rings and Dirichlet () for chains.
Equations (11) are very reminiscent of those defining the integrable Volterra chain. The coupling term involves the sum in the present model, whereas it involves the difference in the Volterra system. The appearance of a difference is however essential for integrability Mikhailov et al. (1987); Yamilov (2006). The present model is therefore not integrable, even in one dimension.
Our goal is to characterize the attractor reached by the dynamics on small systems of increasing sizes, as a function of the initial mass profile. This task soon becomes intractable, except on very small systems, due to the intrinsic complexity of the model. The numbers and of these attractors on rings and chains of nodes are given in Table 2 of Appendix A. These numbers grow exponentially fast with , according to (9), with given by (86).
Ring with
The system consists of two nodes connected by two bonds. Both attractors consist of a single survivor. The dynamical equations (11) read
(12) 
The difference is a conserved quantity. If , the attractor is (meaning that only node 1 survives) and its final amplitude is
(13) 
and vice versa. The survivor is always the node with the larger initial mass.
In the borderline case of equal initial masses, the integrals in (7) are marginally (logarithmically) divergent. Both masses saturate to the universal limit , irrespective of their initial value.
Chain with
The two nodes are now connected by a single bond. The dynamical equations are identical to (12), up to an overall factor of two. Here too, the more massive node is the survivor.
Ring with
The attractors again consist of a single survivor. Although there is no obviously conserved quantity, the attractor can be predicted by noticing that
(14) 
The sign of any difference is therefore conserved by the dynamics. In other words, the order of the masses is conserved. In particular the survivor is the node with the largest initial mass.
Chain with
The central node 2 plays a special role, so that the two attractors are and . There are two conserved quantities, and . If , the attractor is and the final amplitudes read
(15) 
If , the attractor is and .
Ring with
The two attractors are the ‘diameters’ and . The alternating sum is a conserved quantity. If , the attractor is and . If , the attractor is and . The attractor is therefore always the diameter with the larger total initial mass. The individual asymptotic amplitudes cannot however be determined in general.
Chain with
The three attractors are , and . The alternating sum is a conserved quantity. This is the first case where the attractor cannot be predicted analytically in general.
Figure 1 shows the attractors reached as a function of and , for fixed . It is clear that can only be reached for , i.e., above the diagonal, whereas can only be reached for , i.e., below the diagonal. The intermediate pattern is observed in a central region near the diagonal. The form of this region can be predicted to some extent. On the horizontal axis, the transition from to takes place for . Similarly, on the vertical axis, the transition from to takes place for . The central region where is the attractor shrinks rapidly with increasing distance from the origin. This can be explained by considering the dynamics on the diagonal, i.e., in the symmetric situation where and . These symmetries are preserved by the reduced dynamics
(16) 
The reduced attractor is , the full attractor is , and vanishes identically. The reduced dynamics have another conserved quantity, . The asymptotic amplitude becomes exponentially small as increases. The width of the central green region is expected to follow the same scaling law, i.e., to become exponentially narrow with distance from the origin, in agreement with our observation.
Ring with
There are five attractors consisting of two survivors, obtained from each other by rotation: , , , and . There is no obviously conserved quantity, and the attractor cannot be predicted analytically in general.
Figure 2 shows the attractors as a function of and , for fixed , and . The five attractors meet at point P . If launched at P, the system is driven to the unique symmetric solution where all masses converge to the universal limit . A linear stability analysis around the latter solution reveals that its stable manifold is threedimensional, in agreement with the observation that its intersection with the plane of the figure is the single point P.
This is the first case which manifests one of the most interesting features of the model, that of ‘winning against the odds’ Thyagu and Mehta (2011); Mehta (2013); Aich and Mehta (2014). Numerical simulations with two structureless distributions of initial masses, uniform (uni) and exponential (exp), yield the following observations:

The probability that the node with largest initial mass is a survivor is 0.849 (uni) or 0.937 (exp);

The probability that the attractor corresponds to the largest initial mass sum among the possible attractors is 0.791 (uni) or 0.891 (exp);

The probability that the node with smallest initial mass is a survivor is 0.018 (uni) or 0.019 (exp).
This small but nonzero probability for the smallest initial mass to survive is the beginning of the complexity associated with the phenomenon of winning against the odds. It happens essentially because more distant nodes can destroy massive intermediaries between themselves and the small masses concerned, thereby letting the latter survive ‘against the odds’. For larger sizes, the problem soon becomes intractable. First, the number of attractors grows exponentially fast with ; second, larger system sizes make for increased interaction ranges for a given node. This makes it increasingly probable to have winners against the odds, making it more and more difficult to predict attractors based only on initial mass distributions.
The above algorithmic complexity goes hand in hand with the violation of the Edwards hypothesis as well as other specific outofequilibrium features of the attractors, including a superexponential spatial decay of various correlations Luck and Mehta (2005). This phenomenon, put forward in zerotemperature dynamics of spin chains de Smedt et al. (2002), is reminiscent of the behaviour of a larger class of fully irreversible models, exemplified by random sequential adsorption (RSA) Evans (1993).
Iv Approximate analytical treatment
Despite the complexity referred to above, we show here that some ‘onebody observables’ can be predicted by an approximate analytical approach. Our techniques are based on the inhomogeneous meanfield theory and rely on the assumption that the statistical properties of a node only depend on its degree PastorSatorras and Vespignani (2001). Such ideas have been successfully applied to a wide class of problems on complex networks (see Barrat et al. (2008); PastorSatorras et al. (2015) for reviews). The thermodynamic limit is implicitly taken; also, the embedding graph is replaced by an uncorrelated random graph whose nodes have probabilities to be connected to neighbors, i.e, to have degree . In this section, we use this framework to evaluate the survival probability of a node, given its degree and/or initial mass.
iv.1 Survival probability of a typical node
We consider first the simplest observable – the survival probability of a typical node, irrespective of initial mass or degree.
From the reduced dynamical equations (6), the initial decay rate of the amplitude of node is seen to be:
(17) 
If node has degree , the above expression is the sum of the initial masses of the neighboring nodes.
From a modelling point of view, this suggests a decimation process in continuous time, where nodes are removed at a rate given by their degree at time . The initial graph is entirely defined by the probabilities for a node to have degree . Its subsequent evolution during the decimation process is characterized by its timedependent degree distribution, i.e., by the fractions of initial nodes which still survive at time and have degree . The latter quantities start from
(18) 
at the beginning of the process () and converge to
(19) 
at the end of the process (). Indeed, as in the original model, survivors are isolated, so that their final degree is zero. The amplitude is the quantity of interest, as it represents the survival probability of a typical node. Our goal is to determine it as a function of the probabilities .
The obey the dynamical equations
(20)  
The first line corresponds to the removal of a node of degree at constant rate , while the second line describes the dynamics of its neighbors. The removal of one neighboring node adds to the fraction at a rate proportional to while it depletes it at a rate proportional to .
The timedependent quantity is the rate at which a random neighbor of a given node is removed at time , which is, consistent with the above, given by the average degree of a random neighbor of the node at time . This rate can be evaluated as follows (see e.g. PastorSatorras et al. (2015)). The probability that a node of degree has a neighbor of degree at time is, for an uncorrelated network, independent of , and given by . (The shift from to is related to the ‘inspection paradox’ in probability theory (see e.g. Feller (1966)).) The average degree of a random neighbor of the node at time is then given by , i.e.,
(21) 
Introducing the generating series
(22) 
we see that obeys the partial differential equation
(23) 
with initial condition
(24) 
Hence it is invariant along the characteristic curves defined by
(25) 
This differential equation can be integrated as
(26) 
with
(27) 
We have therefore
(28) 
For infinitely long times, irrespective of , the parameter which labels the characteristic curves converges to the limit
(29) 
which we call the dynamical fugacity of the model.
We are left with the following simple expression for the survival probability of a typical node (see (19))
(30) 
Furthermore, the fugacity can be shown to be implicitly given by
(31) 
where the accent denotes a derivative.
We list a few quantitative predictions for important graphs/networks below:

regular graph. In the case of a regular graph, where all nodes have the same degree , we have . We obtain
(35) (36) For the above results become
(37) (38) 
Geometric graph. In the case of a geometric degree distribution with parameter , i.e.,
(39) we have and . We obtain
(40) (41) 
Generalized preferential attachment (GPA) network. This is a generalization of the BA network, where the attachment probability to an existing node with degree is proportional to Dorogovtsev et al. (2000); Krapivsky et al. (2000); Godrèche et al. (2009), with the offset representing the initial attractiveness of a node. The degree distribution
(45) has a powerlaw tail with a continuously varying exponent . We have, as expected for a tree, , irrespective of , and
(46) where is the Gauss hypergeometric function.
The GPA network can be simulated efficiently by means of a redirection algorithm Krapivsky and Redner (2001, 2002). Every new node is attached either to a uniformly chosen earlier node with probability , or to the ancestor of the latter node with the redirection probability
(47) For (i.e., ), we have a uniform attachment rule, yielding and , so that
(48) (49) For (i.e., ), the model becomes singular. The converge to , while we still have , formally. In this singular limit we have
(50) Finally, the BA network is recovered for (i.e., ).
Figure 3 provides a summary of the above results. The survival probability of a typical node is plotted against the mean degree , for all the above. Turning to the specific example of the GPA networks, we note that (48) and (49) provide lower bounds for the dynamical fugacity and the survival probability respectively, which both increase to their upper bounds (50) as functions of the redirection probability . This trend is reflected in the blue arrowheads in Figure 3, where the blue square corresponds to the value for the BA network.
In order to test the above for a few simple cases, we have measured the survival probability of a typical node in three different geometries: the 1D chain, the 2D square lattice and the BA network. We solved the dynamical equations (6) numerically, with initial masses drawn either from a uniform (uni) or an exponential (exp) distribution. The measured values of are listed in Table 1, together with our approximate predictions (36), (38) and (44).
Geometry  1D  2D  BA 

0.4360  0.3755  0.6767  
0.3679  0.2500  0.5530 
Our numerical results suggest that the survival probability depends only weakly on the mass distribution, as long as the latter is rather structureless. While our predictions are systematically lower than the numerical observations, they do indeed reproduce the global trends rather well.
It is worth recalling here that our analysis relies on a degreebased meanfield approach. Such approximate techniques are expected to give poor results in one dimension, and to perform much better on more disordered and/or more highly connected structures. A systematic investigation of the accuracy of the present approach, including the comparison of various lattices with identical coordination numbers but different symmetries, and of various random trees or networks with the same degree distribution but different geometrical correlations, would certainly be of great interest. While such extensive numerical investigations are beyond the scope of the present, largely analytical work, we hope they will be taken up in future.
iv.2 Degreeresolved survival probability
The above degreebased meanfield analysis can be extended to quantities with a richer structure. In this section we consider the degreeresolved survival probability of a node whose initial degree is given. The degree distribution of this special node obeys the same dynamical equations (20) as those of a typical node:
(51)  
with the specific initial condition . The above dynamical equations along the same lines as in Section IV.1. The generating series
(52) 
obeys the partial differential equation (23), with initial condition
(53) 
Using again the method of characteristics, we get, instead of (30), the following simple behavior for the degreeresolved survival probability:
(54) 
where the dynamical fugacity is given by (31). The form of this suggests a simple physical interpretation: the fugacity measures the tendency of a given node to ‘escape’ annihilation. More quantitatively, is the price per initial neighbor which a node has to pay in order to survive forever. By averaging the expression (54) over the initial degree distribution , we recover the result (30) for the survival probability of a typical node.
We now introduce the survival scale
(55) 
so that our key result (54) reads
(56) 
This representation makes it clear that the survival scale corresponds to the degree of the most connected survivors.
We use this to probe the survival statistics of highly connected graphs. Here, the survival scale is large as a consequence of a large mean degree ; correspondingly, there is a decay in the survival probability of a typical node, as shown in the two examples below:
Figure 3 shows that these differing trends for the ER and geometric graphs are already evident even for low connectivity . This is because of the interesting circumstance that the behavior of the degree distribution for relatively small degrees () determines both the growth of the survival scale as well as the decay of the survival probability in highly connected graphs.
Assuming this regime is described by a scaling form
(59) 
governed by an exponent , we obtain after some algebra
(60)  
(61) 
with
(62)  
(63) 
The exponents and which enter the power laws (60), (61) are always smaller than one (the value corresponding to ), when we find and . This is e.g. the case for the ER graph (see (57)).
In order to test our key prediction (54), the dependence of the survival probability of a node of the BA network on its degree was computed numerically for an exponential distribution of initial masses. Figure 4 shows a logarithmic plot of against degree . We find excellent qualitative agreement with our predictions of exponential dependence, although the measured slope, corresponding to is larger than our analytical prediction of (see (44)).
iv.3 Massresolved survival probability
Here, we extend our analysis to the massresolved survival probability of a node whose initial mass is given. Since this will only enter through the reduced mass , defined as the dimensionless ratio
(64) 
we denote the massresolved survival probability by .
Along the lines of Section IV.1, we derive the following dynamical equations for the degree distribution of the special node:
(65)  
with initial condition (see (18)). To compute the rate at which a random neighbor of the special node is removed, we argue as follows. This rate is a product of the probability of finding a random neighbor of degree , given by , and the rate of removal of this node. The latter is nothing but its effective degree: while the initial masses of its typical neighbors can be taken to be , the special node has an initial mass of , so that the effective degree of this node is . The average degree of a random neighbor of the special node at time is then clearly given by , i.e.,
(66) 
After some algebra, the massresolved survival probability reduces to:
(67) 
where the massdependent fugacity is
(68) 
This formula can be recast as
(69) 
where is the fugacity given by (31). We have, consistently, for .
The massresolved survival probability is an increasing function of the reduced mass . We compute this explicitly for the classes of random graphs and networks considered in Section IV.1.

ER graph. In this case, we have
(70) with
(71) 
regular graph. In this case, we have
(72) with
(73) in the generic case (), whereas
(74) for .

Geometric graph. In this case, we have
(75) with
(76)
The dependence of the survival probability of a node on its initial mass was computed numerically for the 1D chain, the 2D square lattice and the BA network, with an exponential mass distribution. Figure 5 shows plots of the measured values of against . The dashed curves show the prediction (67), rescaled so as to agree with the numerics for . In the three geometries considered, the analytical prediction reproduces the overall mass dependence of reasonably well. The observed dependence is slightly more pronounced than predicted on the 1D and 2D lattices, while the opposite holds for the BA network.
Last but by no means least, there is a striking manifestation of (super)universality in the regime of large reduced masses, when (68) simplifies to
(77) 
The massresolved survival probability then goes to unity according to the simple universal law
(78) 
This key result in the largemass limit is one of the strongest results in this paper: all details of the structure and embeddings of networks disappear from the survival probability of a node, leaving only a simple dependence on its mass and the mean connectivity .
iv.4 Degree and massresolved survival probability
Finally, our analysis can be generalized to the full degree and massresolved survival probability of a node whose initial degree and reduced mass are given.
The degree distribution of this special node obeys
(79)  
where the rate is given by (66), and with the specific initial condition . After some algebra along the lines of the previous sections, we find:
(80) 
This last result encompasses all the previous ones, including the expression (54) for the degreeresolved survival probability and the expression (67) for the massresolved survival probability .
The detailed numerical evaluation of degree and massresolved data on networks is deferred to future work, although we expect the levels of agreement to to be similar to those obtained in Sections IV.2 and IV.3. Our point of emphasis here is the simplicity of (80), which shows that the dynamical fugacity is absolutely the right variable to highlight an intrinsic universality in this problem. In showing that the result (30) is resolvable into components of degree and mass, it also points the way towards a deeper understanding of the ‘independence’ of these parameters in the survival of a node.
Finally, and no less importantly, the result (80) again manifests (super)universality in the regime of large reduced masses (). As a consequence of (77), the degree and massresolved survival probability goes to unity according to the simple law
(81) 
The beauty of this result (as well as its analogue (78)) lies in the fact that it is exactly what one might expect intuitively; it suggests that the probability () that a node of mass and degree might not survive, is directly proportional to its degree and inversely proportional to its mass. In everyday terms, the lighter the node, and the more wellconnected it is, the more it is likely to disappear. The emergence of such startling, intuitive simplicity in an extremely complex system is a testament to an underlying elegance in this model.
V Discussion
The problem of finding even an approximate analytical solution to a model which contains multiple interactions is very challenging. In the context of predatorprey models, the (meanfield) LotkaVolterra dynamical system and the full Volterra chain are among the rare examples which are integrable. Most other nonlinear dynamical models with competing interactions are not, and are quite simply intractable analytically.
The model inspired by black holes Luck and Mehta (2005), on which this paper is based, shows how competition between local and global interactions can give rise to nontrivial survivor patterns, and to the phenomenon referred to as ‘winning against the odds’. That is, a given mass can win out against more massive competitors in its immediate neighborhood provided that they in turn are ‘consumed’ by evermore distant neighbors. When it was found numerically Thyagu and Mehta (2011); Mehta (2013); Aich and Mehta (2014) that such survivor distributions seemed to exhibit somewhat surprising features of universality, it was natural to ask the question: could one find the reasons for such behaviour, in the sense of characterizing these distributions at least approximately from an analytic point of view? An additional motivation was found in the work of the Barabási group Wang et al. (2013) on citation networks, where the authors put forward a universal scaling form for the ‘survival’ of a paper in terms of its citation history.
The blackhole model as defined in its original cosmological context Majumdar (2003); Majumdar et al. (2005) had an explicit time dependence due to the presence of a ubiquitous ‘fluid’, as well as a threshold below which even isolated particles did not survive. Neither of these attributes was necessary for the behaviour of most interest to us, namely the multiplicity of attractors (which in our case involve survivor distributions), and their nontrivial dependence on the initial mass profile, as a result of multiple interactions. One of our major achievements in this paper has been the construction of a much simpler model (without the unnecessary complications referred to above), which still retains its most interesting features from the point of view of statistical physics.
Once derived and established, this simple model was the basis of our investigations of universality in survivor distributions. The exact characterisation of attractors as a function of the initial data becoming rapidly impossible, we were led to think of approximate analytical techniques. Our choice of the inhomogeneous (or degreebased) meanfield theory was motivated by our emphasis on random graphs and networks in earlier numerical work Thyagu and Mehta (2011); Mehta (2013); Aich and Mehta (2014). This approach was embodied in an effective decimation process. Some of the analytical results so obtained were robustly universal, including the exponential falloff (54) of the survival probability of a node with its degree, or the asymptotic behaviours (78) and (81) in the largemass regime.
Our approach led us to introduce the associated concept of a dynamical fugacity, key to unlocking the reason behind the manifestation of universality in diverse survivor distributions. Physically, this signifies the tendency of a typical node in a network to escape annihilation, which we illustrate via a simple argument. Every time an agent encounters other, potentially predatory, agents it pays a price in terms of its survival probability: as the probability of each such encounter is independent of the others, the ‘cost’ to the total probability is multiplicative in terms of the number of predators encountered. The dynamical fugacity is then nothing but the cost function per encounter (i.e., per neighbor), so that the survival probability of the original agent depends exponentially on the number of its competitors. In a statistical sense, this depends only on the degree distribution of the chosen network, leading to the emergence of a universal survival probability for a given class of networks. The asymptotic survival probability of very heavy nodes becomes ‘superuniversal’, in the sense of losing all dependence on different geometrical embeddings. Its form also has an appealing simplicity, as the complement of the ratio (degree/mass) of a given node; the heavier the node and the more isolated it is, the longer it is likely to survive.
In conclusion, we have used inhomogeneous meanfield theory to formulate and solve analytically an intractable problem with multiple interactions. While our analytical solutions are clearly not exact (due to the technical limitations of meanfield theory) they are nevertheless the only way to date of understanding the behaviour of the exact system. In particular, and importantly, our present analysis strongly reinforces the universality that has indeed been observed in earlier numerical simulations of this problem Thyagu and Mehta (2011); Mehta (2013); Aich and Mehta (2014). That such universal features emerge in a highly complex manybody problem with competing predatory interactions, is nothing short of remarkable.
Acknowledgements.
It is a pleasure to thank Olivier Babelon and Alfred Ramani for interesting discussions.Appendix A Numbers of attractors and complexity in one dimension
In this appendix we investigate the attractor statistics of the onedimensional problem by means of the transfermatrix formalism. We describe an attractor as a sequence of binary variables or spins:
(82) 
From a static viewpoint, attractors are defined as patterns obeying the constraints listed in Section II. They can therefore be identified with sequences which avoid the patterns 11 and 000. The last two symbols of such a sequence may therefore be 00, 01 or 10. The numbers , and of attractors of length of each kind obey the recursion
(83) 
where the transfer matrix reads
(84) 
Its characteristic polynomial is , and so the CayleyHamilton theorem implies the recursion
(85) 
Hence all the numbers grow exponentially with , in agreement with (9). The complexity reads
(86) 
with being the largest eigenvalue of , i.e., the largest root of .
The total numbers of attractors on chains and on rings of nodes obey recursions derived from (85), i.e.,
(87) 
with two different sets of initial conditions. The sequences and are listed in the OEIS OEIS (), respectively as entries A001608 and A000931, together with many combinatorial interpretations and references. The first few terms are listed in Table 2.
0  2  3  2  5  5  7  10  12  17  22  29  
1  2  2  3  4  5  7  9  12  16  21  28 
The transfermatrix approach can be generalized in order to determine the densitydependent static complexity , characterizing the exponential growth of the number of attractors with a fixed density of survivors. Introducing a static fugacity conjugate to the number of survivors, the transfer matrix becomes
(88) 
whose characteristic polynomial is . The reader is referred to de Smedt et al. (2002); Godrèche and Luck (2005) for details. We obtain after some algebra
(89) 
Figure 6 shows a plot of this quantity against survivor density in the allowed range (). The expression (89) reaches a maximum equal to (see (86)) when equals the mean static density of survivors:
(90) 
References
 Murray (1989) J. D. Murray, Mathematical Biology (Springer, Berlin, 1989).
 Berryman (1992) A. A. Berryman, Ecology 73, 1530 (1992).
 Szabó and Fáth (2007) G. Szabó and G. Fáth, Phys. Rep. 446, 97 (2007).
 Majumdar (2003) A. S. Majumdar, Phys. Rev. Lett. 90, 031303 (2003).
 Majumdar et al. (2005) A. S. Majumdar, A. Mehta, and J. M. Luck, Phys. Lett. B 607, 219 (2005).
 Luck and Mehta (2005) J. M. Luck and A. Mehta, Eur. Phys. J. B 44, 79 (2005).
 Thyagu and Mehta (2011) N. N. Thyagu and A. Mehta, Physica A 390, 1458 (2011).
 Aich and Mehta (2014) S. Aich and A. Mehta, Eur. Phys. J. Special Topics 223, 2745 (2014).
 Mehta (2013) A. Mehta, in Econophysics of Systemic Risk and Network Dynamics, edited by F. Abergel, B. K. Chakrabarti, A. Chakraborti, and A. Ghosh (Springer, Milan, 2013), New Economic Windows, pp. 141–156.
 Mehta et al. (2005) A. Mehta, A. S. Majumdar, and J. M. Luck, in Econophysics of Wealth Distributions, edited by A. Chatterjee, B. K. Chakrabarti, and S. Yarlagadda (Springer, Milan, 2005), pp. 199–204.
 Jenks (1969) R. D. Jenks, J. Diff. Equations 5, 497 (1969).
 Li (2003) J. B. Li, Int. J. Bif. Chaos 13, 47 (2003).
 Thouless et al. (1977) D. J. Thouless, P. W. Anderson, and R. G. Palmer, Phil. Mag. 35, 593 (1977).
 Kirkpatrick and Sherrington (1978) S. Kirkpatrick and D. Sherrington, Phys. Rev. B 17, 4384 (1978).
 Stillinger and Weber (1982) F. H. Stillinger and T. A. Weber, Phys. Rev. A 25, 978 (1982).
 Kirkpatrick and Wolynes (1987) T. R. Kirkpatrick and P. G. Wolynes, Phys. Rev. A 35, 3072 (1987).
 Franz and Virasoro (2000) S. Franz and M. A. Virasoro, J. Phys. A 33, 891 (2000).
 Edwards (1994) S. F. Edwards, in Granular Matter: An Interdisciplinary Approach, edited by A. Mehta (Springer, New York, 1994).
 Dean and Lefèvre (2001) D. S. Dean and A. Lefèvre, Phys. Rev. Lett. 86, 5639 (2001).
 Berg and Mehta (2001) J. Berg and A. Mehta, Europhys. Lett. 56, 784 (2001).
 Prados and Brey (2001) A. Prados and J. J. Brey, Phys. Rev. E 66, 041308 (2001).
 de Smedt et al. (2002) G. de Smedt, C. Godrèche, and J. M. Luck, Eur. Phys. J. B 27, 363 (2002).
 Godrèche and Luck (2005) C. Godrèche and J. M. Luck, J. Phys. Cond. Matter 17, S2573 (2005).
 Mikhailov et al. (1987) A. V. Mikhailov, A. B. Shabat, and R. I. Yamilov, Russian Math. Surveys 42, 1 (1987).
 Yamilov (2006) R. Yamilov, J. Phys. A 39, R541 (2006).
 Evans (1993) J. W. Evans, Rev. Mod. Phys. 65, 1281 (1993).
 PastorSatorras and Vespignani (2001) R. PastorSatorras and A. Vespignani, Phys. Rev. Lett. 86, 3200 (2001).
 Barrat et al. (2008) A. Barrat, M. Barthélemy, and A. Vespignani, Dynamical Processes on Complex Networks (Cambridge University Press, Cambridge, 2008).
 PastorSatorras et al. (2015) R. PastorSatorras, C. Castellano, P. van Mieghem, and A. Vespignani, Rev. Mod. Phys. 87, 925 (2015).
 Feller (1966) W. Feller, An Introduction to Probability Theory and its Applications (Wiley, New York, 1966).
 Erdös and Rényi (1959) P. Erdös and A. Rényi, Publ. Mathematicae 6, 290 (1959).
 Erdös and Rényi (1960) P. Erdös and A. Rényi, Publ. Math. Inst. Hungar. Acad. Sci. 5, 17 (1960).
 Barabási and Albert (1999) A. L. Barabási and R. Albert, Science 286, 509 (1999).
 Barabási et al. (1999) A. L. Barabási, R. Albert, and H. Jeong, Physica A 272, 173 (1999).
 Dorogovtsev et al. (2000) S. N. Dorogovtsev, J. Mendes, and A. Samukhin, Phys. Rev. Lett. 85, 4633 (2000).
 Krapivsky et al. (2000) P. L. Krapivsky, S. Redner, and F. Leyvraz, Phys. Rev. Lett. 85, 4629 (2000).
 Godrèche et al. (2009) C. Godrèche, H. Grandclaude, and J. M. Luck, J. Stat. Phys. 137, 1117 (2009).
 Krapivsky and Redner (2001) P. L. Krapivsky and S. Redner, Phys. Rev. E 63, 066123 (2001).
 Krapivsky and Redner (2002) P. L. Krapivsky and S. Redner, J. Phys. A 35, 9517 (2002).
 Wang et al. (2013) D. Wang, C. Song, and A. L. Barabási, Science 342, 127 (2013).
 (41) OEIS, The OnLine Encyclopedia of Integer Sequences, URL http://oeis.org.