Yiping Ma

Sebastian Gonçalves

Sylvain Mignot

Jean-Pierre Nadal

Mirta B Gordon
email: mirta.gordon@imag.fr

Cycles of cooperation and free-riding in social systems

Keywords: Collective Systems, Complex Adaptive Systems, Social Interactions, Cycles

Basic evidences on non-profit making and other forms of benevolent-based organizations reveal a rough partition of members between some pure consumers of the public good (free-riders) and benevolent individuals (cooperators). We study the relationship between the community size and the level of cooperation in a simple model where the utility of joining the community is proportional to its size. We assume an idiosyncratic willingness to join the community ; cooperation bears a fixed cost while free-riding bears a (moral) idiosyncratic cost proportional to the fraction of cooperators. We show that the system presents two types of equilibria: fixed points (Nash equilibria) with a mixture of cooperators and free-riders and cycles where the size of the community, as well as the proportion of cooperators and free-riders, vary periodically.

Introduction

The aim of this paper is to study the relationship between the growth of a group -or community -and the quality of cooperation within this group. This concerns cooperative organizations, usually studied by anthropologists, volunteer organizations like non-profit associations and charities, informal 1 groups in the internet such as knowledge communities, open-source software communities, and many others.

Generally, contributing to the aim of the community bears some cost. But free-riding on the cooperators' efforts may also be a burden because such behavior is contrary to the social norms generally accepted in such communities and may be felt as a failure of personal morality. Thus, even if membership is voluntary, the members of such associations face a social dilemma, because they may benefit of their participation without contributing to the social good (or the common cause). Actually, a polymorphic configuration where cooperators and free riders coexist seems to be a stable form of organization. Clearly, if all the individuals shirk, everybody gets worse and the community may even disappear.

There is a vast literature on social dilemmas. Olson [START_REF] Olson | The Logic of Collective Action : Public Goods and the Theory of Groups[END_REF] and early work on the Prisoner Dilemma, a game theoretical formulation of this type of problems, predicted that selfish individuals have no interest to cooperate in the production of a public good, and as a result, at equilibrium all the players should be free-riders. However, many experimental economic studies (see [START_REF] Ledyard | Public goods: A survey of experimental research[END_REF] for an extensive survey) do not support this claim. Although a stable cooperation is rarely attained in finitely repeated public good games, there is substantial cooperation on the average, at least in the first periods. More generally, reciprocity, reputation-based altruism and sanctions are mechanisms that may help high and stable cooperation levels to emerge among selfish individuals. They give raise to social norms, which are standards of behaviors based on shared beliefs of how individuals ought to behave. Individuals may value approval from their peers to a generous contribution from their part. They may also value negatively their own free-riding behavior when peers contribute. All these mechanisms have proved efficient in increasing the cooperation level in experiments. Altruistic enforcement may be provided by the players after each period [START_REF] Fehr | Altruistic punishment in humans[END_REF] or by external (third party) players [START_REF] Fehr | Social norms and human cooperation[END_REF]. In finitely repeated games [START_REF] Gaechter | Collective action as a social exchange[END_REF], they generally contribute to increase the cooperation rates but do not succeed to achieve complete cooperation. Although these experiments do not explain why some individuals are willing to bear the cost of punishing, they show that such individuals exist.

The experimental results are usually presented as averages of the individuals' behaviors over several groups playing under similar experimental conditions. However, this average may hide quite complex individual behaviors. In a public goods experiment Hichri and Kirman [START_REF] Hichri | The emergence of coordination in public good games[END_REF] showed that, despite very smooth trends of average contributions, individual groups behaviors may have a very large volatility, with oscillations that become completely smeared-out after averaging. Thus, the actual behavior of individuals in public goods games may be much more complex than usually described. Another evidence is the indirect reciprocity game [START_REF] Greiner | Indirect reciprocity in cyclical networks. an experimental study[END_REF], where there is an asymmetry in the social connexions built-in by the oriented ring of in-teractions. Although the aggregate behavior seems almost monotonic along the periods (besides the break at which the protocol of the game -the decision dynamics-changes from parallel to sequential), the individual behavior (the amount invested at each decision conditional to the received amount) is indeed highly oscillatory. However, little attention has been paid to these oscillatory phenomena in the literature of public goods so far.

Theoretically, besides some simple games, which present cycles clearly related to their payoffs structures (like Rock, Paper and Scissors and similar games [START_REF] Gintis | Game Theory Evolving[END_REF]) there are few examples of oscillatory behavior in large systems in the literature. Huberman and Glance [START_REF] Huberman | Beliefs and cooperation[END_REF] have shown that in a public goods problem with large populations, cycles and even chaos may be the asymptotic behavior of the system when individuals adopt particular decision strategies. In particular, if there are opportunistic individuals having a non-monotonic probability of cooperation (cooperate when there are few cooperators but free-ride if the fraction of cooperators increases above a threshold) cycles in the fraction of cooperators may exist. More recently [START_REF] Imhof | Evolutionary cycles of cooperation and defection[END_REF] it has been shown that finite populations may oscillate between three types of strategies when playing the prisonner's dilemma. These are cycles in the strategy space, but do not necessarily give raise to cycles in the decisions or actual payoffs of the players. However, such cycles disappear in the limit of very large systems. In fact, large systems of interacting agents with the usual best response dynamics may present a cyclic behavior at equilibrium whenever the interactions are non-symmetric. Iori and Koulovassilopoulos [START_REF] Iori | Patterns of consumption in a dischoice model with asymmetric interactions[END_REF] have shown that a simple model of consumers with social interactions similar to the one proposed by Durlauf [START_REF] Durlauf | Statistical mechanics approaches to socioeconomic behavior[END_REF] present oscillations in the consumption level provided that the interactions are sufficiently asymmetric. This generalizes well known results in the statistical physics literature, where it is known that the attractors of systems with completely asymmetrical interactions (i.e. the influence of individual i on individual j and that of j on i have opposite signs) may be cycles of order 2 or 4, depending on the decisions dynamics. Moreover, even with symmetric interactions, if all the individuals make their decisions at the same time (parallel dynamics) cycles of order two cannot be excluded (although they are seldom observed in computer simulations).

In this paper we raise the question of the types of equilibria that may appear in the case of heterogeneous populations, when the heterogeneity is related to both the willingness to join the community, and the strength of the moral burden felt by free-riders. Besides the possible fixed point Nash equilibria, we also look for the possible existence of oscillatory behavior, i.e. cycles. More specifically, intuition suggests that a community is more likely to be stable if the amount of cooperation within the group is high. However, the conditions under which a polymorphic community with both free-riders and cooperators may exist in equilibrium, is still an open question. In particular, is there some threshold of cooperation level beyond which such a community may be stabilized? Another question is whether, from a collective point of view, it is better to have a large community with a large proportion of free-riders or a relatively small one but with a majority of cooperators.

We consider the model proposed in [START_REF] Gordon | Cooperation and free-riding with moral cost[END_REF], which is an extension of the basic economic model of binary choices with externality [START_REF] Durlauf | Statistical mechanics approaches to socioeconomic behavior[END_REF][START_REF] Durlauf | How can statistical mechanics contribute to social science[END_REF][START_REF] Nadal | Multiple equilibria in a monopoly market with heterogeneous agents and externalities[END_REF][START_REF] Gordon | Discrete choices under social influence: Generic properties[END_REF]. It considers a population of interacting individuals with idiosyncratic preferences that have to decide whether to join a peers organization and cooperate, whether to join it without cooperating (free-riding) or whether not to join it at all. Cooperators bear a fixed cost, but free-riding has an idiosyncratic cost due to the social disapproval [START_REF] Gaechter | Collective action as a social exchange[END_REF] experienced by a free-rider when facing cooperators.

The paper is organized as follows: section 2 summarizes the details of the model. Section 3 is devoted to the study of the different possible equilibria in the limit of an infinite population under parallel and random sequential dynamics. In section 4 the results are applied to a particularly interesting case. We show that beyond the usual Nash equilibria, that correspond to fixed points, cycles are also expected. The theoretical analysis is compared with computer simulations, showing that cycles persist in finite size systems. The results are discussed in section 6, and we conclude in section 7.

The model of choice with social interactions

We consider a population of N agents in which each individual i must choose one among the following three possibilities:

• to join the community and cooperate (s i = 1)

• to join the community and free-ride (s i = -1)

• not to join the community (s i = 0)

In the following, we denote η c and η f the population fraction of cooperators and of free-riders respectively,

η c = 1 N N k=1 δ s k ,1 (1a)
η f = 1 N N k=1 δ s k ,-1 , (1b)
where δ s,s ′ denotes the Kronecker delta (= 1 if s = s ′ , 0 otherwise). The total fraction of individuals in the population that belong to the organization is η c + η f . We assume that each agent i has a private estimate H i of the value of the community, that determines his idiosyncratic willingness to join (IWJ) it.

Besides, the community value depends on its social composition, increasing proportionally to the fraction of individuals that join it. We assume that this social component takes the simple form of a badwagon effect, given by (J + G)η c + Jη f with J ≥ 0 and G ≥ 0. The term (J + G)η c means that the social benefit produced by each cooperator is (J + G)/N , while it is only J/N per free-rider. With these hypothesis, free-riders also help making the community attractive, although to a lesser extent than cooperators. Therefore, the value of the community for an individual i is

H i + (J + G)η c + Jη f .
Cooperators bear a fixed cost C ≥ 0 that we assume constant, the same for everyone. Free-riders in contrast support a cost proportional to the number of cooperators, weighted by an idiosyncratic factor X i ≥ 0. This cost may be interpreted as a moral burden: X i reflects the importance given by i to the disapproval of cooperating peers.

