Consensus and diversity in multi-state noisy voter models

TL;DR

This paper analyzes a multi-opinion noisy voter model with all-to-all interactions, revealing a noise-driven transition between consensus and mixed opinion states, with exact and approximate analytical characterizations of the stationary distribution.

Contribution

It introduces an effective death-birth process to analytically study the stationary distribution and transition dynamics in multi-state noisy voter models.

Findings

Identifies a noise-driven transition between consensus and mixed states.

Provides exact solutions for homogeneous rates and approximations for heterogeneous rates.

Derives mean switching times between consensus states.

Abstract

We study a variant of the voter model with multiple opinions; individuals can imitate each other and also change their opinion randomly in mutation events. We focus on the case of a population with all-to-all interaction. A noise-driven transition between regimes with multi-modal and unimodal stationary distributions is observed. In the former, the population is mostly in consensus states; in the latter opinions are mixed. We derive an effective death-birth process, describing the dynamics from the perspective of one of the opinions, and use it to analytically compute marginals of the stationary distribution. These calculations are exact for models with homogeneous imitation and mutation rates, and an approximation if rates are heterogeneous. Our approach can be used to characterize the noise-driven transition and to obtain mean switching times between consensus states.

Equations

$$X_i+X_j$$

$$X_j$$

$${\partial }_{t}P\left(\mathbf{n},t\right)=\sum _{i=1}^m\sum _{j=1}^m\left({\mathcal{E}}_i^{-1}{\mathcal{E}}_j-1\right)\left[T_{j\to i}(\mathbf{n})P\left(\mathbf{n},t\right)\right],$$

$$T_{j\to i}(\mathbf{n})=r_{ji}\frac{n_i{n}_j}{N}+{ϵ}_{ji}{n}_j.$$

$$T_i^{+}(\mathbf{n})\equiv \sum _{j\ne i}{T}_{j\to i}(\mathbf{n}),$$

$$T_i^{-}(\mathbf{n})\equiv \sum _{j\ne i}{T}_{i\to j}(\mathbf{n}).$$

$$T_i^{+}(\mathbf{n})$$

$$T_i^{-}(\mathbf{n})$$

$$P_i(n_i)=\frac{{\prod }_{k=1}^{n_i}\frac{T_i^{+}(k-1)}{T_i^{-}(k)}}{1+{\sum }_{k=1}^N{\prod }_{\ell =1}^k\frac{T_i^{+}(\ell -1)}{T_i^{-}(\ell )}}.$$

$$t_j^{0\to N}=\sum _{k=0}^{N-1}\sum _{\ell =0}^k\frac{1}{T_j^{+}(\ell )}\prod _{{\ell }^{\prime }=\ell +1}^k\frac{T_j^{-}({\ell }^{\prime })}{T_j^{+}({\ell }^{\prime })}.$$

$${\tau }_i=\frac{1}{{\sum }_{j\ne i}(t_j^{0\to N}{)}^{-1}}.$$

$$\tau =\sum _{i=1}^m{p}_i{\tau }_i,$$

$${\overleftrightarrow{\tau }}_i\equiv {\tau }_i+t_i^{0\to N}.$$

$$p_i=\frac{({\overleftrightarrow{\tau }}_i{)}^{-1}}{{\sum }_{j=1}^m({\overleftrightarrow{\tau }}_j{)}^{-1}}.$$

$$\tau =\sum _i\left\{\frac{({\overleftrightarrow{\tau }}_i{)}^{-1}}{\left[{\sum }_j({\overleftrightarrow{\tau }}_j{)}^{-1}\right]\left[{\sum }_{j\ne i}(t_j^{0\to N}{)}^{-1}\right]}\right\},$$

$$T_i^{+}(n_i)=$$

$$T_i^{-}(n_i)=$$

$$N_L=\frac{\left[1+{\tilde{\epsilon }}_{\mathrm{tot}}\right]+\sqrt{{\left[1+{\tilde{\epsilon }}_{\mathrm{tot}}\right]}^2-4\tilde{\epsilon }}}{2\tilde{\epsilon }},$$

$$N_R=\frac{\left[1+\tilde{\epsilon }\right]+\sqrt{{\left[1+\tilde{\epsilon }\right]}^2-4{\tilde{\epsilon }}_{\mathrm{tot}}}}{2{\tilde{\epsilon }}_{\mathrm{tot}}}.$$

$$\tau =\frac{t^{0\to N}}{m-1}.$$

$$T_i^{+}(\mathbf{n})$$

$$T_i^{-}(\mathbf{n})$$

$${⟨r⟩}_{-i}(\mathbf{n})=\frac{1}{N-n_i}\sum _{j\ne i}{n}_j{r}_j,$$

$$r_{-i}^{\ast }\equiv \frac{1}{1-x_i^{\ast }}\sum _{j\ne i}{x}_j^{\ast }{r}_j,$$

$${\dot{x}}_i=f_i(\mathbf{x})+\frac{1}{\sqrt{N}}{\xi }_i(t),$$

$$f_i(\mathbf{x})=\sum _j[{\mathcal{T}}_{j\to i}(\mathbf{x})-{\mathcal{T}}_{i\to j}(\mathbf{x})].$$

$$⟨{\xi }_i(t){\xi }_j(t^{\prime })⟩=B_{ij}(\mathbf{x})\delta (t-t^{\prime }),$$