The utility or surplus U i of individual i depends on his own choice as well as on the choices of the others:

U i (s i = +1|η c , η f ) = H i + (J + G)η c + Jη f -C, (2a)
U i (s i = -1|η c , η f ) = H i + (J + G)η c + Jη f -X i η c , (2b)
U i (s i = 0|η c , η f ) = 0. (2c)
We assume that individuals have a myopic behavior, i.e. at each time step every agent makes the choice which maximizes his utility, estimated from the current observed values of the fractions of cooperators and free-riders. This model is different from the usual games considered in experimental economics of public goods. In the latter, an amount proportional to the individual contributions (Cη in our model) is equally distributed among the members. Such models are suited for explaining cooperation in small groups whose members are remunerated according to the collective output. Here we consider large organizations whose value is strongly related to their collective action, but whose production is not redistributed among the members. This is the case of non-profit organizations and other kinds of communities cited in the introduction. The value of the community is not related to the costs bore by its members, but to their number: larger communities are more attractive.

It is convenient to write the surplus of i as follows:

U i (s i |s -i ) = (A i + B i ) δ s i ,1 + A i δ s i ,-1 (3)
where s -i represents the choices of the other agents. A i is the surplus of joining the community being free-rider, and B i is the bonus for cooperating:

A i = H i + (J + G -X i )η c + Jη f , (4)
B i = X i η c -C. (5)
Then, the best response of agent i given the choices of the other agents may be written as follows:

s i = 1 ⇐⇒ A i + B i > 0 and B i > 0; ( 6a
)
s i = -1 ⇐⇒ A i > 0 and B i < 0; ( 6b
)
s i = 0 otherwise (6c)
Notice that, by eq. (6a) we may have A i < 0 compensated by a positive B i , making thus profitable for i to join the community provided that he cooperates. This may be the case for highly moral individuals, i.e. having large weights X i .

One can note that in this model the social interactions are not symmetric because of the idiosyncratic weights X i . More precisely, the interaction between cooperators is symmetric: if i and j are both cooperators, they contribute equally to their respective utilities, by (J + G)/N . But whenever one free-rider at least is involved, interactions are non-symmetric. For example, if i cooperates and j free-rides, the contribution of i to j's utility is (J + G -X j )/N while the contribution of j to i's utility is J/N . Due to this dissymmetry there is no general result ensuring that the dynamics converges to fixed point solutions. These are guaranteed only for symmetric interactions (see [START_REF] Iori | Patterns of consumption in a dischoice model with asymmetric interactions[END_REF] for a discussion in a case of social systems). Thus, cyclic attractors and even chaos may exist. Notice also that some agents may have X i < J + G, and others X i > J + G: that is, there may be a mixture of positive and negative influences. In this case, reminiscent of spin-glass systems in physics, there may be a very large number of equilibria for some range of parameters. Although multiple equilibria do appear as the generic situation in our model, we will see that their number remains small, at least for the particular case of global interactions considered here.

In our analysis we assume that the idiosyncratic terms H i and X i are randomly distributed among the population, with averages H and X and variances σ H and σ X respectively. For the analysis it is convenient to work with dimensionless variables. To this end we divide all the parameters entering the utility by a quantity that plays the role of unit of measurement1 , that we denote β. Thus we write

H i = β(h + y i ) (7)
where H ≡ βh stands for the population's average willingness-to-join the community. Thus, y i is a dimensionless random variable; its probability density function (pdf) f Y (y) has zero mean and variance σ H /β. We denote the parameters of the model in units of β with small characters:

h = H β , j = J β , g = G β , c = C β , x i = X i β , (8)
and correspondingly,

a i = A i /β, b i = B i /β.
Hereafter we consider the equilibria corresponding to an iterative best response dynamics: starting from an arbitrary initial condition each agent chooses the strategy that maximizes his surplus. Agents may either make their decisions simultaneously (parallel dynamics) or one after the other in an arbitrary order (random sequential dynamics). In order to get some insight about the problem, it is useful to consider a plane whose axes are the quenched random variables x (abscissas) and y (ordinates), where each individual i is represented by a point according to his values (x i , y i ). Then, for any given values of η c and η f , the lines

y=xη c -h-(j+g)η c -jη f do not join (a i <0 ; a i +b i <0) free-ride (a i >0 ; b i <0 ) y x cooperate (b i >0 ; a i +b i >0) x=c/η c y=c-h-(j+g)η c -jη f
y = y m ≡ c -h -(j + g)η c -jη f , (9a)
y = y m -c + xη c , (9b)
x = x m ≡ c η c , (9c)
represent the marginal individuals whose utilities are at the boundaries between different optimal strategies. These lines partition the (x, y) plane into three regions, represented on figure 1. The optimal strategy of individuals whose (x i , y i ) values lie in the upper-right region, defined by

x i > x m or b i > 0, ( 10a
)
y i > y m or a i + b i > 0, (10b)
is to cooperate. The best strategy for those whose (x i , y i ) satisfy

x i < x m or b i < 0, ( 11a
)
y i > y m -c + x i η c or a i > 0, (11b)
is to free-ride. Individuals with (x i , y i ) in the lower-right region do not join the community.

Mean field dynamics

We are interested in the temporal evolution of the system within an iterated game setting, in which individuals have to choose their best strategies s i (t) based on the available information. We assume that each time an agent has to make a decision, he has the exact information of the global proportions of cooperators and free-riders at the preceding outcome, η c (t-1) and η f (t-1), and uses these quantities to estimate his surplus [START_REF] Ledyard | Public goods: A survey of experimental research[END_REF]. Such a dynamics is called Cournot best reply in economics literature. In simulations, starting from initial guesses η c (0) and η f (0), the updating is said to be in parallel if all the individuals in the population first determine their best strategies based on the preceding outcome, and make their decisions simultaneously afterwards. At the opposite, in random sequential updating, a single individual selected at random is asked to make his decision at each time step. The latter dynamics simulates systems where the individuals make their decisions without any temporal correlation. In order to compare the time scales of both dynamics, it is usual to consider that N sequential time steps are equivalent to one parallel update. Intermediate updating schemes may be implemented, but here we only consider these two extreme cases, that are standard in economics and in statistical physics.

Referring back to figure 1, since the boundary lines depend on the values of η c (t) and η f (t), they will shift in the course of updating according to the perceived proportions of cooperators and free-riders. Individuals whose values of y i and x i lie close to the boundaries are susceptible to small changes of η c and η f : their strategies may change with time, and in turn induce changes in those of the others.

Hereafter we study how the system reaches its stable states upon successive updates, that is, the path followed by a point representative of the system's state in the plane (η c , η f ). Such paths should either end up in a stable fixed point or get trapped in other types of attractors if they exist. In the following we consider separately the parallel and sequential updating schemes, since the corresponding dynamic equations are different.

Let us introduce the complementary cumulative functions defined by:

G χ (ζ) ≡ 1 -F χ (ζ) = ∞ ζ f χ (ξ)dξ, (12)
where χ ∈ {X, Y }, and F χ (ζ) is the cumulative distribution of f χ . Using the fact that y is independent of x, the probabilities p and q that an individual i taken at random will cooperate or free-ride respectively at the next time step given values of η c (t) and η f (t) are

p(s i (t + 1) = 1) ≡ p(η c (t), η f (t)) = ∞ xm(t) f X (x)dx ∞ ym(t) f Y (y)dy = G X (x m (t)) G Y (y m (t)), (13a)
q(s i (t + 1) = -1) ≡ q(η c (t), η f (t)) = xm(t) -∞ ∞ ym(t)-c+xηc(t) f X (x)f Y (y)dydx = xm(t) 0 f X (x) G Y (y m (t) -c + xη c (t)) dx, (13b)
where y m (t) and x m (t) are given by equations (9c) and (9a) respectively with η c = η c (t) and η f = η f (t).

Fixed point equations

We consider the limit of very large populations, namely N → ∞, for which η c and η f are unbiased estimators of the probabilities of cooperate and freeride respectively. Since at a fixed point the probabilities of the possible strategies are stationary, they satisfy equations ( 13) with p = η c , q = η f :

η c = p(η c , η f ), (14a)
η f = q(η c , η f ). ( 14b
)
Notice that the trivial solution η c = 0 satisfies (14a), corresponding to a state where nobody contributes to the public good. The community survives because when η c = 0 free-riders do not suffer from any moral burden. The model in this case reduces to the model with (symmetric) social interactions j considered in [START_REF] Gordon | Discrete choices under social influence: Generic properties[END_REF]. Here η f plays the role of the fraction of customers, and may present a single or two fixed points depending on whether the interaction strength j is smaller or larger, respectively, than a critical value j B = 1/f Y (y max ) where y max is the maximum of f Y , the IWJ distribution. Notice however that Nash equilibria without cooperators are a drawback of our keeping the model simple enough to be studied analytically. The main question, namely whether polymorphic organizations with coexistence of cooperators and free-riders are stable, is meaningful outside the region η c = 0.

Parallel dynamics equations: two-dimensional map

In parallel dynamics, one assumes that the agents update their strategies simultaneously, starting from an arbitrary initial configuration. The perceived η c and η f used by all the agents to estimate their utilities are the same. The boundary lines between the three regions of figure 1, which depend on η c (t) and η f (t), partition the population at time t according to their estimated utilities: all the individuals whose idiosyncratic parameters lie in region I will choose to cooperate, those in region II to free-ride and those in region III not to join. In this way, using equations ( 13), the dynamics can be formulated as the following two-dimensional deterministic map Ω :

(η ′ c , η ′ f ) = (p(η c , η f ), q(η c , η f )) in the simplex S := {(η c , η f ) : η c ≥ 0, η f ≥ 0, η c + η f ≤ 1}, namely: η c (t + 1) = p(η c (t), η f (t)), (15a)
η f (t + 1) = q(η c (t), η f (t)). ( 15b
)
Thus, the fixed points ( 14) are the fixed points of the parallel dynamics.

The results presented in this section shed light on their nature and allow to determine if there are other kinds of equilibria.

In two-dimensional maps or flows, a fixed point is called a sink if all the points in its neighborhood converge to it, a source if these points diverge from it and a saddle if in one direction the map or flow converges, while in the other directions it diverges. Together these three types are called hyperbolic, which is the only kind of equilibria possibly encountered in a structurally stable system, i.e. a system stable with respect to small variations of the parameters (in our case h, j, g, c and d). For our two-dimensional map, the nature of an equilibrium point (η c0 , η f 0 ) is determined by the following Jacobian matrix

J m = ∂p(η c0 , η f 0 )/∂η c ∂p(η c0 , η f 0 )/∂η f ∂q(η c0 , η f 0 )/∂η c ∂q(η c0 , η f 0 )/∂η f .
Denoting the two eigenvalues of J m as ω 1 and ω 2 , then the fixed point is a sink if

|ω 1 | < 1 and |ω 2 | < 1, a source if |ω 1 | > 1 and |ω 2 | > 1 and a saddle if (|ω 1 | -1)(|ω 2 | -1) < 0 ([17] §1.4).
There are two types of curves emanating from each saddle point x 0 , namely the stable manifold, defined as

W s (x 0 ) = {x ∈ R 2 |Ω t (x) → x 0 as t → ∞} (16)
and the unstable manifold, defined as

W u (x 0 ) = {x ∈ R 2 |Ω -t (x) → x 0 as t → ∞} (17)
Ω -1 denoting a backward iterate of the map. The boundaries of the basins of attraction are usually formed by stable manifolds. Moreover, if a transversal intersection between the stable and unstable manifolds of the same saddle point exists, the map will exhibit chaotic behavior.

Random sequential dynamics equations: two-dimensional flow

In sequential dynamics, one updates one agent chosen at random in each time step. It is usual to consider N (the number of agents in the system) successive updates (called one Monte-Carlo step) as being comparable to one step of parallel updating, i.e., the time step between two successive individual updatings is τ = 1/N . In each individual update, the expected displacement of the pair (η c , η f ) is

η c (t + τ ) = η c (t) + (1/N )(p(η c (t), η f (t)) -η c (t)), (18a)
η f (t + τ ) = η f (t) + (1/N )(q(η c (t), η f (t)) -η f (t)), (18b)
where p and q are defined in [START_REF] Gordon | Cooperation and free-riding with moral cost[END_REF], and the factor 1/N is the probability of selecting the agent that is updated. Taking the continuous limit N → ∞, we have the following set of differential equations:

dη c dt = p(η c , η f ) -η c , (19a)
dη f dt = q(η c , η f ) -η f , (19b)
with the time unit being one Monte-Carlo step. Now the system evolves as an autonomous system in a planar phase space, i.e. as a two-dimensional flow in S. Clearly, the fixed points of the sequential dynamics are the same as the equilibria [START_REF] Durlauf | How can statistical mechanics contribute to social science[END_REF].

The three generic types of equilibria, namely source, sink and saddle exist in a structurally stable two-dimensional flow as well. The nature of an equilibrium point (η c0 , η f 0 ) is now determined by the Jacobian matrix of the flow

J f = ∂p(η c0 , η f 0 )/∂η c -1 ∂p(η c0 , η f 0 )/∂η f ∂q(η c0 , η f 0 )/∂η c ∂q(η c0 , η f 0 )/∂η f -1 .
Denoting the eigenvalues of

J f as ω ′ 1 and ω ′ 2 , then the point is a sink if ω ′ 1 < 0 and ω ′ 2 < 0, a source if ω ′ 1 > 0 and ω ′ 2 > 0 and a saddle if ω ′ 1 ω ′ 2 < 0 ([17] §1.2∼1.3
). Note that the condition for an equilibrium to be a saddle gives the same inequality as in the parallel case. Therefore, the set of saddle points in sequential updating coincides with that in parallel updating. However, the set of sinks or sources are not necessarily the same in both cases, since the governing inequalities are quite different. The definitions of the stable and unstable manifolds are analog to those in the two-dimensional map with the discrete time step replaced by the continuous time variable t. According to the Poincaré-Bendixson Theorem ([17] Theorem 1.8.1), the two-dimensional flow system will never go into chaos.

4 Equilibria for particular distributions f X and f Y To go further with the analysis of the preceding section we need to specify the pdfs of the idiosyncratic parameters. In this paper we present the most interesting of the cases we have studied, in which the x i follow a uniform distribution of finite width d:

f X (x) = 1 d for 0 ≤ x ≤ d, (20a)
f X (x) = 0 otherwise. ( 20b
)
and the y i are distributed according to:

f Y (y) = 1 4 cosh 2 (y/2) . ( 21
)
The cumulative function corresponding to [START_REF] Alligood | Chaos: An Introduction to Dynamical Systems[END_REF] is the logistic distribution.

The complementary functions [START_REF] Durlauf | Statistical mechanics approaches to socioeconomic behavior[END_REF] are:

G X (z) = (1 -z)Θ(1 -z), (22a)
G Y (z) = 1/[1 + exp(z)], (22b)
where Θ is the Heaviside function.

Fixed points

The fixed point equations ( 14) may be written as follows:

η c = (1 - ρ η c )Θ(1 - ρ η c ) G Y (c -Z) (23a)
η f = ρ cη c log[1 + (e c -1)G Y (c -Z)], (23b)
where ρ ≡ c/d and Z ≡ h + j(η c + η f ) + gη c . Since f X has a bounded support, equation (23a) vanishes if 0 ≤ η c (t) ≤ ρ. Thus, if η c (0) ≤ ρ, the system evolves towards the trivial equilibrium without cooperators. In parallel dynamics, equation (15a) gives η c (1) = 0: the state with no cooperators is reached after a single parallel update. Afterwards, on the axis η c = 0, η f evolves according to equation (15b),

η f (t + 1) = G Y (-h -jη f (t)). ( 24
)
The equilibrium value of η f is given by (14b), obtained by replacing η f (t) and η f (t + 1) by η f in the above equation. As already mentioned (see [START_REF] Gordon | Discrete choices under social influence: Generic properties[END_REF]) there is a critical value of j, j B , such that for j < j B , there is a single equilibrium. In that case all the points η c (0) < ρ will eventually be mapped to it. If j > j B (14b) has three solutions. One of them, η f u , is an unstable fixed point separating the basin of attraction of the two others, that are stable. To summarize, when η c (0) < ρ, after a first time step the system is mapped to the axis η c = 0 and then evolves, following (24), either to one fixed point (if j < j B ) or to one of two fixed points (if j > j B ) depending on whether η f (0) > η f u or η f (0) < η f u . In the case of the logistic distribution f Y considered here, j B = 4 [START_REF] Nadal | Multiple equilibria in a monopoly market with heterogeneous agents and externalities[END_REF]. Under random sequential updating, since p(η c (0), η f (0)) = 0 for 0 ≤ η c (0) ≤ ρ, equation (19a) shows that η c vanishes exponentially fast, the fixed points being the same as for the parallel dynamics. These solutions have been described here for completeness since, as already stated, they are not within the scope of the questions addressed by the model Notice that there are no fixed points with η f = 0: since G Y ≥ 0, the right hand side of (23b) vanishes only if c = 0, i.e. full cooperation may exist only if cooperation is costless.

When η c (0) > ρ, calling G the value taken by G Y (c -Z), equation (23a) gives η c = (1 -ρ ηc ) G, which can be solved for non-vanishing η c in terms of G:

η c = η ± c [G] ≡ 1 2 G {1 ± [1 - 4ρ G ] 1/2 } ( 25
)
Later we shall see that both the + and the -branches may give stable equilibria, though the parameter range for which a stable equilibrium exists for the -branch is much narrower than for the + branch. From (25) one gets that the Nash equilibria with cooperators satisfy η c ≥ 2ρ (equality can occur when j = g = 0): the fraction of cooperators is larger than a threshold given by the cooperation cost relative to the average weight (x = d/2) of the social disapproval for free-riding. Equation (23b) can also be parameterized in terms of G by defining

η ≡ j(η c + η f ) + gη c j + g . (26)
Introducing ( 25) and (23b) into (26) we obtain an equation for η:

η = η 1 [G] ≡ η ± c [G] + j j + g 1 η ± c [G] 1 d log[1 + (e c -1)G] (27) Inverting (22b) for z = c -h -(j + g)η, we obtain η in terms of G: η = η 2 [G] ≡ 1 j + g {c -h -log 1 -G G }. ( 28
)
The possible solutions η are then obtained by the intersects of the curves given by equations ( 20) and ( 21). Equilibria with parallel dynamics.

Notice that the above analysis allows to determine the fixed points of the system, but doesn't give any hint about the existence of cycles. We have thus numerically computed the stable and unstable manifolds for each saddle point, with the same set of values for j, g, c and d as in figure 2, while changing h -c from -4 to -3 in steps of 0.02. We summarize here the main results, leaving the details to appendix A.1.