$$B_{ij}(\mathbf{x})=\left\{\begin{array}[]{cc}\sum_{k\neq i}\left[{\cal T}_{i\to k}(\mathbf{x})+{\cal T}_{k\to i}(\mathbf{x})\right]&\mbox{if\leavevmode\nobreak\ }i=j\\ \nobreak\leavevmode\hfil\\ -\left[{\cal T}_{j\to i}(\mathbf{x})+{\cal T}_{i\to j}(\mathbf{x})\right]&\mbox{if\leavevmode\nobreak\ }i\neq j\end{array}\right..$$

$$f_i(\mathbf{x})=\sum _{j\ne i}\left[x_i{x}_j(r_{ij}-r_{ji})+{\epsilon }_{ji}{x}_j-{\epsilon }_{ij}{x}_i\right],$$

$$B_{ij}(\mathbf{x})=-x_i{x}_j\left(r_{ij}+r_{ji}\right)-{\epsilon }_{ij}{x}_i-{\epsilon }_{ji}{x}_j,$$

$$B_{ii}(\mathbf{x})=-\sum _{j\ne i}{B}_{ij}(\mathbf{x}).$$

$$\frac{\partial {\mathcal{P}}_i(x_i)}{\partial t}=$$

Full text

I Introduction

The Voter Model (VM) was initially introduced as a model for spatial conflict Clifford and Sudbury (1973); Liggett (1985). In its most basic variant, the model describes a population of voters with pairwise interactions. Each individual can be in one of two states. In an interaction, one individual copies the state of another. If the population is finite, this process continues until all individuals have reached the same state; no further dynamics are then possible.

This model has been studied in the context of different applications. For example, variants of the VM have been used to describe autocatalytic reactions Frachebourg and Krapivsky (1996), herding in financial markets Kirman (1993), opinion dynamics Castellano et al. (2009), and the evolution of language Castelló et al. (2006); Castellano et al. (2009). In the context of biological evolution, the VM is closely related to the Moran process with neutral selection Granovsky and Madras (1995). In each of these applications the states an individual can take represent different properties. They can be opinions on a given issue, trading behaviors in a model of a financial market, or different species in a biological context.

Although the models used in these applications share common characteristics, different research communities have focused on different aspects of the VM and its applications. For the purpose of opinion dynamics, for example, one may be interested in whether or not a population of agents reaches consensus and, if they do, how long this takes. In economics one may ask how copying of trading behaviour leads to herding, and in the context of linguistics one might be interested in understanding how features of languages spread and organise in space Kauhanen et al. (2018). Applications of the VM include the use to model actual elections Fernández-Gracia et al. (2014).

The VM is also of interest from a point of view of statistical mechanics. The model is similar to spin models in traditional statistical physics, except that it is defined dynamically and not through an energy function. As such, the VM has become one of the most studied models of non-equilibrium statistical physics. The two-state VM in two dimensions has been shown to exhibit logarithmic coarsening Dornic et al. (2001). This is in contrast for example with Glauber dynamics for the Ising model, in which coarsening is driven by surface tension and growth laws are algebraic. Multi-state voter models, on the other hand, can show algebraic coarsening Dall’Asta and Galla (2008). The traditional VM presents an interesting class of dynamics, with two absorbing states and $\mathbb{Z}_{2}$-symmetry. It has been shown to define its own universality class Dornic et al. (2001); Al Hammal et al. (2005), and field theories have been devised to study its critical point Dornic et al. (2005). A further aspect attracting attention in the physics community has been the coupled dynamics of the voter model and interaction networks between agents. Such co-evolutionary processes have, been shown to lead to fragmentation transitions Holme and Newman (2006); Vazquez et al. (2008), enriching the number and nature of the absorbing states of the VM.

One interesting variant in the class of voter dynamics is the so-called ‘noisy voter model’ Granovsky and Madras (1995). The terminology might at first sound surprising – the process of opinion changes through interaction in the original VM is already stochastic. However, it allows for absorbing states, in which no further dynamics can occur.

What is meant by the term ‘noisy’ in this context is that, in addition to the interaction with other individuals, random changes of opinion can occur; they do not require interaction with another agent, and are sometimes also described as ‘mutations’. These random changes drive the system away from consensus; there are no absorbing states.

As a consequence, two effects compete in the noisy voter model: a drive towards consensus through interaction, and spontaneous opinion changes promoting coexistence. This leads to a transition between a regime in which the system is mostly ordered (i.e., all agents are of the same opinion), and another regime in which both opinions are represented in the system. This phenomenon is known as ‘noise-induced bistability’ and has been investigated in the context of autocatalytic reactions, surface-reaction models, decision making of insects, and biochemical reactions Fichthorn et al. (1989); Togashi and Kaneko (2001); Ohkubo et al. (2008); Biancalani et al. (2012, 2014); Houchmandzadeh and Vallade (2015); Saito et al. (2016). Recently, transitions of this type have also been studied in noisy voter models on complex networks Carro et al. (2016); Peralta et al. (2018a, b).

Most of the existing work on models of this type focuses on the case in which individuals can take two different states; the transition is then between a phase in which the stationary distribution of individuals is unimodal and another in which it is bimodal; hence the term ‘noise-induced bistability’.

The aim of the present work is to generalise the model to the case of multiple opinion states. In the absence of mutation there are then multiple absorbing states. As we demonstrate, the inclusion of spontaneous state changes leads to noise-induced multi-stability. We investigate this numerically and analytically. To do this, we compute marginals of the stationary distribution of the model, as well as switching times between the different consensus states. Our analysis focuses on the case of a ‘mean-field’ geometry, that is, a population in which all pairs of individuals can interact at all times.

The remainder of the paper is organized as follows. In Sec. II we define the model and introduce the general notation. In Sec. III we derive an effective master equation for the dynamics of a single species and use it to describe marginals of the stationary distribution of the multi-species model; closed-form analytical expressions are obtained. This approach is exact if all species have the same imitation and mutation rates, and an approximation otherwise. We also estimate switching times between different consensus states. Sec. IV focuses on the model with homogeneous imitation and mutation rates. We obtain analytical predictions for its phase diagram and test these predictions against simulations. In Sec. V we carry out a similar analysis for a multi-state model with heterogeneous rates. Sec. VI contains a discussion of our results and an outlook on future work.

II Definition of the model and noise-induced multi-stability

II.1 Model definition and notation

We consider a population of $N$ individuals; each can be of one of $m$ types ($m\in\mathbb{N}$). In the context of social dynamics these would represent different opinions; in evolutionary biology they may stand for different species. Our aim is not to study a specific application, but the general structure that emerges from multi-state noisy voter models; we will therefore use the terms opinions, species and types interchangeably.

The population is unstructured, that is, any individual can interact with any other member of the population at all times. We write $n_{i}$ for the number of individuals of species $i$, where ${i\in\{1,\dots,m\}}$; these variables vary as the population evolves. At each time the state of the system is fully specified by the vector ${\mathbf{n}=(n_{1},\dots,n_{m})}$. The size of the population is constant, i.e., we have $\sum_{i=1}^{m}n_{i}=N$ at all times; the state space is a simplex in $m$-dimensions.

The dynamics of the model are defined by the following imitation and mutation reactions:

[TABLE]

The notation $X_{i}$ represents an individual of species $i$, and the first reaction describes the process in which an individual of opinion $j$ imitates an individual of opinion $i$. We write $r_{ji}$ for the corresponding imitation rate; the exact interpretation of these coefficients will be made more precise below in Eq. (3). In the context of evolutionary dynamics reactions of this type represent a death-birth event: an individual of type $j$ dies, and is replaced by the offspring of an individual of species $i$. This offspring is also of type $i$. To keep the model general, we allow this rate to depend on both species involved, $j$ and $i$.

The second reaction in Eq. (II.1) describes spontaneous changes of opinion, i.e., an individual of type $j$ turns into an individual of type $i$. The rate with which this process occurs is $\varepsilon_{ji}$, where we again include possible dependence on $j$ and $i$. In the context of population dynamics it may be more realistic to incorporate mutation in the reproduction events; that is, one could consider models in which the offspring of an individual does not necessarily belong to the same species as the parent. The effect of both ways of introducing random state changes is similar though: they allow departure from states in which the entire population is made up of individuals of one single species. We will refer to instances of the second process in Eq. (II.1) as ‘mutation events’.

The dynamics are defined in continuous time, and can be described by the master equation

[TABLE]

where $P(\mathbf{n},t)$ is the probability of finding the population in state $\mathbf{n}$ at time $t$. We have written ${\cal E}_{i}$ for the creation operator for individuals of species $i$; it is defined by its action ${\cal E}_{i}f(n_{1},\dots,n_{i},\dots,n_{m})=f(n_{1},\dots,n_{i}+1,\dots,n_{m})$ on functions $f(\mathbf{n})$ of the state of the population. The quantity $T_{j\to i}(\mathbf{n})$ is the rate with which individuals of type $j$ are converted into individuals of type $i$ if the system is currently in state $\mathbf{n}$; these rates are given by

[TABLE]

We note that this is the overall rate for such events in the population, and not a per capita rate. As it is common practice, the scaling with the population size, $N$, is such that $T_{j\to i}={\cal O}(N)$. This ensures that ${\cal O}(N)$ events occur per unit time; in other words, time is measured in generations, rather than individual events.

Among the states of the population there are states of complete consensus, in the language of opinion dynamics. This is the case when $n_{i}=N$ for one opinion $i$, and $n_{j}=0$ for all $j\neq i$. In the context of evolutionary dynamics these are monomorphic states; one species has taken over the entire population. These consensus states are not absorbing when mutation rates $\varepsilon_{ji}$ are non-zero. As we will discuss next, it is the interplay between imitation, mutation and the associated intrinsic noise that generates the interesting behaviour of the model.

II.2 Noise-induced multistability

Simulations of the model can be carried out efficiently using the standard Gillespie algorithm Gillespie (1976, 1977). Different types of outcomes are illustrated in Fig. 1. We show the time evolution of the fraction of individuals in each species, $x_{i}=n_{i}/N$, in individual simulations of models with $m=2$ and $m=5$ species. In panels (a)–(c) the population is of size $N=50$, and the population is seen to visit states of complete consensus relatively frequently. These are not permanent; instead, the population switches between different consensus states. In panels (d)–(f) the population size is $N=500$. Consensus may occasionally be reached, but other than that the system is mostly found in states of diversity, where different opinions are represented at fluctuating frequencies. Panels (c) and (f) show the multi-modal and unimodal stationary distributions for a model with $m=3$.

III Reduction to effective single-opinion dynamics

We proceed to characterise the multi-state noisy VM analytically. To do this, we first derive an effective death-birth process for the individuals of a particular species. We then study the shape of the stationary distribution resulting from these effective dynamics. We also use the single-species dynamics to obtain information about typical switching times between consensus states.

III.1 Effective single-species master equation

We focus on a particular species $i$ and the dynamics of the number of individuals $n_{i}$ of this species. Events in the population dynamics convert individuals of one species into individuals of another. Many of these events will not involve species $i$, and are therefore not relevant for the effective dynamics of species $i$. In the following, we only consider events which either increase or decrease $n_{i}$.

The overall rate with which individuals are converted to species $i$ is

[TABLE]

where $T_{j\to i}(\mathbf{n})$ is the conversion rate from $j$ to $i$, as defined in Eq. (3). Similarly, the rate with which individuals of species $i$ are converted into any other type is

[TABLE]

We note that these rates are in general not specified by $n_{i}$ alone; instead they depend on the entire state vector $\mathbf{n}$. The rates in Eqs. (4) and (5) by themselves, therefore, do not describe a well-defined stochastic process for species $i$.

Inserting the definitions from Eq. (3) into the expressions in Eqs. (4) and (5), one finds

[TABLE]

Our analysis in following sections is based on the assumption that we can formulate a closed death-birth process for species $i$, with birth and death rates $T_{i}^{\pm}(n_{i})$, dependent only on $n_{i}$. As we will describe below, an exact description of this form can be obtained for the case of homogeneous rates across opinions (see Sec. IV). If the imitation or mutation rates are heterogeneous this approach constitutes an approximation (see Sec. V).

III.2 Stationary distribution for individual opinions

At long times the $m$-species distribution $P(\mathbf{n},t)$ becomes time-independent; we write $P_{\rm st}(\mathbf{n})$ for the stationary distribution. In principle, this stationary distribution can be obtained from the linear set of equations for $P_{\rm st}$, obtained by setting the right-hand side of the master equation (2) to zero. In practice, it is very difficult to evaluate this analytically. However, as we will describe next, we can obtain the marginal distributions $P_{i}(n_{i})=\sum_{\mathbf{n}_{-i}}P_{\rm st}(\mathbf{n})$ in some cases; the notation $\sum_{\mathbf{n}_{-i}}\dots$ in this expression indicates a sum over all variables $n_{1},\dots,n_{m}$, except $n_{i}$.

We assume that we can formulate a closed death-birth process for species $i$, with birth and death rates $T_{i}^{\pm}(n_{i})$. The stationary distribution of this process, is then the marginal $P_{i}$, and it is given by

[TABLE]

This expression for the marginal stationary distribution can be obtained by standard methods from the backward master equation (see for example Ewens (2004); Traulsen and Hauert (2009)). We will test these predictions against simulations for the homogeneous and heterogeneous noisy voter models in Sections IV and V. In the figures and in Appendix A we express the marginals of the stationary distribution in terms of ${x_{i}=n_{i}/N}$; we then use the notation ${{\cal P}_{i}(x_{i})}$.

III.3 Noise-induced transition

As discussed above, noisy voter models show a transition between parameter ranges with unimodal and multi-modal stationary distributions. This transition is noise induced, and occurs as the system size is varied. As a consequence, it is useful to identify a critical population size $N_{c}$ which separates the two types of outcomes. In the case of only two species ($m=2$) this is particularly straightforward, as there exists a population size for which the stationary distribution becomes flat. This population size has been used to define the transition point in refs. Biancalani et al. (2014); Houchmandzadeh and Vallade (2015); see also Carro et al. (2016).

In the multi-state case ($m>2$) the situation is more complicated; the marginal distributions for the $\{x_{i}\}$ do not assume a flat shape for any population size, as we will discuss below. In order to characterise the transition, we therefore consider the shape of the distribution $P_{i}(n_{i})$ at the left and right boundaries of phase space, i.e., near $n_{i}=0$ and $n_{i}=N$, respectively. Specifically we find the population sizes $N_{L}$ and $N_{R}$ for which $P_{i}(0)=P_{i}(1)$ and $P_{i}(N-1)=P_{i}(N)$, respectively. These can be used to separate the unimodal and multi-modal regimes. We will discuss this in more detail for the homogeneous and heterogeneous multi-state models in Secs. IV and V.

III.4 Switching times between states of consensus

The system is said to have reached consensus on opinion $i$ if $n_{i}=N$, implying $n_{j}=0$ for all $j\neq i$. We refer to this as consensus state $i$. If the population is not in a consensus state, then at least two of the $n_{j}$ are non-zero, i.e., several opinions are represented in the population; we call these situations mixed states. As shown in Fig. 1, the system can transition between consensus states, with intermediate periods in mixed states. We will refer to this as ‘switching’ between consensus states, and we proceed to calculate the typical time between such switches.

In order to define switching times, we first introduce the concept of an ‘arrival’ at a consensus state. We say that an arrival at consensus state $i$ occurs at time $t$ if the system transitions from a mixed state into consensus state $i$ at time $t$, and if the last consensus state the population visited before $t$ was not $i$. This is illustrated in Fig. 2.

We define the mean switching time $\tau$ as the average time that elapses between two subsequent arrivals. By definition, these arrivals are at two different consensus states. This is again illustrated in Fig. 2. The quantity $\tau^{-1}$ is the mean number of arrivals per unit time.

In order to compute $\tau$, we generalized the procedure used in ref. Houchmandzadeh and Vallade (2015). First, we consider the mean time it takes a single species $j$ to reach state $n_{j}=N$ starting from $n_{j}=0$; we denote this time by $t_{j}^{0\to N}$. This is a standard hitting-time problem, and we find Norris (1998); Traulsen and Hauert (2009)

[TABLE]

To proceed, we now assume that the population is currently in consensus state $i$. The mean time it takes any of the species $j\neq i$ to reach the state $n_{j}=N$ is then given by $t_{j}^{0\to N}$. This expression is exact, provided closed-form processes can be formulated for the individual species.

To be able to describe the switching time of the multi-state model, we now treat the processes for each of the species $j\neq i$ as independent. We associate each process with a ‘clock’, and say that the clock for species $j$ ‘ticks’ when $n_{j}=N$. We assume that the rate with which this happens is constant in time for each process, and given by $1/t_{j}^{0\to N}$. This is an approximation. First, due to the coupling of the different species, the processes for the different types $j$ are in general not independent. Second, the distribution of waiting times until state $n_{j}=N$ is reached is not exponential for one-step Markov processes of the type above; instead, the hitting time can be written as the sum of multiple exponential random variables Karlin and McGregor (1959); Keilson (1979).

Proceeding nevertheless with this approximation, the total rate for any of the clocks $j\neq i$ to tick (i.e., any species $j\neq i$ to take over the population) is given by $\sum_{j\neq i}(t_{j}^{0\to N})^{-1}$. As a consequence, we find the mean waiting time until any of the consensus states $j\neq i$ is reached as

[TABLE]

This approximates the average time for the system to reach any consensus state $j\neq i$, if it is currently in consensus state $i$.

The mean time between two subsequent arrivals can then be written as

[TABLE]

where $p_{i}$ denotes the proportion of the number of arrivals at consensus state $i$ among the total number of arrivals. In other words, of all arrivals the system makes, a fraction $p_{i}$ occurs at consensus state $i$. Evidently, one has ${\sum_{i=1}^{m}p_{i}=1}$.

We now proceed to estimate the coefficients $p_{i}$. To do this, we first consider the typical time between two arrivals at the same consensus state $i$. This can be approximated as

[TABLE]

This expression can be understood by first assuming that the population is in consensus state $i$. The term $\tau_{i}$ in Eq. (11) is the mean time that elapses until the system reaches any of the other consensus states, $j\neq i$. At that point there will be no individuals of type $i$ ($n_{i}=0$). The second term, $t_{i}^{0\to N}$, is the average time required to reach consensus state $i$ again.

The quantity $p_{i}$ in turn is proportional to the number of arrivals at $i$ per unit time, given by (\overset{\text{\scriptsize\leftrightarrow}}{\tau_{i}})^{-1}. With appropriate normalization, we find

[TABLE]

To summarise, the overall switching time is obtained as

[TABLE]

using the expressions in Eqs. (8), (9), and (11).

IV Homogeneous model

We first consider the fully homogeneous set-up. If the rates are homogeneous across species, i.e., if $r_{ij}\equiv r$ and $\varepsilon_{ij}\equiv\varepsilon$, the expressions in Eq. (III.1) reduce to

[TABLE]

In particular, the $T_{i}^{\pm}$ are now fully specified by $n_{i}$ alone. Species $i$ undergoes a well-defined death-birth process with rates given by Eq. (IV). All species $j\neq i$, can be lumped together into one ‘other species’ and there are $N-n_{i}$ individuals of this other species. We note that no approximations have been made in deriving this result. As far as the marginal process for any particular species $i$ is concerned, the dynamics effectively reduces to a two-species process, for which analytical descriptions are available (see e.g. Ewens (2004)). We will use these tools to calculate stationary distributions for individual species, as well as switching times between consensus states.

IV.1 Marginals of the stationary distribution

For uniform rates, the marginal probability distributions $P_{i}(n_{i})$ are identical across species and can be obtained by using the rates from Eqs. (IV) in Eq. (7). These closed-form expressions can be evaluated numerically.

Results can be found in Fig. 3. In panel (a) we show the case of $m=2$ species, already discussed in refs. Biancalani et al. (2014); Houchmandzadeh and Vallade (2015). This is shown mainly for comparison.

Given that ${n_{2}=N-n_{1}}$, the model with two species is, by construction, described by a single death-birth process. As seen in the figure, simulation results and theory are in agreement. For $N>N_{c}$ the resulting distribution is unimodal; for $N<N_{c}$ a bimodal shape is observed. At the critical population size $N=N_{c}$ the distribution is flat; this transition point is $N_{c}=r/\varepsilon$ Biancalani et al. (2014); Houchmandzadeh and Vallade (2015). The stationary distribution is symmetric about $x=1/2$ for all $N$.

In panel (b) of Fig. 3 we show the model with $m=5$ opinion states. Again, simulations and theoretical predictions are in agreement. As in the case of two states, a transition from a unimodal to a bimodal marginal distribution is observed as the population size is decreased. Bimodality in the marginal single-species distributions indicates multi-modality in the state space of all $m$ species, as illustrated in Fig. 1(c) for $m=3$. In contrast with the two-species case, the marginal distributions for the abundances of individual species never becomes flat in Fig. 3(b). We also note that the marginal distributions do not exhibit any particular symmetry, despite the fact that the imitation and mutation rates are homogeneous across species.

IV.2 Phase diagram

As mentioned in Sec. III.3, we use the population sizes at which ${P_{i}(0)=P_{i}(1)}$ or ${P_{i}(N-1)=P_{i}(N)}$ to characterise the transition between the unimodal and multi-modal regimes. The condition for the left edge, $P_{i}(0)=P_{i}(1)$, translates into ${T_{i}^{+}(0)=T_{i}^{-}(1)}$; using Eqs. (IV) one finds the physical solution

[TABLE]

where we have written $\widetilde{\varepsilon}_{\rm tot}=(m-1)\varepsilon/r$, and $\widetilde{\varepsilon}=\varepsilon/r$.

Similarly, requiring $P_{i}(N-1)=P_{i}(N)$ at the right edge yields $T_{i}^{+}(N-1)=T^{-}(N)$, from which we obtain

[TABLE]

For $m=2$ this reduces to the result in refs. Biancalani et al. (2014); Houchmandzadeh and Vallade (2015), i.e., ${N_{L}=N_{R}=r/\varepsilon}$.

The resulting expressions for $N_{L}$ and $N_{R}$ only depend on the relative mutation rate $\varepsilon/r$. This is due to the fact that, up to a re-scaling of time, the reaction rates in Eq. (IV) only depend on the ratio of $\varepsilon$ and $r$. We also note that Eqs. (15) and (16) are symmetric with respect to exchanging $\widetilde{\varepsilon}$ and $\widetilde{\varepsilon}_{\rm tot}$. This is a consequence of the same symmetry in the pairs $\{T^{+}(0),T^{-}(N)\}$ and $\{T^{-}(1),T^{+}(N-1)\}$. We also notice that ${N_{L}\approx\widetilde{\varepsilon}^{-1}}$ and ${N_{R}\approx(\widetilde{\varepsilon}_{\rm tot})^{-1}}$ when ${\varepsilon/r\ll 1}$.

The phase diagram resulting from Eqs. (15) and (16) is depicted in Fig. 4. We use $r=1$ in both panels. In panel (a) we show the population sizes $N_{R}$ and $N_{L}$ as a function of the number of species, $m$, at a fixed mutation rate $\varepsilon$. The upper sold line indicates $N_{L}$ and the lower solid line is $N_{R}$. Results from simulations are shown as markers; as seen in the figure, theoretical predictions are in agreement with simulations.

We note that, in Fig. 4(a), $N_{L}$ is nearly flat as a function of the number of species $m$. This is in line with the behaviour of $N_{L}\approx r/\varepsilon$ for small mutation rates, and can be further understood as follows. The shape of the stationary state distribution near $n_{i}=0$ is largely determined by the frequency of mutation events that generate new individuals of type $i$. This mutation rate is given by $N\varepsilon$; if $n_{i}=0$, each of the $N$ individuals in the population mutate into $i$ with rate $\varepsilon$ [see Eq. (IV)]. This is the case independently of the number of species $m$. Therefore, we expect $N_{L}$ to be approximately independent of $m$, as seen in Fig. 4(a).

Near $n_{i}=N$, however, the shape of the distribution is determined by the rate with which individuals of type $i$ mutate into any other species. If $n_{i}=N$, each of the $N$ individuals mutates into any one of the species $j\neq i$ with rate $\varepsilon$; the total rate of mutation events out of species $i$ is then $N\varepsilon(m-1)$ [see again Eq. (IV)]. For fixed $\varepsilon$, this total rate increases as either $N$ or $m-1$ are increased. Therefore, we expect $N_{R}$ to be a decreasing function of $m$, which is confirmed in the phase diagram in Fig. 4(a). This is again consistent with the behaviour $N_{R}\approx r/[(m-1)\varepsilon]$ for small mutation rates.

In panel (b) of Fig. 4 we show the phase diagram for the choice $\varepsilon=10^{-2}/(m-1)$, i.e., for a fixed value of $\widetilde{\varepsilon}_{\rm tot}=10^{-2}$. Recalling that an individual of any species undergoes mutations to any of the other $m-1$ species with rate $\varepsilon$, the total mutation rate for any individual is constant with this choice, irrespective of the number of species $m$. The model shows the same three phases as in Fig. 4(a), but the shape of the phase lines changes. We now find that $N_{L}$ is nearly independent of $m$, whereas $N_{R}$ is a decreasing function of $m$. This is consistent with the interpretation given above. The shape of the distribution near the left edge is largely determined by ${T_{i}^{+}(0)=10^{-2}N/(m-1)}$. An increased value of $m$ now requires a larger value of $N$ to keep this rate constant. The shape of the distribution on the right edge is governed by ${T_{i}^{-}(N)=10^{-2}N}$; this quantity is now independent of the number of species, $m$.

Similar to Houchmandzadeh and Vallade (2015), our approach so far has not involved any continuous approximation of the discrete state space of the population. The analysis of ref. Biancalani et al. (2014), on the other hand, is based on a diffusion approximation to the dynamics of the two-state noisy voter model. We have extended this approach to the multi-state case. The dashed and dotted lines in both panels of Fig. 4 are the predictions for $N_{L}$ and $N_{R}$ obtained from the diffusion approximation. The mathematical details are discussed in the Appendix.

The shape of $P_{i}(n_{i})$ in the different regimes of Fig. 4 can be summarised as follows [see also Fig. 3(b)]: For large populations (${N>N_{L}}$) the single-species distribution is unimodal (light blue shading in Fig. 4), taking its maxi-mum strictly in the interior of the interval ${0<n_{i}<N}$. For intermediate sizes (${N_{R}<N<N_{L}}$, unshaded) the distribution function is monotonously decreasing in $n_{i}$, taking its maximum value at ${n_{i}=0}$. Finally, for small populations (${N<N_{R}}$) the distribution for $n_{i}$ is multi-modal (dark purple shading), with maxima at ${n_{i}=0}$ and ${n_{i}=N}$.

IV.3 Switching times

In the homogeneous case, the times $t_{j}^{0\to N}$ in Eq. (8) are uniform across species; we write $t^{0\to N}$ for their common value. As a consequence, we find $\tau_{i}=t^{0\to N}/(m-1)$ in Eq. (9), and therefore the mean switching time is also given by

[TABLE]

Theoretical predictions and results from simulations are compared in Fig. 5. We show data for different choices of the number of species $m$; as before we use $r=1$. To allow a better comparison, we plot the switching time as a function of $\widetilde{\varepsilon}_{\rm tot}=(m-1)\varepsilon$. This ensures that the total mutation rate any one species experiences is comparable across the different values of $m$.

The diagram indicates that the switching time becomes smallest at intermediate values of the mutation rate. This can be understood as follows. For small mutation rates, the population will reside at consensus states for relatively long times; escape from these states occurs through mutation. These long sojourn times at the consensus state mean that the switching time, i.e., the typical time between arrivals, will also be large. On the other hand, for large mutation rates the population spends most of its time in the interior of phase space. This leads to long periods in which the population does not visit any of the consensus states. As a consequence, the switching time will also be large. The shortest switching times are thus seen in the intermediate range of mutation rates, when neither of these effects dominate.

V Heterogeneous model

We now turn to instances of the model in which the imitation and mutation rates are not homogeneous across species. We limit the discussion to examples in which $r_{ji}=r_{j}$ and $\varepsilon_{ji}=\varepsilon_{j}$. The approach we develop can be extended to cover more general cases; we briefly comment on this in Appendix B.

The quantity $r_{j}$ describes the rate with which members of species $j$ imitate any other species, and $\varepsilon_{j}$ is the rate with which individuals of type $j$ spontaneously change to any other type. In the context of opinion dynamics, large values of $r_{j}$ and $\varepsilon_{j}$ thus describe beliefs that are only weakly held. If $r_{j}$ and $\varepsilon_{j}$ are low, on the other hand, then $j$ describes a strong view.

V.1 Analytical approximation

With the above choice of reaction rates, the transition rates $T_{i}^{\pm}$ in Eq. (III.1) can be written as

[TABLE]

where

[TABLE]

and similarly for $\left\langle{\varepsilon}\right\rangle_{-i}(\mathbf{n})$. These quantities are weighted averages over the imitation and mutation rates of all species other than $i$; $\left\langle{r}\right\rangle_{-i}$ for example represents the mean rate with which opinion $i$ is imitated by any other species. We note that $\left\langle{r}\right\rangle_{-i}(\mathbf{n})$ and $\left\langle{\varepsilon}\right\rangle_{-i}(\mathbf{n})$ depend on the variables $n_{j}$ ($j\neq i$), i.e., on the entire state of the population. As a consequence, no closed process can be formulated for individual species for the heterogeneous model.

We therefore resort to an approximation. We replace all $\{n_{j}\}$ with $\{Nx_{j}^{*}\}$ in Eq. (19); $x_{j}^{*}$ represents the proportion of individuals of species $j$ at the fixed point of the deterministic rate equations, i.e., the dynamics in the limit of infinite populations. These can be obtained from a Kramers–Moyal expansion of the master equation (2), as detailed in Appendix A. We then define

[TABLE]

and similarly for $\varepsilon^{*}_{-i}$. Our approximation consists in replacing the stochastic quantities $\left\langle{r}\right\rangle_{-i}(\mathbf{n})$ and $\left\langle{\varepsilon}\right\rangle_{-i}(\mathbf{n})$ with $r_{-i}^{*}$ and $\varepsilon^{*}_{-i}$, respectively, in Eq. (V.1). In this way we find an approximate, but closed, effective process for species $i$. This can then be used to obtain marginals of the stationary distribution and switching times, following the procedure described in Sections III.2 and III.4.

V.2 Marginals of the stationary distribution

In order to test the accuracy of the approximations described above, we consider an example with $m=5$ species. We choose the imitation rates as ${r_{j}=r[1-\delta+(j-1)\frac{2\delta}{m-1}]}$, for $j=1,\dots,m$, and simi-larly ${\varepsilon_{j}=\varepsilon[1-\delta+(j-1)\frac{2\delta}{m-1}]}$ for the mutation rates. The coefficients $\varepsilon$ and $r$ represent the mean rates (across species) and are model parameters. This choice indicates that the rates $\{r_{j}\}$ and $\{\varepsilon_{j}\}$ are equally spaced in the intervals $[(1-\delta)r,(1+\delta)r]$ and $[(1-\delta)\varepsilon,(1+\delta)\varepsilon]$, respectively. The parameter $\delta$ therefore describes the relative spread of these rates, and quantifies the degree of heterogeneity. We note that species $j=1$ has the lowest rates $\varepsilon_{j}$ and $r_{j}$ among all species (it represents a strong opinion), and species $j=m$ the highest (it describes a weak opinion). We restrict our analysis to relatively small values of $\delta$; this avoids situations in which the population spends most of its time in one single consensus state (the one of the strongest opinion). In such cases, it is difficult to measure stationary distributions and switching times in simulations.

The marginals of the resulting stationary distribution are shown in Fig. 6. In panel (a) the population is relatively small ($N=50$), panel (b) shows an intermediate case ($N=250$), and in panel (c) $N=900$.

As the population size is increased, the marginal distributions change shape. Each marginal is bimodal for very small populations [panel (a)]; in particular, the slope of each curve is negative near $n_{i}=0$, and positive near $n_{i}=N$. For larger populations the gradient near $n_{i}=N$ changes sign when a critical size, $N_{R,i}$, is reached [panel (b)], similar to what was found in the homogeneous model. Notably though the population size at which this happens can differ across the species, as indicated by the subscript $i$ in $N_{R,i}$. If the population size is larger still, the gradient of the marginal distribution near $n_{i}=0$ will also change sign [panel (c)]; this occurs at population sizes $N_{L,i}$, which can again vary across the different species. Mathematically, $N_{L,i}$ and $N_{R,i}$ can be defined using conditions similar to those in Sec. IV.2; $N_{L,i}$ is the population size for which $T_{i}^{+}(0)=T_{i}^{-}(1)$, and $N_{R,i}$ is defined as the population size for which $T_{i}^{+}(N-1)=T_{i}^{-}(N)$.

V.3 Phase diagram

In Fig. 6(a) the population size is such that $N<N_{R,i}$ for all $i$, so that all marginals have bimodal shape. In panel (b) we have $N_{R,i}<N<N_{L,i}$ for all $i$, and in panel (c) the population is sufficiently large so that $N>N_{L,i}$ for all $i$.

The resulting phase diagram for the model with heterogeneous rates is shown in Fig. 7. It exhibits the three different phases outlined above.

For a given value of the mean mutation rate, the marginal distributions for all species are unimodal if $N>N_{L,i}$ for all $i$. On the other hand, if $N<N_{R,i}$ for all $i$ then all marginals are bimodal. It is useful to define $N_{L}^{\rm min}\equiv\mbox{min}_{i}\,N_{L,i}$, and $N_{L}^{\rm max}\equiv\mbox{max}_{i}\,N_{L,i}$, and similarly for $N_{R}^{\rm min}$ and $N_{R}^{\rm max}$. We also define $N_{L}^{\rm hom}$ and $N_{R}^{\rm hom}$ as the corresponding quantities for the homogeneous model, with $r_{i}\equiv r,\varepsilon_{i}\equiv\varepsilon$ for all $i$ [see Eqs. (15) and (16)].

In Fig. 7, we plot these quantities for different choices of $m$. The white unshaded region shows population sizes for which $N_{R}^{\rm min}<N<N_{L}^{\rm max}$. Simulation results for $N_{R}^{\rm min}$, $N_{R}^{\rm max}$, $N_{L}^{\rm min}$, and $N_{L}^{\rm max}$ are also shown, and confirm the theoretical predictions.

V.4 Switching times

In Fig. 8 we demonstrate the accuracy of the analytical approximation in Eq. (13) for the switching times in the model with heterogeneous rates. Panel (a) shows the times $t_{i}^{0\to N}$ for the different species $i$ in the model with $m=5$. Imitation and mutation rates were chosen as described in Sec. V.2. As seen in the figure, the overall agreement is reasonable, although not quantitatively accurate. This is not surprising, as the approximation of Eq. (20) is a relatively severe intervention. Nevertheless, the approximation accurately predicts the qualitative shape of $t_{i}^{0\to N}$ as a function of the mutation rate, as well as the ordering across species. We note that the different curves in panel (a) show species $i=1,\dots,5$ from bottom to top. This is in-line with intuition: individuals who are of opinion $i=5$ are easily convinced of other opinions and are likely to change their views spontaneously. It is therefore difficult for this opinion to spread in the population; as a consequence $t_{5}^{0\to N}$ is large. Individuals who are of opinion $i=1$, on the other hand, do not easily change to any other opinion; this leads to a relatively small time $t_{1}^{0\to N}$.

In Fig. 8(b) we show the resulting switching time $\tau$ for the model with heterogeneous rates. Results are reported for models with different numbers of species, $m$, and as a function of $(m-1)\varepsilon$ as in previous figures. In each case the switching time shows a minimum as a function of the mean mutation rate, as in the homogenous model.

It is interesting to note that the switching time is found to decrease at a fixed mutation rate $(m-1)\varepsilon$ when the number of species increases. To understand this, we observe that the spacing between rates for different species is proportional to $\delta/(m-1)$; see the definitions at the beginning of Sec. V.2. At fixed $\delta$ species $i=1$ is thus the strongest opinion by a large margin when $m$ is small. In the extreme case there are only two opinions; one is a very strong view and the other very weak in comparison. The population visits consensus on the stronger opinion frequently, and remains there for relatively long times. This leads to long switching times. When there are more species available (larger values of $m$) and $\delta$ is kept fixed, the opinion spectrum becomes more finely spaced and the advantage of the most strongly held view is reduced. This leads to a reduction in switching times.

A similar effect occurs when the parameter $\delta$ is varied keeping all other model parameters fixed. This parameter controls the spread of the rates $r_{j}$ and $\varepsilon_{j}$. We find in simulations that switching times increase with increasing spread $\delta$. This is again due to an increased advantage of the most strongly held opinion.

VI Summary and discussion

We have introduced a multi-state noisy voter model and studied its behaviour in a population with all-to-all interaction. It can be seen as a stylised model of imitation and mutation processes in opinion dynamics or evolution, for example. The model shows a noise-driven transition between a state of diversity for large populations and a mostly ordered state in small populations. In the diverse state multiple opinions (or species) are present in the population at most times, and the resulting stationary distribution is unimodal, with a peak in the interior of state space. The ordered regime occurs for small population sizes; the system is frequently found in consensus states, in which all individuals are of the same type. Switches between these consensus states are possible, as the dynamics do not permit absorbing states for non-zero mutation rates.

The transition between the two types of behaviour is a consequence of the balance of mutation effects, driving the system towards mixed states, and intrinsic noise due to the discreteness of the population. Similar to what was reported in ref. Biancalani et al. (2014) for the model with two states, this noise is strongest in the interior of phase space, and drives the system away from mixed states. The overall amplitude of this noise scales as $N^{-1/2}$; as a consequence consensus states are more likely in small populations, when the noise is strong.

While this general behavior is well-known in the two-species noisy voter model, our analysis focuses on the case of multiple species, and specifically includes models with heterogeneity in the properties of the different species. We find that the transition persists, and we devise a method with which to calculate marginals of the multi-species stationary distribution. Our approach is based on formulating a reduced stochastic process for single species. This reduction is exact if imitation and mutation rates are homogeneous across species; in the heterogeneous model it is an approximation. We use the effective processes for individual species to obtain analytical results for marginals of the stationary distribution of the model. This allows us to identify the population sizes $N_{L}$ and $N_{R}$, which characterise the transition from a unimodal to a multi-modal stationary distribution. The effective dynamics for single species can also be used to approximate the typical switching times between different consensus states.

There are some limitations to our work. The most severe restriction is that, in its present form, the approach requires all-to-all interaction; this simplifies the analysis greatly. It is an obvious challenge to extend the results we have obtained to multi-state noisy voter models on lattices and more complex graphs. For example, it would be interesting to see if and how our ideas can be combined with pair-approximation methods used recently for noisy voter models in Peralta et al. (2018b). A further limitation of our approach consists in its inability to capture correlations between species. The method focuses on the marginal dynamics of individual species, and so by construction correlations are not part of the analysis. In order to include these, one next step may be to consider the effective dynamics of pairs of species, while lumping together the remaining species. This would generate an effective three-species model; it is not immediately clear how to proceed with the analytical calculation of joint distributions or switching times. However, progress might be possible in the homogeneous case and exploiting symmetries, similar for example to what has been done in the context of cyclic dynamics in Cremer et al. (2008); Mobilia (2010).

Despite these limitations, we believe our work can contribute to future research in several ways. For example, the method of decoupling multi-species dynamics into several effective single-species processes is likely to be applicable more widely. This includes other models of opinion dynamics and language evolution. In many of these instances evolution is ‘neutral’, i.e., no species has an intrinsic advantage over any other. In such homogeneous models, an approach based on effective single-species processes can be expected to produce accurate results. On the other hand, there are also situations in which heterogeneity is relevant, including evolutionary processes in biology. The method we have put forward can then be useful as an approximation.

Acknowledgements

FHA thanks Consejo Nacional de Ciencia y Tecnología (CONACyT, Mexico) for support.

Appendix A Diffusion approximation

Insight into the model can be gained by approximating the discrete degrees of freedom by continuous variables, and the stochasticity in the dynamics as Gaussian noise. This is based on an expansion in the inverse size of the population, known as the diffusion approximation in mathematical biology Hartl and Clark (2007). In this Appendix we briefly summarise the outcome of this approximation for the multi-state noisy voter model.

A.1 Kramers–Moyal expansion and stochastic differential equation

Formally, the approximation results in a set of coupled stochastic differential equations (SDEs) for the variables $x_{i}=n_{i}/N$. These can be obtained by carrying out a Kramers–Moyal expansion of the master equation (2) in the limit of large, but finite population sizes $N$. Alternatively, these SDEs can also be written down using Kurtz’ theorem Kurtz (1978), or Gillespie’s derivation of the chemical Langevin equation for reaction systems Gillespie (2000).

For the multi-state noisy VM, as defined by Eqs. (2) and (3), one obtains

[TABLE]

where the $\{\xi_{i}(t)\}$ are zero-average Gaussian noise variables, and where

[TABLE]

We have introduced the notation ${\cal T}_{j\to i}(\mathbf{x})\equiv T_{j\to i}(N\mathbf{x})/N$. The $f_{i}(\mathbf{x})$ are commonly referred to as ‘deterministic drift’ Risken (1996). They describe the flow of the dynamics in the limit of infinite populations, $N\to\infty$. In this limit intrinsic stochasticity reduces to zero, due to the pre-factor $N^{-1/2}$ in Eq. (21). The noise variables $\xi_{i}(t)$ are uncorrelated in time, but correlated across components. Writing

[TABLE]

one has, for example based on Kurtz’ theorem Kurtz (1978),

[TABLE]

See also Bladon et al. (2010) for further details in a different context. We note that $\sum_{j}B_{ij}=0$ for all $i$, reflecting the fact that the total number of individuals in the population is constant.

In explicit form, one finds the drift terms for the multi-state noisy VM as follows,

[TABLE]

and the correlation matrix (across components) for the variables $\xi_{i}(t)$ has entries

[TABLE]

for $i\neq j$, and

[TABLE]

These results are valid both in the case of homogeneous and heterogeneous imitation and mutation rates across species.

A.2 Homogeneous rates

We now focus on individual species in the model with homogeneous rates. The Fokker-Planck equation governing the marginal distribution for species $i$ is given by

[TABLE]

where $f_{i}(x_{i})$ and $B_{ii}(x_{i})$ are found to depend only on $x_{i}$ in the homogeneous model. This reduction is similar to the one observed for the transition rates $T_{i}^{\pm}$ in Sec. IV. We note that Eq. (A.2) can be obtained from a direct Kramers–Moyal expansion of the master equation describing the death-birth process for species $i$, with rates as in Eq. (IV).

In explicit form we have

[TABLE]

and

[TABLE]

In the homogeneous model these are the same for all species.

The stationary solution of Eq. (A.2) is found as

[TABLE]

where $B_{ii}^{\prime}(y)$ is the derivative of $B_{ii}(y)$ with respect to $y$, and where $C_{i}$ is a constant ensuring normalisation.

We note that the derivative of the marginal stationary distribution with respect to $x_{i}$ is given by

[TABLE]

The conditions used to define $N_{L}$ and $N_{R}$ translate into ${\cal P}_{i}^{\prime}(0)=0$ and ${\cal P}_{i}^{\prime}(1)=0$, respectively. We therefore find

[TABLE]

within the diffusion approximation.

From Eq. (33) one then obtains

[TABLE]

Results from Eqs. (34) are shown as dashed lines in Fig. 4.

In the limit of small mutation rates, $\varepsilon\ll r$, we find

[TABLE]

This reproduces the asymptotic behaviour of $N_{L}$ and $N_{R}$ in Eqs. (15) and (16).

The results given in this section are for the case of homogeneous rates; one can then formulate a closed Fokker-Planck equation for the marginal distribution of single species [Eq. (A.2)]. If the imitation and mutation rates vary across species, this is no longer possible in exact form. However, an approximate closed Fokker-Planck equation can be formulated for ${\cal P}_{i}$, following the principles of Sec. V.1. The approximation consists of replacing $x_{j}$, $j\neq i$, with the fixed point values $x_{j}^{*}$ of the deterministic dynamics for infinite populations. The analog of Eqs. (33) can then be obtained following the steps described above.

A.3 Further criterion to characterise the transition

In this section we briefly discuss a further method to identify the transition between the regimes with unimodal and multi-modal stationary distributions. This is motivated by the procedure used in Biancalani et al. (2014). We restrict the discussion to models with homogeneous imitation and mutation rates across species.

The solution in Eq. (31) can be written in two equivalent forms,

[TABLE]

where we have dropped the index $i$ for simplicity, and where

[TABLE]

The constant $D$ in Eq. (36) ensures normalisation. We stress that both forms of ${\cal P}$ in Eq. (36) describe the same mathematical expression. The purpose of giving both representations will become clear below. We also note that the only dependence of ${\cal P}(x)$ on $N$ is in $\gamma$, and the only dependence on $x$ is in $q(x)$ and $B(x)$.

In ref. Biancalani et al. (2014) a similar form of the stationary distribution is obtained for the model with two species; our result reduces to this known case when $m=2$. We note that the dependence on $x$ in the numerator of Eq. (36) vanishes for $m=2$, as $\alpha=0$ in this case. In ref. Biancalani et al. (2014) the phase transition for the two-species model was identified as the population size $N$ for which the stationary distribution is flat; for $m=2$, this is equivalent to setting the exponent of $B(x)$ in Eq. (36) to zero, leading to $\gamma=1$, i.e., $N_{c}=r/\varepsilon$.

We propose a similar criterion to characterise the transition for the multi-state model: we find the population size at which the exponent $1-\gamma(1\pm\alpha)$ in the denominator of Eq. (36) vanishes, notwithstanding the fact that a residual dependence on $x$ remains through $q(x)$.

Proceeding on this basis, we use

[TABLE]

to characterise the population size associated with the transition between the unimodal and multi-modal regimes. We note that solving Eq. (38) for $N$ leads to two different solutions, one for each choice of the sign in front of $\alpha$. We label these as $N_{\pm}$, and find

[TABLE]

which for $\varepsilon/r\ll 1$ simplifies to

[TABLE]

To see this note that $A\approx 1$, and $\alpha\approx\frac{m-2}{m}$ in this limit. Eqs. (40) reproduce the asymptotic behaviour of the expressions for $N_{L}$ and $N_{R}$ in Eqs. (34) for $\varepsilon/r\ll 1$. Results from Eqs. (39) are shown as dotted lines in Fig. 4.

Appendix B General heterogeneous model

We briefly address the model with general heterogeneous imitation and mutation rates, $r_{ij}$ and $\varepsilon_{ij}$.

Using the transition rates in Eq. (3), the effective birth and death rates for species $i$ in Eq. (III.1) can be written as

[TABLE]

where

[TABLE]

The averages $\rho_{i}(\mathbf{n})$, $\mu_{i}(\mathbf{n})$ and $\lambda_{i}(\mathbf{n})$ are similar to those in Eq. (19). The only complication is a second index of the object which is being averaged. The object $\nu_{i}$ is a uniform average of $\varepsilon_{ij}$ over all species $j\neq i$, and not dependent on the state $\mathbf{n}$ of the population.

As in the main text, one can proceed by replacing $n_{i}/N$ with the values $x_{i}^{*}$ at the deterministic fixed point. This approximation results in closed definitions for death-birth processes for individual species $i$. These can then be analysed following the same steps as in Sections III.2 and III.4.