Besides the trivial fixed point(s) with η c = 0, for any h -c ≥ -3.4201 there is a fixed point O + corresponding to a relatively large organization having a proportion of cooperators larger than the proportion of free-riders. This occurs for cooperation costs relatively large compared to the average willingness to join (h -c is negative). The value of j considered is large enough to compensate negative values of h i -c on a large fraction of the population, producing a bandwagon effect. For -3.3504 < h -c < -3.1146, this fixed point O + coexists with another non-trivial attractors (see figure 3). For -3.3504 < h -c < -3.2726, the attractor is a limit cycle C S . It shrinks (at h -c = -3.2726) to a second fixed point O ′ + which corresponds to an organization with qualitatively similar proportions of cooperators and free riders, and with fewer members than at the equilibrium O + . This second fixed point O ′ + exists for -3.2726 < h -c < -3.1146. For larger values (h -c > -3.1146) it disappears, leaving only the fixed point O + with a very large basin of attraction.

Figures 4 and5 present the results for two values of h -c in the regions where we expect coexistence of the two non-trivial attractors. On the left hand side are represented the different manifolds predicted by the mean field theory (details lefts to Appendix A.1). On the right hand side we present the dynamical paths obtained under parallel dynamics updating of a system with N = 1000 agents with the same parameters, starting from different initial conditions. The agreement with the analytical results is very good, although there are some finite size effects (see the figure's legend). One of the most striking result of the simulations are the winding trajectories that turn around over long times before reaching the corresponding attractor. This means that an oscillatory behavior might be observed in actual systems even before convergence to the cycle attractor. Notice also that these oscillations may exist even if the system is expected to reach a Nash equilibrium.

Equilibria with random sequential dynamics

We have also computed the stable and unstable manifolds for the flow, the mean field approximation for sequential dynamics. The results, presented in appendix A.2 are summarized on figure 3. They are similar to those with parallel dynamics, although for the parameters considered (j = 5, g = 1, c = 4.5 and d = 72) the system only presents fixed points. Cycles exist, but The non trivial attractors are a stable cycle C s (its basin of attraction is limited by the dashed line corresponding to the stable manifold of H + , with its two ends at the line η f = 0). The point O + is a fixed point. Its basin of attraction is the region bounded by the stable manifold of H -between the lines η c + η f = 1 and η f = 0, excluding the basin of attraction of C s .

Right. Trajectories of the simulated system for four different initial conditions, that converge respectively to a trivial fixed point (squares), to a cycle of length 41 (empty triangles), to a fixed point (full circles) not predicted by the mean field equations (due to finite size effects) and to the non-trivial fixed point O + (empty circles). 3, with h -c = -3.2. The stable limit cycle C s has merged with the source X + through a Hopf bifurcation to become a sink O ′ + , whose basin of attraction is the stable manifold of the saddle H + . Right. Trajectories of the simulated system (parallel dynamics) for five different initial conditions. One of them converges to a trivial fixed point (squares), one (circles) to a cycle of length 13 (instead of the predicted non-trivial fixed point O ′ + ), the three others (triangles) converge to the non-trivial fixed point O + , which has a large basin of attraction.

for other parameter values, as discussed in the appendix. Figure 6 presents an example of evolution of the same system under parallel and sequential dynamics, showing the dramatic difference between both types of updating rules.

Simulations results and phase diagrams

The theoretical results of the preceding sections predict the behavior of an infinite system. On figure 7 we present the phase portrait obtained through parallel dynamics with systems of different sizes, for h -c = -3.3 and the same set of parameters as in figure 4. Starting from different initial conditions, we determined the attractors in the phase space (η c , η f ). As N decreases, more and more fixed points invade the region where cycles exist, and many different cycles appear. This multiplicity of attractors is due to missing values of the IWJ in the finite population. However, the mean field predictions, and in particular the existence of cycles, remain qualitatively correct.

In order to get deeper insight on the phase diagram of the model, we studied the attractors of large systems as a function of the different parameters. Results as a function of h -c, summarized on figure 8, are in very good agreement with the mean field predictions, as may be seen through the position of the theoretical boundaries λ n .

Figures 9 present separately the η c values at the fixed points and at the cycles for systems of different sizes, ranging from N = 1 000 to N = 100 000, as a function of h -c. Consistently with the above results, the size effects are important, but do not modify qualitatively the phase diagram.

As may be seen on figure 10 (left), the fraction of the phase space that

Discussion

In order to understand the existence of the limit cycles, we analyzed in details the dynamics of a particular system with N = 100 individuals. We isolated on figure 11 the agents involved in the limit cycle, identified by the values of their parameters x i and y i . As expected, these values are close to the marginal values (9), i.e. they lie close to the phase boundaries defined in figure 1. Notice that these boundaries depend on the values of η c and η f , so that their positions change on time. In the course of the cycle, some individuals oscillate between states s i = 1 and s i = 0, others between s i = 1 and s i = -1, and others between s i = -1 and s i = 0. Some even go successively through the three possible states, as may be seen in the inset of figure 11. that oscillate between cooperating or not to participate at all (outsiders), in black those that oscillate between cooperating or free-riding, and in blue those that either free-ride or do not participate. Notice (see inset) that some individuals may adopt alternatively each of the three possible states.

y i x i 1 0 0 1 -1 1 1 -1 -1 0 0 -1 16 20 24 1.8 2.0 2.2
Intuitively, the dynamics of the cycle is as follows: assume that at some instant t (like at t = 26 in figure 12) the system has proportions η c (t) and η f (t). Then, as more cooperators enter the organization, the cost of freeriding becomes too high for some agents in the organization (those with relatively large x i ). These will leave the organization, decreasing thus the total number of members. As a consequence, the value of the organization decreases for everybody and both the fractions of cooperators and of free-riders at the following periods will decrease. At some point, certain individuals may join the organization as free-riders, because the social pressure due to cooperators has decreased. As a result, the total number of members increases making it worth for some cooperators to join the organization, initiating a new oscillation. As may be seen on figure 12, the cycle is actually closed (in the mathematical sense) only after many oscillations, when the fractions of cooperators and free-riders repeat themselves exactly. It is interesting that the total fraction of members also oscillates, but with a different phase, revealing the complexity of the dynamics. In the computer simulations, the length of the cycles depend very strongly on the particular realization of the random variables y i and x i . Figures 13 show examples for different system sizes. Notice that the length of the cycles depend on the system's details (the particular realization of the random variables x i and y i ) even for systems with 1 million agents. Finally, let us point out that our results remain qualitatively unchanged if G = 0. In contrast, the existence of cycles requires J > 0, as well as the two idiosyncratic components, namely, the willingness to join H i and the moral burden weight X i . In fact, it is straightforward to show that for J = 0 the dynamics converges necessarily to fixed points.

Conclusion

In this paper we studied in details a model of a charity organization, or any other organization whose value is proportional to the number of its members [START_REF] Phan | Adoption and cooperation in communities: mixed equilibrium in polymorphic populations[END_REF]. Individuals may not participate, and those joining the organization may or not cooperate. Thus, the size of the organization is an emergent property, not fixed a priori.

We investigated whether such a voluntary organization may survive, and under what conditions. Similar situations have been studied within the public goods literature. According to the standard game theoretical assumption of selfish individuals, such an organization is expected to be exclusively composed of free-riders. However, daily life experience and experimental economics results have shown that some level of cooperation often exists, a fact attributed to moral feelings. The model considers heterogeneous individuals; they have idiosyncratic preferences for joining the community and they also weight idiosyncratically the social disapproval for free-riding. The surplus of the community members is proportional to their number, be they cooperators or not, but cooperators contribute more than free riders to the community value. Upon joining, cooperators bear a fixed cost while free-riders suffer a cost due to cooperators' disapproval, that they weight idiosyncratically.

We have shown that in this model there are three kinds of equilibrium attractors. There is a trivial Nash equilibrium for all the parameters of the model, in which the organization is exclusively composed of free-riders. Thus, the population is partitioned into two categories: the outsiders (those that do not join the community) and the free-riders. Such an equilibrium between only two strategies exists because in the model the community has a value for its members even if nobody cooperates. In this extreme case there is no cost for free-riding, and the parameter X i does not play any role. The system's properties are the same as those of simple bandwagon models where the individuals have to make a binary choice (to buy or not a good at a given price, participate or not to a riot, etc). It has been shown [START_REF] Gordon | Discrete choices under social influence: Generic properties[END_REF] that for some values of the pertinent parameters (h -the population average willingness to join-, and j -the strength of the bandwagon effect -), there may be one or two fixed points. Coexistence of two possible collective states (Nash equilibria), having either a large or a small fraction of free-riders, arise for values of h smaller than, and j larger than, threshold values that depend on the distribution of the idiosyncratic component H i . Besides this trivial equilibrium, the model considered here presents one or two Nash equilibria with finite fractions of cooperators and free-riders.

An interesting result of this model is that besides the Nash equilibria, it exhibits limit cycles with corresponding fractions of cooperators, of free-riders (and their sum -the fraction of the population that joins the organization-) oscillating through time. This phenomenon has its root in the strong asymmetry of the interaction between agents, introduced by the idiosyncratic weights that free-riders attach to social disapproval. Although we studied in details the case where these weights X i are uniformly distributed, we verified that the same kind of oscillations exist in the case of a bimodal distribution2 . The fixed points in this case may be found analytically.

We studied in details the dynamics of the cycles under parallel and random sequential updating both analytically and numerically. Assuming a global neighborhood (i.e. every agent's surplus depends on the fractions of cooperators and free-riders in the whole population), we obtained results which are exact in the asymptotic limit N → ∞. The two-dimensional dynamical maps and flows have been studied for a large range of parameters. Although qualitatively similar, the phase diagrams of maps and flows corresponding to the same parameters are not identical. The fact that a rich class of bifurcations can occur within a narrow parameter region is in itself truly remarkable. A natural question to ask is what happens in intermediate regimes, where subsets of agents make decisions simultaneously, while others make them at independent times.

Numerical simulations of finite size systems with parallel dynamics qualitatively agree with the analytical results. The dynamical behaviour of the simulated systems show winding trajectories in very large regions of the phase diagram. Oscillations in the fraction of cooperators and free-riders are thus expected in such systems before the onset of the cycle and even when the attractor is a fixed point. We are currently performing simulations with sequential dynamics, where the notion of cycle is not well defined, due to the randomness in the order in which the agents' decisions are updated. Preliminary results also exhibit the above mentionned oscillatory behavior, but the analysis of the data is much more cumbersome.

Future work should focus on the dynamics for interactions restricted to closer neighborhoods, which probably requires more sophisticated techniques from the theory of dynamical systems. Another interesting question is how the phase diagram is modified in the case where the agents learn from past actions and can form expectancies. The algorithms used to compute the stable and unstable manifolds can be found in [START_REF] England | Computing onedimensional stable manifolds and stable sets of planar maps without the inverse[END_REF] and Chapter 10 of [START_REF] Alligood | Chaos: An Introduction to Dynamical Systems[END_REF], respectively. The results classified according to the qualitative dynamical features are shown in figures [START_REF] Durlauf | How can statistical mechanics contribute to social science[END_REF].

A Appendix

Hereafter we consider a system with parameters j = 5, g = 1, c = 4.5, d = 72, and provide the increasing values of h -c, denoted by λ i , at which the successive bifurcations appear, together with the qualitative descriptions of their nature (figure 3 of section 4 summarizes the results).

• λ 1 = -4.0971. For h -c < λ 1 there exists no intersection between η 1 (G) and η 2 (G) and all the points in S are eventually mapped to the trivial fixed point, which is always a solution of the equations. At h -c = λ 1 , η 2 (G) begins to intersect η - 1 (G) at two points, one of which is a saddle (H -) and the other is a source (X -). The process of simultaneous creation or elimination of a pair of equilibrium points is called a saddle-node bifurcation. Note that this bifurcation does not alter the overall behavior of the system, since no new attractor is created. Thus, only the fixed points η c = 0, η f > 0 exist, and since j > j B , the system may flow to either of the two fixed points for η f depending on the initial conditions.

• λ 2 = -3.4201. Another saddle-node bifurcation happens on η + 1 (G) at h-c = λ 2 , with a saddle (H + ) and a sink (O + ) created simultaneously. Now the stable manifolds of H + and H -possess a common end emanating from X + , and together they divide S into two regions. Points outside the region bounded by the stable manifolds should converge to the trivial fixed point η c = 0, η f > 0 as before, and those inside to O + .

• λ 3 = -3.3617. For λ 2 < h -c < λ 3 , the unstable manifold of H + approaches the stable manifold of H -until they coincide at h -c = λ 3 (the common manifold is called a separatrix between H + and H -). Beyond this value, the unstable manifold of H + folds back and converges to O + while the stable manifold of H -hits the simplex boundary η f = 0. Now the basin boundary is determined solely by the stable manifold of H -, causing the basin of attraction of O + to experience a sudden expansion.

• λ 4 = -3.3504. For λ 3 < h -c < λ 4 , the stable and unstable manifolds of H + approach each other until they coincide at h -c = λ 4 . Beyond this value, the former ends up at the boundary η f = 0 of S and the latter folds back into a stable limit cycle C s around X + . At this transition C s is introduced into the system as a new attractor with its own basin of attraction delimited by the stable manifold of H + .

• λ 5 = -3.2726. For λ 4 < h -c < λ 5 , C s shrinks and eventually merges with X + at h -c = λ 5 to produce a sink O ′ + . The process of transition between a sink and a source with the simultaneous appearance or disappearance of a limit cycle is called a Hopf bifurcation.

• λ 6 = -3.1146. For λ 5 < h -c < λ 6 , O ′ + and H + approach each other until they are annihilated at h -c = λ 6 . We are left with H -and O + , the stable manifold of H -serving as the basin boundary. This topology persists for arbitrarily larger values of h -c. Summarizing the analysis above, we can identify three types of transition:

1. saddle-node bifurcation; 2. Hopf bifurcation; 3. separatrix bifurcation between two saddle points, or one saddle point with itself.

Type 1 can be directly determined from the number of intersections between η 1 (G) and η 2 (G) (cf. figure 2), or from the fact that J m has an eigenvalue +1. Due to the lack of symmetry in our system, transcritical and pitchfork bifurcations never occur ([17] §3.4). In case that J m has an eigenvalue -1, one can have a period-doubling bifurcation ([17] §3.5). This never happens in our map, since one can easily verify that

J m (1, 1) > 0, J m (2, 2) > 0, ∆(J m ) > 0, (A-1)
which excludes the possibility that J m has a negative eigenvalue. Types 2 and 3 cannot be observed from the intersections between η 1 (G) and η 2 (G) and must be determined numerically. Type 2 remains a local bifurcation with J m having a pair of conjugate complex eigenvalues with unit modulus. In contrast, Type 3 is a global bifurcation with the area of the basin of attraction experiencing a sudden jump or, in physicist's terms, a first-order transition. Numerically we have found no evidence for homoclinic intersection leading to chaos, but the question as to whether chaos exists in our two-dimensional map remains open.

Figures 15 are results of numerical simulations corresponding to the same parameters as in figure 14.

A.2 Phase structure of the flows

We have computed the stable and unstable manifolds as for the map. The results are shown in Figures 16 (see figure 3 of section 4 for a summary). Hereafter we provide the bifurcation values of h -c, denoted by µ i , together with the qualitative descriptions of the nature of the bifurcation:

• µ 1 = λ 1 = -4.0971. As in the parallel case, a saddle node bifurcation occurs at h -c = µ 1 and produces a saddle H -and a source X -. • µ 2 = -3.6172. A Hopf bifurcation occurs at h -c = µ 2 and produces an unstable limit cycle C u together with a sink O -from the source X -. The basin of attraction of O -is the area surrounded by C u . Note that we've been able to choose h -c such that the sink remains in the negative branch, which implies that variation of η with respect to h-c does not always qualify to determine the stability of the fixed point.

• µ 3 = -3.4980. A separatrix bifurcation of the stable and unstable manifolds of H -eliminates C u .

• µ 4 = λ 2 = -3.4201. A saddle node bifurcation occurs at h -c = µ 4 and produces a sink O + and a saddle H + .

• µ 5 = λ 6 = -3.1146. A saddle node bifurcation occurs at h-c = µ 5 and annihilates the saddle H + and the source O ′ + . The topology persists for larger h -c. Thus, we have the same types of bifurcation as in the parallel case, and Peixoto's Theorem ( [START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems, and bifurcations of vector fields[END_REF] Theorem 1.9.1) indicates that these are the only types of bifurcation expected in a two-dimensional flow system. However, let's note that the overall topology of the dynamics depends strongly on their order of occurrence. For example, if the separatrix bifurcation happens before the Hopf bifurcation, a stable limit cycle can be created just like in the parallel case (see the figure at the bottom (right) of figures 16).

Figure 1 :

 1 Figure 1: Boundaries between regions corresponding to the different individual choices, as a function of the values of the (dimensionless) quenched random variables y i and x i .

η 1 [Figure 2 :

 12 Figure 2: Fixed points. Solution by curve intersection for distributions f X and f Y

Figure 3 :

 3 Figure 3: Regions of existence of the different stable equilibria under parallel and random sequential dynamics (see appendices A.1 and A.2 for the values of the λ n and the µ n ).

Figure 4 :

 4 Figure 4: Left. Stable and unstable manifolds of the theoretical saddle points in the plane (η c , η f ) for the same parameter values as in figure 3, with h -c = -3.3.The non trivial attractors are a stable cycle C s (its basin of attraction is limited by the dashed line corresponding to the stable manifold of H + , with its two ends at the line η f = 0). The point O + is a fixed point. Its basin of attraction is the region bounded by the stable manifold of H -between the lines η c + η f = 1 and η f = 0, excluding the basin of attraction of C s . Right. Trajectories of the simulated system for four different initial conditions, that converge respectively to a trivial fixed point (squares), to a cycle of length 41 (empty triangles), to a fixed point (full circles) not predicted by the mean field equations (due to finite size effects) and to the non-trivial fixed point O + (empty circles).

Figure 5 :

 5 Figure 5: Left. Theoretical equilibria for the map in the plane (η c , η f ), for the same parameter values as in figure3, with h -c = -3.2. The stable limit cycle C s has merged with the source X + through a Hopf bifurcation to become a sink O ′ + , whose basin of attraction is the stable manifold of the saddle H + . Right. Trajectories of the simulated system (parallel dynamics) for five different initial conditions. One of them converges to a trivial fixed point (squares), one (circles) to a cycle of length 13 (instead of the predicted non-trivial fixed point O ′ + ), the three others (triangles) converge to the non-trivial fixed point O + , which has a large basin of attraction.

Figure 6 :

 6 Figure 6: Examples of trajectories for parallel and random sequential dynamics starting from the same initial conditions, for a system with the same parameter values as in figure 3, and h -c = -3.3.

Figure 7 :

 7 Figure 7: Attractors and their basins of attraction for different system sizes, obtained by starting from all the possible initial conditions (with a mesh of 0.01). Parameters are: j = 5, g = 1, d = 72, h = 1.2. The white region is the basin of attraction of the cycles.

Figure 8 :

 8 Figure 8: Attractors for parallel dynamics with N = 100 000, for 20 simulated systems, with the same parameters as in the preceding figures, i.e. j = 5, g = 1, c = 4.5 and d = 72. The points in this and the following figures are obtained as follows: we draw 20 systems at random with probabilities f X and f Y . The points in the figures correspond to the final states obtained for each system, starting from each possible initial point in the plane η c , η f laying on a grid of mesh 0.1 satisfying the condition η c + η f ≤ 1. That corresponds to 45 different initial conditions for each of the 20 systems, i.e. a total of 900 simulations.

Figure 9 :Figure 10 :

 910 Figure 9: Fixed points (left) and cycles (right) obtained through simulations of systems of different sizes.

Figure 11 :

 11 Figure 11: Parameters (x, y) space. Individuals involved in a cycle: in red those

20

 20

Figure 12 :

 12 Figure 12: Cycles. Composition of the system for N = 1 000.

Figure 13 :

 13 Figure 13: Cycles. Total fraction of community members as a function of t for different system sizes (left) and for three different systems with N = 10 6 .

A. 1

 1 Phase structure of the maps

Figure 14 :

 14 Figure 14: The stable and unstable manifolds for increasing values of h -c. For h -c = -3.7 only the trivial fixed points with η c = 0 exist. For h -c = -3.4, H + and O + have been created from a saddle node bifurcation on η+ 1 (G). At h -c = -3.36 a separatrix bifurcation between the stable manifold of H -and the unstable manifold of H + has introduced a first order transition in the basin of attraction. At h -c = -3.3 (reproduced from section 4) a separatrix bifurcation between the stable and unstable manifolds of H + has introduced a stable cycle C s with its own basin of attraction (the interior of the line C s ). At h -c = 3.2 (reproduced from section 4) the stable limit cycle C s has merged with the source X + through a Hopf bifurcation to become a sink O ′ + . At h -c = -3.1 a final saddle node bifurcation has annihilated O ′ + and H + .

Figure 15 :Figure 16 :

 1516 Figure 15: Examples of trajectories of a simulated system (N = 1000), starting from different initial conditions, corresponding to the maps of figure 14.

We may choose any convenient unit, like H, σH , homogeneous to a utility.

Note added in proof: A similar model for an asset market with two kinds of agents, that corresponds to a bimodal distribution of the Xi, restricted to binary decisions[START_REF] Belitsky | Stability analysis with applications of a two-dimensional dynamical system arising from a stochastic model for an asset market[END_REF] was shown to present cycles.

Acknowledgements

M.B.G. and J-P.N. are members of CNRS. S.G. thanks support from FAPERGS and CNPq (Brazil) and also from CNRS during his stay in Grenoble.


D Papo

J M Buldú

Brain synchronizability, a false friend

Synchronization plays a fundamental role in healthy cognitive and motor function. Two issues are still poorly understood. First, how synchronization depends on the interplay between local dynamics, coupling and topology? Second, how prone to synchronization a network is, given its topological organization? To investigate the synchronizability of both anatomical and functional brain networks various studies resorted to the Master Stability Function (MSF) formalism, an elegant tool which allows analysing the stability of synchronous states in a dynamical system consisting of many coupled oscillators. Here, we argue that brain dynamics does not fulfil the formal criteria under which synchronizability is usually quantified and, perhaps more importantly, this measure refers to a global dynamical condition that never holds in the brain (not even in the most pathological conditions), and therefore no neurophysiological conclusions should be drawn based on it. We discuss the meaning of synchronizability and its applicability to neuroscience and propose alternative ways to quantify brain networks synchronization.

Introduction

Consider a network in which each node is a dynamical system, e.g. an oscillator, and the links are couplings between these nodes. Can these oscillators synchronize with each other creating a coherent state and, if so, under what circumstances is this state stable? Given a particular dynamical system and coupling scheme, the Master Stability Function (MSF) formalism [START_REF] Pecora | Master stability functions for synchronized coupled systems[END_REF][START_REF] Boccaletti | Complex networks: structure and dynamics[END_REF][START_REF] Arenas | Synchronization in complex networks[END_REF] allows relating the stability of the fully synchronized state to the spectral properties of the underlying matrix of connections, and assessing which network structures can maintain complete synchronization of the whole network.

At the macroscopic scales of typical non-invasive neuroimaging techniques, brain activity can be thought of as the collective dynamics of a set of coupled dynamical units. Synchronization among these units has been suggested to be a basic mechanism of healthy brain functioning [START_REF] Varela | The brainweb: phase synchronization and large-scale integration[END_REF]. Thus, at first glance, the problem above may seem to apply to brain activity, justifying the use of the MSF formalism to quantify brain network synchronizability. However, all formalisms are created to address very specific questions and come with their own set of formal and theoretical assumptions, the compliance with which ultimately decides whether they can be used in a given context.

In the remainder, we argue that some essential characteristics of the brain render the MSF framework difficult to apply to neuroscience, review some misunderstandings about the synchronizability construct, and propose alternative ways to understand synchronization in brain networks.

Brain synchronizability

The use of synchronizability, initially designed to study theoretical models, rapidly extended to the analysis of real datasets and, in the context of neuroscience, to quantify the ability of anatomical [START_REF] Chavez | Dynamics of excitable neural networks with heterogeneous connectivity[END_REF][START_REF] Zhao | Competition between intracommunity and inter-community synchronization and relevance in brain cortical networks[END_REF][START_REF] Ton | Structure-function discrepancy: inhomogeneity and delays in synchronized neural networks[END_REF][START_REF] Phillips | Graph theoretic analysis of structural connectivity across the spectrum of Alzheimer's disease: the importance of graph creation methods[END_REF][START_REF] Tang | Developmental increases in white matter network controllability support a growing diversity of brain dynamics[END_REF] and functional [START_REF] De Haan | Disruption of functional brain networks in Alzheimer's disease: what can we learn from graph spectral analysis of resting-state magnetoencephalography?[END_REF][START_REF] Bassett | Adaptive reconfiguration of fractal small-world human brain functional networks[END_REF][START_REF] Reijneveld | The application of graph theoretical analysis to complex networks in the brain[END_REF][START_REF] Stam | Graph theoretical analysis of complex networks in the brain[END_REF][START_REF] Schindler | Evolving functional network properties and synchronizability during human epileptic seizures[END_REF][START_REF] Deuker | Reproducibility of graph metrics of human brain functional networks[END_REF][START_REF] Van Wijk | Comparing brain networks of different size and connectivity density using graph theory[END_REF][START_REF] Jalili | EEG-based functional networks in schizophrenia[END_REF][START_REF] Van Dellen | MEG network differences between low-and high-grade glioma related to epilepsy and cognition[END_REF][START_REF] Tahaei | Epilepsy synchronizability of EEG-based functional networks in early Alzheimer's disease[END_REF][START_REF] Bialonski | Assortative mixing in functional brain networks during epileptic seizures[END_REF][START_REF] Lehnertz | Evolving networks in the human epileptic brain[END_REF][START_REF] Niso | What graph theory actually tells us about resting state interictal MEG epileptic activity[END_REF][START_REF] Khambhati | Virtual cortical resection reveals push-pull network control preceding seizure evolution[END_REF] brain networks to synchronize. For example, Tang and co-workers [START_REF] Tang | Developmental increases in white matter network controllability support a growing diversity of brain dynamics[END_REF] investigated how the human brain's anatomical organization evolves from childhood to adulthood by measuring changes in the synchronizability parameter, and proposed that during the course of development human brain anatomy evolves towards an organization that limits synchronizability [START_REF] Tang | Developmental increases in white matter network controllability support a growing diversity of brain dynamics[END_REF]. The authors suggested that as the brain evolves towards its mature state, it reduces its ability to synchronize, and that this reduction would help increase brain controllability. Furthermore, a few studies focused on the effects of different pathologies on brain synchronizability, such as epilepsy [START_REF] Schindler | Evolving functional network properties and synchronizability during human epileptic seizures[END_REF][START_REF] Van Dellen | MEG network differences between low-and high-grade glioma related to epilepsy and cognition[END_REF][START_REF] Tahaei | Epilepsy synchronizability of EEG-based functional networks in early Alzheimer's disease[END_REF][START_REF] Bialonski | Assortative mixing in functional brain networks during epileptic seizures[END_REF][START_REF] Lehnertz | Evolving networks in the human epileptic brain[END_REF][START_REF] Niso | What graph theory actually tells us about resting state interictal MEG epileptic activity[END_REF][START_REF] Khambhati | Virtual cortical resection reveals push-pull network control preceding seizure evolution[END_REF], Alzheimer's disease [START_REF] Phillips | Graph theoretic analysis of structural connectivity across the spectrum of Alzheimer's disease: the importance of graph creation methods[END_REF][START_REF] De Haan | Disruption of functional brain networks in Alzheimer's disease: what can we learn from graph spectral analysis of resting-state magnetoencephalography?[END_REF], or schizophrenia [START_REF] Jalili | EEG-based functional networks in schizophrenia[END_REF], showing statistically significant changes in the synchronizability parameter in association with these diseases. Interestingly, epilepsy was associated with an increased synchronizability during interictal activity [START_REF] Lehnertz | Evolving networks in the human epileptic brain[END_REF], while it decreased during ictal activity [START_REF] Schindler | Evolving functional network properties and synchronizability during human epileptic seizures[END_REF]. Functional networks synchronizability has been reported to decrease the electroencephalographic (EEG) activity of schizophrenic patients [START_REF] Jalili | EEG-based functional networks in schizophrenia[END_REF]. Studies using magnetoencephalography (MEG) showed that synchronizability values depend on the frequency band considered when constructing functional networks [START_REF] Bassett | Adaptive reconfiguration of fractal small-world human brain functional networks[END_REF].

While changes in synchronizability clearly exist, is this particular metric measuring what it is supposed to measure?

The Master Stability Function formalism

The meaning and scope of the synchronizability construct should be understood in the MSF theoretical context it is predicated upon.

Given a group of N diffusively coupled dynamical systems whose dynamics in isolation follows = F x , the evolution of the whole system is given by the equation:

= -∑ [ x -x ], = 1, … , [1]
where is the m-dimensional state vector of the ith oscillator, the coupling strength, x a vectorial output function and the elements of the adjacency matrix A, with and = 1 if system (i.e., node) i and j are connected and zero otherwise. Since oscillators are coupled in a diffusive way, i.e., with a linear coupling proportional to the difference between their state variables, we can introduce into equation [START_REF] Pecora | Master stability functions for synchronized coupled systems[END_REF] the Laplacian matrix defined as = -, where D is a diagonal matrix whose elements ! are the degree (i.e., the number of neighbours) of each node [START_REF] Boccaletti | Complex networks: structure and dynamics[END_REF]. Equation [START_REF] Pecora | Master stability functions for synchronized coupled systems[END_REF] then reads

= -∑ ! x , = 1, … , [2]
and the coupling term is basically dependent on the coupling strength and the elements of the Laplacian matrix ! [24].

Here we consider the Laplacian matrix to be symmetric (which ensures that its eigenvalues are real), as it simplifies the definition of the MSF and, in turn, it is the most common situation in the neuroscience literature using the MSF methodology. For identical systems with the same coupling function x , the synchronized state is a solution of " = # , 2 with = $ = ⋯ = = " . Next, a linear stability analysis around the synchronization manifold can be carried out. This involves checking whether the Lyapunov exponents corresponding to phase space directions transverse to the synchronization manifold are all negative. To do so, a deviation δ = -" from the synchronous state " is introduced and the linear stability of the whole system evaluated, obtaining an m-dimensional variational equation, m being the number of variables of the dynamical system x [1-3]:

' = [( " -)( " ]', [ 3
]
where J is the Jacobian operator and the independent variable ) is defined as ) * = σ, * , , * being the non-zero eigenvalues of the Laplacian matrix. Finally, the stability of the synchronization manifold can be obtained by studying the parametrical behaviour of the maximum Lyapunov , -./ ν as function of ). 1 ) = , -./ ν with ν ∈ [0, ∞ is referred to as the master stability function [START_REF] Pecora | Master stability functions for synchronized coupled systems[END_REF][START_REF] Boccaletti | Complex networks: structure and dynamics[END_REF]. Importantly, the MSF as a pure function of ν is independent of the network of connections between the dynamical systems. However, considering a specific network of connections determines particular values of ) * , which can be related to the MSF, ultimately determining the conditions for the stability of the synchronization manifold of that specific network. Overall, the MSF tells how (i) the dynamical system, through and (ii) the network topology, through the second term on the right side of equation [START_REF] Boccaletti | Complex networks: structure and dynamics[END_REF] concur in determining the stability of the synchronization manifold. Note that, the term synchronizability refers to the stability of the global synchronization state. The synchronization manifold is stable when all ) * 5 associated with the non-zero eigenvalues of the Laplacian matrix lie in a region in which the MSF is negative. However, different dynamical systems with different coupling functions lead to qualitatively different MSFs (see Fig. 1a for details), which can be classified as [START_REF] Boccaletti | Complex networks: structure and dynamics[END_REF]: class I (always positive), class II (always negative above a threshold) and class III (negative only within a specific region). Interestingly, in the context of brain networks, synchronizability is commonly evaluated as if the brain were a class III system, although no proof of it exists. Thus, the lower the ratio 6 between the largest and smallest (non-zero) eigenvalues (i.e. 6 = λ λ $ ⁄ , the more packed the eigenvalues of the Laplacian are and the highest the ability to fall within a window where the MSFs is always negative. In that sense, brain networks' synchronizability is sometimes [START_REF] Tang | Developmental increases in white matter network controllability support a growing diversity of brain dynamics[END_REF] quantified by the dispersion of the eigenvalues of the Laplacian matrix 8 [25].

Synchronizability: some common misconceptions

The meaning of synchronizability and the questions it allows addressing are a matter of frequent confusion and numerous misconceptions.

An important issue is whether synchronizability can be measured when ignoring the characteristics of the dynamics. Stability under perturbations exists when all eigenvalues of the combinatorial Laplacian matrix 9λ : fall within the region of stability due to the fact that the coupling is strong enough to guarantee that the MSF enters the region but weak enough to guarantee that it does not leave this region from the other side. Thus, the value of the coupling strength σ is crucial to determine ) * once the values of , * are obtained from the network structure. Synchronizability is ultimately determined by the sign of the MSF evaluated at points that are indeed given by the spectrum of the Laplacian matrix and an overall coupling strength. The functional form of the MSF crucially depends on the dynamics of the coupled oscillators and the function that couples its state variables to those of other oscillators [START_REF] Nishikawa | Network synchronization landscape reveals compensatory structures, quantization, and the positive effect of negative interactions[END_REF][START_REF] Huang | Generic behavior of masterstability functions in coupled nonlinear dynamical systems[END_REF]. Depending on the shape of the MSF, dynamical systems may never synchronize, always synchronize above a certain coupling strength or synchronize only for coupling strength values within a certain range [START_REF] Boccaletti | Complex networks: structure and dynamics[END_REF][START_REF] Huang | Generic behavior of masterstability functions in coupled nonlinear dynamical systems[END_REF]. While the MSF for various families of dynamical systems is typically convex for generic oscillator systems, its exact shape depends not only on the dynamical systems but also on the kind of coupling between them. Thus, quantifying the synchronizability of anatomical brain networks using a parameter based on the eigenvalues of the Laplacian matrix alone, without information about the underlying dynamical oscillators and their coupling function and strength cannot ensure that the whole system falls within MSF's stability region. In other words, it is not the network structure per se that is synchronizable, but the particular combination of dynamical systems, coupling strength and network structure formed by the connections between these systems. While the eigenvalues of the Laplacian matrix likely contain potentially valuable information of some sort [START_REF] Phillips | Graph theoretic analysis of structural connectivity across the spectrum of Alzheimer's disease: the importance of graph creation methods[END_REF][START_REF] De Haan | Disruption of functional brain networks in Alzheimer's disease: what can we learn from graph spectral analysis of resting-state magnetoencephalography?[END_REF], eigenvalue dispersion of the anatomical network alone without at least some information on node dynamics cannot determine the system class one is dealing with, and conclusions on its MSF are no more than guesses (see Fig. 1a). Furthermore, the stability region may not be unique [START_REF] Huang | Generic behavior of masterstability functions in coupled nonlinear dynamical systems[END_REF], and the MSF classes may not be limited to the three depicted in Fig. 1a, further complicating the interpretation of MSF theory's 6.

Two related important questions are: what are high or low synchronizability values telling us? When can synchronizability values be compared? The bare comparison of synchronizability values across dynamical systems and the characterization of a given topology as being more or less synchronizable than another are potentially problematic: insofar as different dynamical systems haven't necessarily got similar MSFs, the synchronizability parameter of a brain network cannot be compared with others as long as its MSF is unknown.

Perhaps at the root of most other ones, a major problematic issue relates to the frequent confusion between synchronizability and synchronization. The synchronizability parameter does not tell if the system is synchronized or not: a system can be highly synchronizable without being synchronized, and synchronized with a low synchronizability parameter (see Fig. 1b).

Finally, it is worth stressing that the synchronizability construct only applies to anatomical networks. This is because the MSF formalism relies on a structural property, i.e. the connectivity pattern between dynamical units, which should be complemented by the coupling strength. However, the construction of functional networks relies on the reported coordination between brain regions, i.e. a dynamical property. Therefore, functional networks are not the cause of a certain level of synchronizability, but their consequence. ) is related to the (N-1) non-zero eigenvalues λ of the network Laplacian matrix as ) = λ where σ is the coupling strength. The synchronization manifold is stable when all ) lie in a region where the MSF is negative. MSFs can be classified as [START_REF] Boccaletti | Complex networks: structure and dynamics[END_REF]: class I (always positive), class II (always negative above a threshold ) $ ) and class III (negative only within a specific region [) $ , ) ]). Note that, when the MSF has several minima and maxima, the stability region may not be unique [START_REF] Huang | Generic behavior of masterstability functions in coupled nonlinear dynamical systems[END_REF]; these MSF classes are not depicted in the figure . (a The interval [) $ , ) ] contains the eigenvalues of the Laplacian matrix corresponding to a given network. Note that in this particular example only for the class II system are values of the MSF corresponding to the set of ) all negative, indicating that the synchronization manifold is locally stable. (b) Counter-example showing that defining synchronizability parameter as the inverse of the dispersion of the eigenvalues can be misleading: network A has lower dispersion (i.e., higher synchronizability) but lies in the region of the MSF where the synchronization manifold is unstable, while network B, has higher dispersion (and lower synchronizability) but can synchronize. In this qualitative example, we consider the Laplacian matrix to be symmetric.

Why synchronizability should not be used (in neuroscience)

Even discounting the technical issues discussed above, fundamental reasons make the MSF-based synchronizability inapplicable to neuroscience.

Crucially, in its original formulation [START_REF] Pecora | Master stability functions for synchronized coupled systems[END_REF], the MSF applies to diffusively-coupled identical dynamical systems, i.e. all interacting units of the network should have the same variables and internal parameters. However, irrespective of the scale at which it is observed, the brain is dynamically highly heterogeneous, ruling out an application of the MSF. While the MSF formalism can be generalized to heterogeneous systems, this comes at the price of rather restrictive conditions hampering its application to brain data [START_REF] Nishikawa | Network synchronization landscape reveals compensatory structures, quantization, and the positive effect of negative interactions[END_REF][START_REF] Sun | Master stability functions for coupled nearly identical dynamical systems[END_REF][START_REF] Zhang | Identical synchronization of nonidentical oscillators: when only birds of different feathers flock together[END_REF]. No less importantly, the condition of diffusivity, is also grossly violated in the brain.

Perhaps the most fundamental obstacle to the use of the MSF in brain sciences is represented by two issues related to the definition of synchronization. First, while various kinds of synchronization, including phase [START_REF] Varela | The brainweb: phase synchronization and large-scale integration[END_REF], generalized [START_REF] Stam | Synchronization likelihood: an unbiased measure of generalized synchronization in multivariate data sets[END_REF], and relay synchronization [START_REF] Vicente | Dynamical relaying can yield zero time lag neuronal synchrony despite conduction delays[END_REF] have been reported for brain dynamics, and may even coexist [START_REF] Malagarriga | Consistency of heterogeneous synchronization patterns in complex weighted networks[END_REF], synchronizability refers to a specific synchronization mode, complete synchronization. Complete synchronization requires that all dynamical units have exactly the same phase and amplitude once the synchronization manifold is reached, a state that has never been reported in the brain (not even in its most pathological conditions). Second, physics and neuroscience understand synchronization in fundamentally different ways: in the former, synchronization refers to a global and stable state, while in the latter to a local and transient one. While local complete synchronization may be a relevant mechanism or a reasonable modelling representation of functionally segregated regions or circuits, its dynamics is necessarily transient. Brain dynamics has in general a complex phase space geometry, and possibly no stricto sensu attractor at all [START_REF] Rabinovich | Transient dynamics for neural processing[END_REF][START_REF] Tognoli | The metastable brain[END_REF], a scenario that cannot be dealt with using the MSF in its current form.

Towards neurophysiologically plausible alternatives to synchronizability

Several technical systems, e.g. power grids, wireless communication systems, require stable synchronization of their units [START_REF] Kinzel | On chaos synchronization and secure communication[END_REF][START_REF] Tyrrell | Emergent slot synchronization in wireless networks[END_REF][START_REF] Rohden | Self-organized synchronization in decentralized power grids[END_REF][START_REF] Rohden | Impact of network topology on synchrony of oscillatory power grids[END_REF][START_REF] Motter | Spontaneous synchrony in power-grid networks[END_REF]. Synchronizability represents, to good approximation, a construct that can be used to model and regulate their dynamics and function.

However, synchronizability refers to a type of synchronization that the brain does not, certainly should not, and possibly cannot achieve in a stable way. In addition to being incompatible with the dynamical and functional heterogeneity of normal brain functioning, a completely synchronized state represents a considerable loss of complexity, and would likely be associated with an unphysiological energetic cost [START_REF] Torrealdea | Energy efficiency of information transmission by electrically coupled neurons[END_REF][START_REF] Moujahid | Energy and information in Hodgkin-Huxley neurons[END_REF].

Before figuring out possible alternatives to the MSF-based synchronizability, one should perhaps address the following question: why is the MSF framework used although it so evidently at odds with neurophysiological stylized facts? What makes the MSF a convenient tool? While a unique coupling function for all network nodes and some hypotheses on the coupling matrix are convenient mathematical conditions which ensure the existence of an invariant set representing the complete synchronization manifold and considerably simplify the analysis of its stability, using steady state dynamics and complete synchronization dispenses with defining spatial topography and temporal scales of the target process.

To figure out possible alternatives to the MSF-based synchronizability construct one needs to understand both the role played by synchronization within this conceptual framework and the objective pursued by the studies using it and the problems that they aim to address. On the one hand, while in neuroscience synchronization typically refers to transient bivariate coupling between two neuronal ensembles, the synchronization referred to by synchronizability is in fact better thought of as a (steady-state) process on a network. On the other hand, from a teleological view-point, resorting to the MSF formalisms can be understood in terms of the need to address the relationship between structure (or, more precisely, the topology defined on it) and dynamics in complex systems [START_REF] Arenas | Synchronization in complex networks[END_REF][START_REF] Skardal | Optimal synchronization of complex networks[END_REF][START_REF] Menck | How basin stability complements the linear-stability paradigm[END_REF][START_REF] Ponten | The relationship between structural and functional connectivity: graph theoretical analysis of an EEG neural mass model[END_REF]. Furthermore, this relationship must be understood at all meaningful spatial scales of the brain [START_REF] Kelso | Outline of a general theory of behavior and brain coordination[END_REF]. Given an observed dynamics and topological organization, a construct teleologically equivalent to synchronizability may possibly be framed in terms of a networked system's propensity to enter a functionally desirable state or regime, as in [START_REF] Wang | On partial contraction analysis for coupled nonlinear oscillators[END_REF]. But what dynamical states or regimes may represent a valuable target, the distance from which may be used as a neurophysiologically meaningful benchmark?

The true difficulty in finding alternatives to the MSF-based synchronizability is that one loses the uniqueness and taskinvariance of the complete synchronization state and needs to cope with brain dynamics' spatial heterogeneity and temporal multiscaleness, and brain function's translational invariance. To define a valid equivalent of synchronizability will likely require three key ingredients: neurophysiologically plausible and functionally meaningful order parameters describing collective brain activity; mechanisms through which they may emerge; and, no less importantly, those through which they may wane. On the one hand, this should for instance involve considering networks of heterogeneous oscillators and plausible synchronization processes, compatible with metastable dynamics [START_REF] Tognoli | The metastable brain[END_REF][START_REF] Deco | The dynamics of resting fluctuations in the brain: metastability and its dynamical cortical core[END_REF][START_REF] Roberts | Metastable brain waves[END_REF]. On the other hand, the mechanisms through which neural assemblies interact and their role in human brain function at various scales of brain anatomy and dynamics should be better understood at both functional/computational and algorithmic/dynamical levels. These mechanisms are likely context-specific, and various ones may even coexist [START_REF] Malagarriga | Consistency of heterogeneous synchronization patterns in complex weighted networks[END_REF]. As a consequence, the definition of a dynamical target may vary as a function of the putative role of synchrony lato sensu in the target activity. On the one hand, dynamical references may be associated with different dynamical regimes, e.g. cluster synchronization, in which patterns or sets of synchronized elements emerge, chimera-like states, i.e. spatiotemporal patterns with coexisting coherence and incoherence, or partial synchronization, where only some parts of the network synchronize, while others do not [START_REF] Zhou | Hierarchical synchronization in complex networks with heterogeneous degrees[END_REF][START_REF] Abrams | Solvable model for chimera states of coupled oscillators[END_REF][START_REF] Bi | Coexistence of quantized, time dependent, clusters in globally coupled oscillators[END_REF]. Importantly, these dynamical regimes would induce spatial scales. On the other hand, rather than statistically stationary states, what is needed is an analysis of their dynamics, stability, bifurcations, and symmetries [START_REF] Abrams | Solvable model for chimera states of coupled oscillators[END_REF][START_REF] Pecora | Cluster synchronization and isolated desynchronization in complex networks with symmetries[END_REF]. Importantly, a reference regime should also replicate the temporal scales of some (task-specific or task-independent) reference brain activity. The construct may for instance contain predictive information on the properties of and on the conditions under which these clusters form and dissolve. Defining meaningful dynamical target processes and some sort of distance from them to given observed ones, understanding whether and the extent to which these may emerge from interactions between local dynamics and network topology are all highly non-trivial but fundamental questions, finding answers to which will likely keep the neuroscience community busy for some time to come.

Concluding remarks

We have argued that not only is the synchronizability construct an inadequate tool to quantify brain networks' ability to synchronize, but the problem itself to which it is supposed to provide an answer appears to be ill-posed when studying brain dynamics. More generally, the brain differs in many essential ways from the systems (e.g. the electrical power-grid or the Internet) most network theory constructs were originally designed to account for. Neuroscience, a field where network theory has only relatively recently come to the foreground [START_REF] Bullmore | Complex brain networks: Graph theoretical analysis of structural and functional systems[END_REF], has so far mainly borrowed its tools and concepts without inspiring fresh theory, and this has exposed it to the risks inherent in such an application: over-, under-and misuse of existing tools [START_REF] Papo | Functional brain networks: great expectations, hard times, and the big leap forward[END_REF][START_REF] Papo | Beware of the Small-World neuroscientist! Front[END_REF]. The brain's unique properties can help promoting a fundamental reformulation of network neuroscience for benefit of both neuroscience and graph theory.

24. An adaptation to weighted networks can easily be obtained by including the weights of the connections in and replacing the node degree by the node strength (i.e., the sum of the weights of node's links). 25. The dispersion of the eigenvalues of the Laplacian matrix is given the following expression:

1 5 $ = $ -1 ∑ =λ -λ > = $ ? ⁄ ⁄ with λ > = 1 -1 ⁄ ∑ λ ?

Fig 1 .

 1 Fig 1. Master Stability Function Ψ(ν) as a function of the parameter ν. ) is

)

  Qualitative example of ) of a network that would synchronize class I and II dynamical systems, but not class III (assuming thatλ λ $ ; ⁄ ) ) $ . ⁄

Acknowledgements -The author acknowledges financial support from the program Accueil de

Talents of the Métropole Européenne de Lille and from the Labex (laboratory of excellence) DISTALZ (Development of Innovative Strategies for a Transdisciplinary approach to ALZheimer's disease).