Spatial effect on stochastic dynamics of bistable evolutionary games

Let us assume that we have a population composed of N (≫ 1) individuals, n of which follow strategy A and N − n of which follow strategy B. We assume that the population is well-mixed, i.e. all the individuals play games with all the other individuals.

We describe the games performed between the individuals using a payoff matrix

$$\begin{equation} \left( \begin{array}{@{}cc@{}} a & b \\ c & d \end{array} \right), \end{equation} \tag{ 1 }$$

which specifies payoffs of the games in the following manner: The game between two A individuals gives both of them a payoff a. The game between an A individual and a B individual gives them payoffs b and c respectively. The game between two B individuals gives both of them a payoff d. Then, when the number of A individuals is n, the respective average payoffs of A and B individuals are Π_A(n/N) and Π_B(n/N), where

$$\begin{eqnarray} &&\Pi_A(q) := q a + \left( 1 - q \right) b, \qquad \Pi_B(q) := q c + \left( 1 - q \right) d. \end{eqnarray} \tag{ 2 }$$

Here, q := n/N is the proportion of the A individuals in the population (note that we allowed self-interaction for simplicity). In this paper, we consider only a bistable 'coordination game', i.e. we impose the following conditions:

$$\begin{equation} a > c, \ \ b < d. \end{equation} \tag{ 3 }$$

The condition means that, to obtain a higher payoff, it is always best to play the same strategy as the opponent.

As an update rule, we adopt a 'pairwise comparison process' [37, 38]. In this process, two individuals, called a focal individual and a role individual, are selected randomly at a rate λ. The focal individual adopts the strategy of the role individual with probability p(Δ Π) := (1 + w Δ Π)/2, where Δ Π := Π_r − Π_f and Π_f and Π_r are the payoffs of the focal and role individuals respectively. The parameter w controls the effect of the payoffs on the spread of strategies and is hence called 'the intensity of selection'. Thus the strategy yielding a higher payoff is more likely to be imitated by others. We assume w ∈ [0, 1] and 0 < a−c < 1, 0 < d−b < 1 to guarantee p(Δ Π) ∈ [0, 1].

To avoid fixation, we introduce mutations (i.e. the possibility that imitation fails): If the role individual is A (resp. B), the focal individual adopts strategy B (resp. A) with a small probability μ_A (resp. μ_B). We assume that the mutation rates are so small that the dynamics of the expectation value of n remains bistable (see section 2.2).

The parameters used in the present paper are summarized in table 1.

We define the transition rates for the processes, in which the number of A individuals increase or decrease by one, as λ W₊(n/N) and λ W₋(n/N), respectively, where W_± are the dimensionless transition rates. By adopting the update rule described above, W_± are given as follows:

$$\begin{eqnarray} &&\fl W_{+}(q) := (1- \mu_A)q \left(1-q \right)\frac{1 + w \left[ \Pi_A(q) - \Pi_B(q) \right]}{2} + \frac{\mu_B}{2}\left(1-q \right)^2, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} &&\fl W_{-}(q) := (1- \mu_B)q \left(1-q \right)\frac{1 + w \left[ \Pi_B(q) - \Pi_A(q) \right]}{2} + \frac{\mu_A}{2} q^2 . \end{eqnarray} \tag{ 5 }$$

The system is then described by a continuous time Markov process with these transition rates. Note that these transition rates depend on payoff matrix elements only through a − c and d − b because

$$\begin{equation} \Pi_A(q) - \Pi_B(q) = (a-c+d-b)q - (d-b) . \end{equation} \tag{ 6 }$$

The master equation for P(n, t), the probability that the number of A is n at time t, can be written as

$$\begin{eqnarray} &&\fl \partial_t P(n,t) = \lambda W_{+}\left(\frac{n-1}{N}\right)P(n-1,t) + \lambda W_{-}\left(\frac{n+1}{N}\right)P(n+1,t) \nonumber\\ &&- \left[\lambda W_{+}\left(\frac{n}{N}\right) + \lambda W_{-}\left(\frac{n}{N}\right)\right] P(n,t). \end{eqnarray} \tag{ 7 }$$

Let us first consider the dynamics of the expectation value of n, which we denote by 〈n(t)〉. In the limit N → ∞, the stochastic fluctuations can be neglected and 〈n(t)〉 obeys the deterministic equation

$$\begin{equation} \frac{\rm d}{{\rm d}t} \langle n(t) \rangle = \lambda \left [ W_{+}(\langle n(t) \rangle/N) - W_{-}(\langle n(t) \rangle/N) \right]. \end{equation} \tag{ 8 }$$

It is convenient to introduce the new variable q := n/N, which represents the proportion of A in the population and rewrite the equation as, with dimensionless time τ := λ t/N,

$$\begin{eqnarray} &&\frac{\rm d}{{\rm d}\tau} \langle q(\tau) \rangle = W_{+}(\langle q(\tau) \rangle) - W_{-}(\langle q(\tau) \rangle). \end{eqnarray} \tag{ 9 }$$

For the model considered here, W₊(q) − W₋(q) is a cubic function of q:

$$\begin{eqnarray} &&\fl W_{+}(q)-W_{-}(q) = q(1-q)\left\{w \left(1-\frac{\mu_A + \mu_B}{2} \right)\left[(a-b-c+d)q + b-d\right] + \frac{\mu_B - \mu_A}{2}\right\} \nonumber\\ &&+ \frac{\mu_B}{2} (1-q)^2 - \frac{\mu_A}{2}q^2. \end{eqnarray} \tag{ 10 }$$

Figure 1 shows a graph of W₊(q) − W₋(q) and the flow of q determined by (9).

Zoom In Zoom Out Reset image size

The fixed points of the deterministic equation (9) are determined by the condition W₊(q) − W₋(q) = 0, i.e. the transition rates of both processes balance. For coordination games (a > c, b < d) without mutation (μ_A = μ_B = 0), it can be easily shown that the system has three fixed points: two stable fixed points q₁ = 0 (= everyone adopts strategy B) and q₃ = 1 (= everyone adopts strategy A) and one unstable fixed point q₂ = (d − b)/(a − c+d − b) (= mixture of strategy A and B). The positions of the fixed points shift slightly when small mutations are introduced. However, the system remains bistable, i.e. it has two stable fixed points (q₁ and q₃) and one unstable fixed point (q₂) in the range [0, 1] (see figure 1). Hereafter, we assume that the mutation rates μ_A and μ_B are sufficiently small so that the system remains bistable.

As discussed above, the expectation value of n exhibits bistability. In addition, as can be seen from the definition of the model, fixation to one strategy is impossible because of mutations. Thus, although the system remains at states q₁ or q₃ for an extremely long time, large fluctuations can occasionally carry the system from one state to the other. Therefore, states q₁ and q₃ are metastable and have a long but finite lifetime.

To demonstrate this feature, we performed a Monte Carlo simulation of the stochastic process described by the master equation (7) based on the Gillespie algorithm [39]. We show an example of stochastic time evolution of the number of A individuals for N = 50 in figure 2, which shows that while the system stays around the two stable fixed points (q₁ and q₃) for a long time, there are rare transitions between these two metastable states due to stochastic fluctuations. These transitions are expected to be very rare when N is large, because fluctuations around metastable states are suppressed.

Zoom In Zoom Out Reset image size

In this paper, we evaluate the lifetimes of these metastable states, which are defined to be the mean waiting times until the system escapes from the given metastable state and undergoes transition to the other metastable state. In the following section, we mainly focus on the lifetime of the metastable state q₃, because the lifetime of the metastable state q₁ can be calculated in the same manner.

2.3.1. Path integral expression.

The lifetime of the metastable state q₃ is calculated from P_{q3, wm}(T), the probability that the system stays around q₃ from τ = 0 to τ = T (≫ 1) without ever visiting the other metastable state q₁. P_{q3, wm}(T) decays exponentially with T and the inverse of the decay rate gives a lifetime of q₃.

P_{q3, wm}(T) can be approximated by the probability that q(T) = q₃ and q(τ)≠ q₁ (∀τ ∈ [0, T]) given that q (0) = q₃. We express the latter probability using the path integral formulation of stochastic processes [30, 40], in which the probability of given paths can be expressed as a summation, with some weight, over the paths. With this technique, we obtain (see appendix A for the derivation)

$$\begin{eqnarray} &&P_{q_3,{\rm wm}}(T) \simeq \int^{\prime}_{q(0) = q_3, q(T) = q_3} \mathcal{D} q \mathcal{D} p \ {\rm e}^{-N S[q(\cdot),p(\cdot)]} , \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} &&S[q(\cdot),p(\cdot)] := \int_{0}^{T} {\rm d}\tau [ p(\tau) \partial_{\tau} q(\tau) - H_0(q(\tau),p(\tau)) ], \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} &&H_0(q,p) := ( {\rm e}^p - 1 ) W_{+}(q) + ({\rm e}^{-p} - 1)W_{-}(q) , \end{eqnarray} \tag{ 13 }$$

where W_± are given by (4) and (5), p is a new variable conjugate to q and $\int^{\prime}_{q(0) = q_3, q(T) = q_3} \mathcal{D}q \mathcal{D}p$ represents the (restricted) summation over all the paths (q, p) satisfying q (0) = q(T) = q₃ and q(τ)≠ q₁ (∀τ ∈ [0, T]).

The meaning of this expression is that the probability density of each path q (·), p (·) is exp(−NS[q (·), p (·)]) and that the desired probability is a weighted sum over all paths satisfying the conditions. Note that the present formalism is an analog of the path integral in quantum mechanics. From this analogy, S and H₀ are called an action and a Hamiltonian, respectively. The path integral expression is useful for the analysis of N ≫ 1 cases as discussed in the next section.

2.3.2. Semiclassical approximation.

We will now evaluate the path integral expression (11) under the assumption N ≫ 1. We adopt a semiclassical approximation, in which the path integral is approximated by contribution from the stationary paths (q, p) of action S (i.e. δS[q, p] = 0) and the fluctuations around them because, when N is large, only the stationary paths contribute dominantly to the path integral. This approximation is an extension of the steepest descent method in evaluating integrals and an analog of the WKB (or semiclassical) approximation in quantum mechanics. It is known that the semiclassical approximation is suitable for treating rare events [30].

The stationary condition δS = 0 yields differential equations:

$$\begin{eqnarray} &&\frac{{\rm d}q}{{\rm d}\tau} = \frac{\partial H_0}{\partial p} = {\rm e}^p \ W_{+}(q) - {\rm e}^{-p} \ W_{-}(q), \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} &&\frac{{\rm d}p}{{\rm d}\tau} = -\frac{\partial H_0}{\partial q} =- ({\rm e}^p -1)W^{\prime}_{+}(q) - ({\rm e}^{-p}-1)W^{\prime}_{-}(q), \end{eqnarray} \tag{ 15 }$$

(where $W^{\prime}_{\pm}(q) := {\rm d}W_{\pm}(q)/{\rm d}q)$ subjected to the boundary conditions q (0) = q(T) = q₃ and the condition q(τ)≠ q₁ (∀τ ∈ [0, T]). Because the equations have the same form as Hamilton's equations of motion in analytical mechanics, the solutions (q, p) satisfying these differential equations are called 'classical' trajectories in the following discussion. Figure 3 shows phase portraits of the flow of these equations. Note that H₀(q, p) is conserved along each trajectory.

Zoom In Zoom Out Reset image size

There are two important types of trajectories. The first type are the trajectories on the horizontal line of p = 0. These trajectories represent the dynamics of the expectation value of q, because the equation of motion (14) for p = 0 coincides with the equation of 〈q〉 (9). Note that H₀(q, p) = 0 for these trajectories.

The second type are the trajectories shown by bold lines in the left panel of figure 3, connecting stable fixed points ((q₁, 0) and (q₃, 0)) and an unstable fixed point ((q₂, 0)). These two trajectories are called 'activation trajectories' [29, 34]. Because the Hamiltonian is always constant on the connected trajectories, H₀(q, p) = 0 holds on the activation trajectories. Hence, the shape of the activation trajectories can be collectively expressed as p = p_a(q), where

$$\begin{equation} p_a(q) := - \log \left[ \frac{W_{+}(q)}{W_{-}(q)} \right], \end{equation} \tag{ 16 }$$

which can be obtained from the condition H(q, p_a(q)) = 0 (see (13)).

In the semiclassical approximation of P_{q3, wm}, classical trajectories should be properly chosen so that the conditions q (0) = q(T) = q₃ and q(τ)≠ q₁ (∀τ ∈ [0, T]) are satisfied. Therefore, the classical trajectories relevant for the evaluation of P_{q3, wm} are restricted to (i) a trivial solution (q(τ), p(τ)) ≡(q₃, 0) and (ii) nontrivial solutions circulating on the closed trajectory through (q₂, 0) and (q₃, 0) (the red closed trajectory shown in the right panel of figure 3), i.e. the closed trajectory composed of an activation trajectory and a p = 0 trajectory. Note that only these trajectories spend a significant amount of time going around and can satisfy the boundary conditions q (0) = q(T) = q₃ for T ≫ 1. These solutions are called 'bounce solutions' [30].

The action of the trivial solution is zero, whereas the action of nontrivial solutions are given by n S_{q3, wm}, where n (= 1, 2, 3, ...) is the number of the rotation and S_{q3, wm} is the action per one cyclic motion on the closed trajectory described above. Using the zero-energy condition H₀(q, p) = 0 and the expression for the activation trajectories p_a(q), S_{q3, wm} can be expressed as

$$\begin{eqnarray} &&S_{q_3, {\rm wm}} = \int_{q_3}^{q_2} p_a(q) {\rm d}q , \end{eqnarray} \tag{ 17 }$$

which is equal to the area shaded in red in the right panel of figure 3.

Summing up all the contributions from the bounce solutions, we obtain

$$\begin{eqnarray} &&{P_{{q_3},{\rm{wm}}}}(T) \simeq \sum\limits_{n = 0}^\infty {\int_0^{{t_2}} {\rm{d}} } {t_1} \cdots \int_0^T {\rm{d}} {t_n}{\left( { - K{{\rm{e}}^{ - N{S_{{q_3},{\rm{wm}}}}}}} \right)^n}\\ &&= \exp \left( { - K{{\rm{e}}^{ - N{S_{{q_3},{\rm{wm}}}}}}T} \right). \end{eqnarray} \tag{ 18 }$$

Here, the prefactor K is a positive value determined by Gaussian integrals around bounce solutions.

Thus, the lifetime of the metastable state q₃, τ_{q3, wm}, can be expressed as

$$\begin{equation} \tau_{q_3, {\rm wm}} \simeq \frac{1}{K} \exp(N S_{q_3, {\rm wm}}). \end{equation} \tag{ 19 }$$

The lifetime of the other metastable state q₁ can be evaluated in the same manner, except that the other closed trajectory (the blue trajectory in the right panel of figure 3) should be used:

$$\begin{eqnarray} &&\tau_{q_1, {\rm wm}} \simeq \frac{1}{K^{\prime}} \exp\left( N S_{q_1, {\rm wm}}\right), \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} &&S_{q_1, {\rm wm}} := \int_{q_1}^{q_2} p_a(q) {\rm d}q , \end{eqnarray} \tag{ 21 }$$

where the prefactor K' is a positive quantity (determined by Gaussian integrals around bounce solutions around q₁).

It can be seen from (19) and (20) that the lifetimes are extremely long (τ_{qi, wm}≫ 1) under the condition N ≫ 1 assumed in this paper. The lifetimes strongly depend on the action S_{qi, wm} because of the presence of a large factor N in the exponent. In this paper, we focus on the exponent of the lifetimes and ignore the weak parameter dependence of the prefactors K and K'.

In this section, we briefly summarize the parameter dependence of the lifetime τ_{qi, wm}(i ∈ {1, 3}).

2.4.1. N dependence.

2.4.2. w dependence.

Figure 4 shows S_{q1, wm} and S_{q3, wm} as the functions of w (the intensity of selection). As can be seen from the figure, the action grows with w. This result indicates that, when N ≫ 1, the lifetimes τ_{qi, wm}∝exp(N S_{qi, wm}) (i ∈ {1, 3}) increase rapidly as natural selection becomes strong.

Zoom In Zoom Out Reset image size

2.4.3. Payoff matrix dependence.

As can be seen from (4)–(6), the stochastic dynamics of the evolutionary game considered here depends on payoff matrix elements only through a − c and d − b. The dependence of S_{q1, wm} and S_{q3, wm} on these two parameters is shown in figure 5 (only the region 0.2 < a−c < 0.8 and 0.2 < d−b < 0.8 is plotted because outside this region there are parameter sets for which the system is not bistable).

Zoom In Zoom Out Reset image size

This dependence can be intuitively explained in the following manner: The larger a − c is, the more advantageous A individuals are in a population dominated by A, compared with B individuals in the same population, and the longer the lifetime of q₃. Almost the same discussion applies to d − b and the lifetime of q₁. Note that when μ_A = μ_B, a symmetry relation S_{q1, wm}(a − c, d − b) = S_{q3, wm}(d − b, a − c) holds.

We consider a 1D array of M patches. Let l be the separation between neighbouring patches and $\tilde{L} := Ml$ be the length of the system. We impose a periodic boundary condition. In each patch, there are N_p (≫ 1) individuals, who take either strategy A or strategy B and change their strategies according to the same evolutionary game as described in section 2. We further assume that the individuals can move between neighbouring patches. This migration process is modelled as the 'swapping' of individuals so that the number of individuals per patch is conserved (see figure 6).

Zoom In Zoom Out Reset image size

Let n_i be the number of A individuals in the ith patch (i ∈ {1, 2, ···, M}). The parameters of the evolutionary rule in each patch are the same as those described in section 2 (see table 1). The transition rates of n_i due to the strategy update can be written as λ W₊(n_i/N_p) and λ W₋(n_i/N_p), where W_± are given in (4) and (5), respectively.

The migration (swapping) processes are defined as follows: A pair of neighbouring patches, i and j, is randomly chosen at a rate σ. One individual is randomly picked up from each patch and they are swapped. The rate of the process (n_i, n_j) → (n_i −1, n_j +1) due to swapping is given by σ W_m(n_j/N_p, n_j/N_p), where

$$\begin{equation} W_{\rm m}(q_i,q_j) := q_i (1-q_j). \end{equation} \tag{ 22 }$$

With these transition rates, the master equation for the probability distribution P (n, t) on the population configuration n = (n₁, n₂, ···, n_M) is written as

$$\begin{eqnarray} &&\fl\partial_t P({\bit n}, t) = \sum_{i=1}^{M} \left\{\lambda W_{+}\left( \frac{n_i-1}{N_p} \right)P({\bit n} - {\bit e}_i, t) + \lambda W_{-}\left(\frac{n_i+1}{N_p}\right)P({\bit n} + {\bit e}_i, t) \right. \nonumber\\ &&\left. - \left[ \lambda W_{+}\left(\frac{n_i}{N_p}\right) + \lambda W_{-}\left(\frac{n_i}{N_p}\right) \right] P({\bit n}, t) \right\} + \sum_{\langle i, j \rangle} \left[ \sigma W_{\rm m} \left(\frac{n_i+1}{N_p}, \frac{n_j-1}{N_p}\right)\right. \nonumber\\ &&\left. \times P({\bit n} + {\bit e}_i - {\bit e}_j, t) \right. \left. - \sigma W_{\rm m}\left(\frac{n_i}{N_p}, \frac{n_j}{N_p} \right) P({\bit n},t) \right], \end{eqnarray} \tag{ 23 }$$

where 〈i, j〉 indicates neighbouring patches and

$$\begin{equation} {\bit e}_i := (0, \cdots, 0, \stackrel{i}{\breve{1}}, 0, \cdots, 0). \end{equation} \tag{ 24 }$$

In this section, we derive a differential equation for the expectation values of n_i and examine its steady solutions and their stability, which play an important role in understanding the transitions between metastable states in the spatial model.

3.2.1. Deterministic equations.

Let 〈n_i〉 be the expectation value of n_i. Then, 〈n_i〉(i ∈ {1, 2, ···, M}) obey

$$\begin{eqnarray} &&\fl \frac{{\rm d} \langle n_i \rangle}{{\rm d}t} = \lambda W_{+}(\langle n_i\rangle / N_p) - \lambda W_{-}(\langle n_i \rangle / N_p) \nonumber\\ &&+ \, \sigma W_{\rm m}(\langle n_{i+1} \rangle /N_p, \langle n_i \rangle /N_p) + \sigma W_{\rm m}(\langle n_{i-1}\rangle /N_p, \langle n_i \rangle /N_p) \nonumber\\ &&-\, \sigma W_{\rm m}(\langle n_i \rangle /N_p, \langle n_{i+1}\rangle /N_p ) - \sigma W_{\rm m}(\langle n_i \rangle /N_p, \langle n_{i-1}\rangle /N_p ). \end{eqnarray} \tag{ 25 }$$

By introducing new variables q_i := n_i/N_p (the proportion of A individuals in the ith patch) and a rescaled time τ := λ t/N_p, we obtain

$$\begin{eqnarray} \frac{{\rm d} \langle q_i \rangle}{{\rm d}\tau} &=& W_{+}(\langle q_i \rangle) - W_{-}(\langle q_i \rangle) + \frac{\sigma}{\lambda} \left( \langle q_{i+1} \rangle + \langle q_{i-1} \rangle - 2 \langle q_{i} \rangle \right) . \end{eqnarray} \tag{ 26 }$$

In this paper, we assume σ≫λ, i.e. migration processes occur sufficiently faster than strategy update processes. Under this assumption, 〈q_i〉 changes smoothly as a function of i and can be expressed by a function of a continuum spatial degree of freedom x:

$$\begin{equation} q(x,\tau) := \langle q_{\frac{x}{l}}(\tau) \rangle \qquad (x \in [0, \tilde{L}]). \end{equation} \tag{ 27 }$$

The second term on the right hand side of (26) is approximated by D∂² q/∂x², where D := σ l²/λ is a diffusion constant. By introducing the dimensionless space variable $\xi := x/ \sqrt{D}$, we arrive at the following reaction diffusion equation:

$$\begin{eqnarray} &&\frac{\partial q(\xi,\tau)}{\partial \tau} = \frac{\partial ^2 q(\xi,\tau)}{\partial {\xi}^2} + W_{+}(q(\xi,\tau)) - W_{-}(q(\xi,\tau)) \qquad \left(\xi \in [0,L] \right), \end{eqnarray} \tag{ 28 }$$

where $L := \tilde{L}/\sqrt{D}$ is the rescaled system size.

3.2.2. Steady solutions.

During the analysis of the well-mixed model discussed in section 2, fixed points and their stability played important roles. In the spatial model, steady solutions will play similar roles.

Let q_s(·) be a steady solution of (28). q_s obeys an ordinary differential equation

$$\begin{equation} \frac{{\rm d}^2 q_s(\xi)}{{\rm d} {\xi}^2} + W_{+}(q_s(\xi)) - W_{-}(q_s(\xi)) = 0. \end{equation} \tag{ 29 }$$

Integrating this equation yields

$$\begin{eqnarray} &&\frac{1}{2} \left( \frac{{\rm d} q_s(\xi)}{{\rm d}\xi} \right) ^2 + V(q_s(\xi)) = {\rm const.} =:E , \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} &&V(q) := \int_{0}^{q} \left[ W_{+}(q^{\prime}) - W_{-}(q^{\prime}) \right] {\rm d}q^{\prime}. \end{eqnarray} \tag{ 31 }$$

By regarding ξ as the 'time', we can describe q_s as a coordinate of a particle moving in one dimension, subjected to potential V. Hence, the trajectories of q_s on the (q, dq/dξ) plane can be expressed by the contours of 'energy' E defined in (30) (dotted lines in the left panel of figure 7). Because of the periodic boundary condition q_s(ξ+L) = q_s(ξ), q_s must be a closed orbit or a point on the (q, dq/dξ) plane. Therefore, the steady solutions must be

• (a)  
  uniform solutions: q_s(ξ) ≡ q_i(i ∈ {1, 2, 3}) (recall that q_i are the solutions of W₊(q) − W₋(q) = 0 defined in section 2.2) or
• (b)  
  non-uniform periodic solutions: closed orbits surrounding q₂.

Zoom In Zoom Out Reset image size

First, we consider the uniform solutions, which always exist regardless of system size L. Because $W^{\prime}_{+}(q) - W^{\prime}_{-}(q)$ is negative for q₁, q₃ and positive for q₂ (see figure 1), two solutions, q_s(ξ) ≡ q₁ and q_s(ξ) ≡ q₃, are linearly stable and q_s(ξ) ≡ q₂ is linearly unstable; the former solution corresponds to the two metastable states, whereas the latter solution corresponds to the 'marginal' state located at the boundary q₂, below which the system goes to q₁ and above which the system goes to q₃.

Next, we examine the properties of the non-uniform solutions. The forms of these solutions depend on the system size L, because the period of the solution must coincide with L:

$$\begin{equation} 2\int_{q_{\rm min}}^{q_{\rm max}} \frac{{\rm d}q}{\sqrt{2(E-V(q))}} =L. \end{equation} \tag{ 32 }$$

Here, q_min and q_max are the minimum and the maximum of q in the trajectory, which is determined by E. The left panel of figure 7 shows the phase portrait for the case of V(q₁) > V(q₃), where the non-uniform solution in the limit of L → ∞ corresponds to the homoclinic orbit starting from q₃ and ending at q₃ (the blue line in the figure). The right panel of figure 7 shows the spatial profile of a non-uniform solution for the same parameter set as the left panel. As can be seen from the figure, the solution represents a 'nucleus' of B individuals surrounded by a region dominated by A individuals. Note that for the opposite case (i.e. V(q₁) < V(q₃)), the form of the non-uniform solution becomes 'upside down' compared to the right panel of figure 7, i.e. it represents a nucleus of A individuals surrounded by a region dominated by B individuals.

It can be shown that the obtained non-uniform solutions are unstable, i.e. even an infinitesimally small perturbation drives the system away from the non-uniform solutions and toward the stable states q₁ or q₃ (see appendix B). Because of such 'critical' behaviour, these non-uniform steady solutions are called 'critical nuclei' [36] and we denote them by q_cn.

As L decreases, both q_min and q_max approach q₂ and at a critical length L_c, they coalesce with q₂. This indicates that critical nuclei do not exist when L is smaller than L_c, which is calculated, by linearizing (29) around q₂, as

$$\begin{equation} L_c = \frac{2\pi}{\sqrt{W^{\prime}_{+}(q_2)- W^{\prime}_{-}(q_2)}}. \end{equation} \tag{ 33 }$$

L_c gives the characteristic length scale of the critical nuclei and plays an important role in considering the spatial effect on the transitions between metastable states, as discussed later in section 4.

As shown in the previous section, the spatial model considered deterministically has two stable steady solutions, i.e. the uniform solutions q_s(ξ) ≡ q₁ and q_s(ξ) ≡ q₃, which we henceforth simply call q₁ and q₃, respectively. If stochasticity is taken into account, these solutions correspond to metastable states: although the system stays at q₁ or q₃ for an extremely long time, it occasionally undergoes transitions from one state to the other as a result of large stochastic fluctuations. In this section, we evaluate the lifetimes of these metastable states. We show the calculation of the lifetime of q₃ (the lifetime of q₁ can be calculated in the same manner).

3.3.1. Path integral expression.

To evaluate the lifetime of a metastable state q₃, we calculate P_{q3, sp}(T), the probability that the system stays around q₃ from τ = 0 to τ = T (≫ 1) without ever visiting the other metastable state q₁. P_{q3, sp}(T) decays exponentially with T and the inverse of the decay rate gives the lifetime of q₃.

The path integral formalism, which was applied to the well-mixed model in section 2, is also applicable to the spatial model (see appendix A for the detail). P_{q3, sp} can be approximated by the probability that q (·, T) = q₃ and q (·, τ)≠ q₁ (∀τ ∈ [0, T]) given that q (·, 0) = q₃. By adopting a continuum description assuming σ≫λ and using a dimensionless spatial coordinate $\xi := jl/\sqrt{D} \ \ ( j \in \left\{1,2, \cdots, M \right\} )$ and the rescaled system size $L := \tilde{L}/\sqrt{D}$, we obtain

$$\begin{eqnarray} &&\fl P_{q_3, {\rm sp}}(T) \simeq \int^{\prime}_{q(\cdot, 0) = q_3, q(\cdot, T) = q_3} \mathcal{D} q\mathcal{D}p \exp \left( -N_p \sqrt{\frac{\sigma}{\lambda}} S[q(\cdot,\cdot), p(\cdot,\cdot)] \right) , \end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray} &&\fl S[q(\cdot,\cdot),p(\cdot,\cdot)] := \int_{0}^{T} {\rm d}\tau \int_{0}^{L} {\rm d}\xi \left\{p(\xi,\tau) \partial_{\tau}q(\xi,\tau) - h [q(\cdot,\cdot),p(\cdot,\cdot)](\xi,\tau) \right\} , \end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray} &&\fl h[q(\cdot,\cdot),p(\cdot,\cdot)](\xi,\tau) := h_0(q(\xi,\tau),p(\xi,\tau)) \nonumber\\ &&- \left\{[\partial_\xi q(\xi,\tau)][\partial_\xi p(\xi,\tau)] -q(\xi,\tau)[1-q(\xi,\tau)][\partial_\xi p(\xi,\tau)]^2 \right\} , \end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray} &&h_0(q,p) := ({\rm e}^p -1)W_{+}(q) + ({\rm e}^{-p}-1)W_{-}(q) , \end{eqnarray} \tag{ 37 }$$

where the prime in the path integral indicates the restriction to paths that satisfy q (·, τ)≠ q₁ (∀τ ∈ [0, T]).

3.3.2. Semiclassical approximation.

The semiclasslcal approximation can be applied in almost the same manner as was performed in section 2.3.2 for the well-mixed model. We consider stationary solutions of the action represented by (35). The stationary condition δS = 0 leads to partial differential equations

$$\begin{eqnarray} \fl \partial_\tau q &=& {\rm e}^{p}W_{+}(q) - {\rm e}^{-p}W_{-}(q) + \left[ \partial^{2}_{\xi}q - 2 q(1-q)\partial^{2}_{\xi} p - 2(1-2q)(\partial_\xi p)(\partial_\xi q) \right] , \end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray} \fl \partial_\tau p &=& -({\rm e}^{p}-1)W^{\prime}_{+}(q) - ({\rm e}^{-p}-1)W^{\prime}_{-}(q) - \left[ \partial^{2}_{\xi} p + (1-2q)(\partial_\xi p)^2 \right], \end{eqnarray} \tag{ 39 }$$

subjected to a boundary condition q (·, 0) = q (·, T) = q₃. Note that when p = 0, the equation for q coincides with (28), i.e. the equation for the expectation value.

We now calculate P_{q3, sp}(T) (T ≫ 1) by considering solutions to the equations of motion (38) and (39), under the boundary condition q(ξ, 0) = q (·, T) = q₃ and the restriction that q never visits q₁. Unlike the well-mixed model, it is not easy to obtain such solutions in the spatial model. However, we can infer qualitative features of these solutions from the analysis of the well-mixed model, where the solution is composed of two trajectories: a trajectory moving from q₃ to q₂ in the p < 0 region and a trajectory moving from q₂ to q₃ in the p = 0 region (see section 2.3.2, figure 3 right panel). Note that the unstable fixed point q₂ plays the role of a 'watershed', i.e. a dividing point between the two stable states q₁ and q₃. In the spatial model, there are two kinds of unstable solution: (a) the uniform solution q(ξ) ≡ q₂ and (b) the critical nucleus. Therefore, it is inferred that there are two kinds of bounce solution:

• (a)  
  a solution which is initially q(ξ) ≡ q₃, then changes into q(ξ) ≡ q₂ and finally returns to q(ξ) ≡ q₃ (see figure 8 left panel), and
• (b)  
  a solution which is initially q(ξ) ≡ q₃, then changes into q = q_cn and finally returns to q(ξ) ≡ q₃ (see figure 8 right panel).

Zoom In Zoom Out Reset image size

We call the former a 'uniform bounce solution' (or solution α) and the latter a 'non-uniform bounce solution' (or solution β). It is clear that the non-uniform bounce solution exists only when L > L_c.

Let $S_{\alpha, q_3}$ and $S_{\beta, q_3}$ be the action of the bounce solutions α and β calculated from (35)–(37). If L < L_c, α is the only possible bounce solution. Hence, P_{q3, sp}(T) is calculated to be

$$\begin{equation} P_{q_3, {\rm sp}}(T) \simeq \exp \left( -K_{\alpha, q_3} {\rm e}^{-N_p\sqrt{\frac{\sigma}{\lambda}}S_{\alpha, q_3}} \ T \right), \end{equation} \tag{ 40 }$$

where $K_{\alpha,q_3}$ is a prefactor determined by the Gaussian integrals around the stationary solution α. Then, τ_{q3, sp}, the lifetime of the metastable state q₃ can be expressed as $\tau_{q_3, {\rm sp}} \simeq K_{\alpha, q_3}^{-1} \exp\left( N_p \sqrt{\sigma/\lambda} S_{\alpha, q_3} \right)$. If L > L_c, there are bounce solutions α and β. By summing up all the contributions from these two bounce solutions, we obtain

$$\begin{eqnarray} &&\fl P_{q_3, {\rm sp}}(T) \simeq \sum_{n=0}^{\infty} \frac{1}{n!} \sum_{k=0}^{n} {}_{n}C_k \int_{0}^{T}{\rm d}t_1 \cdots \int_{0}^{T} {\rm d}t_n \quad \times \left( -K_{\alpha, q_3} {\rm e}^{-N_p \sqrt{\frac{\sigma}{\lambda}}S_{\alpha, q_3}} \right)^k \left( -K_{\beta, q_3} {\rm e}^{-N_p \sqrt{\frac{\sigma}{\lambda}}S_{\beta, q_3}} \right)^{n-k} \end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&\nonumber = \exp\left[ -\left( K_{\alpha, q_3} {\rm e}^{-N_p \sqrt{\frac{\sigma}{\lambda}} S_{\alpha, q_3}} + K_{\beta, q_3} {\rm e}^{-N_p \sqrt{\frac{\sigma}{\lambda}}S_{\beta, q_3}} \right) T \right], \end{eqnarray} \tag{ 42 }$$

where $K_{\alpha, q_3}$ and $K_{\beta, q_3}$ are prefactors determined by the Gaussian integrals around α and β respectively. Then,

$$\begin{equation} \frac{1}{\tau_{q_3, {\rm sp}}} \simeq K_{\alpha, q_3} {\rm e}^{-N_p \sqrt{\frac{\sigma}{\lambda}} S_{\alpha, q_3}} + K_{\beta, q_3} {\rm e}^{-N_p \sqrt{\frac{\sigma}{\lambda}}S_{\beta, q_3}}. \end{equation} \tag{ 43 }$$

Under the condition $N_p\sqrt{\sigma/\lambda} \gg 1$ assumed in this paper, the lifetime is well-approximated by one of the two terms, which has the smaller value of the action. Let γ be either α or β, which corresponds to the smaller value of action (e.g. if $S_{\alpha, q_3} > S_{\beta, q_3}$, then γ = β. In fact, for the case treated in the next section, we obtain γ = β). Then, we obtain $\tau_{q_3, {\rm sp}} \simeq K_{\gamma, q_3}^{-1} \exp \left( N_p \sqrt{\sigma/\lambda} S_{\gamma, q_3} \right)$. Thus, the expression for the lifetime is summarized as

$$\begin{equation} \tau_{q_3, {\rm sp}} \simeq \left\{\begin{array}{@{}ll} \displaystyle\frac{1}{K_{\alpha, q_3}} \exp\left( N_p \sqrt{\frac{\sigma}{\lambda}} S_{\alpha, q_3} \right) & (\textrm{for}L < L_c) \\ \displaystyle\frac{1}{K_{\gamma, q_3}} \exp\left( N_p \sqrt{\frac{\sigma}{\lambda}} S_{\gamma, q_3} \right) & (\textrm{for}L > L_c) \\ \end{array}\right. . \end{equation} \tag{ 44 }$$

To discuss the L dependence of the lifetimes, we numerically calculated the action of the bounce solutions (see appendix C for the detail). If two kinds of bounce solution exist (i.e. if L > L_c), we choose the one with the smaller action. The action obtained is denoted by S_{qi, sp}(i ∈ {1, 3}).

Figure 9 shows the L dependence of the action calculated in this manner for a specific parameter set, for which V(q₁) > V(q₃) holds. If L < L_c, both S_{q1, sp}(L) and S_{q3, sp}(L) increase linearly with the system size L. A different feature appears when the system size exceeds the threshold length L_c, where the non-uniform bounce solutions become important; S_{q3, sp} is now independent of L, whereas S_{q1, sp} continues increasing with L. A detailed discussion of these results and their intuitive meaning is provided in the next section.

Zoom In Zoom Out Reset image size

First, we summarize the results for the well-mixed model (section 2). For the evolutionary game dynamics considered here (coordination games with mutations), there are two metastable states denoted by q₁ and q₃, where q ∈ [0, 1] is the proportion of A individuals in the system; the state q₁ corresponds to the population dominated by strategy B and q₃ corresponds to the one dominated by strategy A. Although the system spends an extremely long time at q₁ or q₃, it can undergo a transition from one state to the other, via an unstable state q₂ (satisfying q₁ < q₂ < q₃). The lifetime of the metastable states q₁ and q₃ are calculated to be

$$\begin{equation} \tau_{q_i, {\rm wm}} \sim \exp(N S_{q_i, {\rm wm}}) \qquad (i \in \left\{1,3 \right\}), \end{equation} \tag{ 45 }$$

where S_{q1, wm} and S_{q3, wm} are the actions of the bounce solution defined by (21) and (17), respectively. Because S_{q1, wm} and S_{q3, wm} are positive and independent of N, the lifetimes τ_{qi, wm} increase exponentially with the population size N.

Next, we summarize the results for the spatial model (section 3). The lifetimes of the metastable states q₁ and q₃ in the spatial model are evaluated to be

$$\begin{equation} \tau_{q_i, {\rm sp}} \sim \exp\left(N_p \sqrt{\frac{\sigma}{\lambda}} S_{q_i, {\rm sp}}(L) \right) \ \ (i \in \left\{1,3 \right\}), \end{equation} \tag{ 46 }$$

where S_{q1, sp} and S_{q3, sp} are the positive functions of the system size $L = M\sqrt{\lambda/\sigma}$ and N_p is the population number per patch. Note that L is proportional to the total population number N_pM, where M is the number of patches. As discussed in section 3.2.2, the spatial model has a characteristic length scale L_c. For small systems (L < L_c), S_{q3, sp} grows linearly with L, whereas for large systems (L > L_c), it becomes independent of L (see the right panel of figure 9). Thus, from (46) one can conclude that the lifetime τ_{q3, sp} grows exponentially with L for L < L_c, while it is independent of L for L > L_c. In other words, while the model is effectively well-mixed when L < L_c, it behaves differently from the well-mixed model when L > L_c. We stress that this difference is quite drastic due to the large factor $N_p \sqrt{\sigma / \lambda} \ (\gg 1)$ appearing in (46). On the other hand, S_{q1, sp} has a linear dependence on L for both L < L_c and L > L_c (see the left panel of figure 9). This implies that there is no qualitative difference in the system size dependence of τ_{q1, sp} between the two regions L < L_c and L > L_c, which is in clear contrast to what is observed for τ_{q3, sp}.

Note that the difference between the lifetimes of q₁ and q₃ originates from that of V(q₁) and V(q₃), where V is defined in (31). If V(q₁) < V(q₃) (the opposite to the above situation), the roles of q₁ and q₃ are reversed so that the lifetime of q₁ takes a constant value for L > L_c while the lifetime of q₃ keeps growing (the intuitive reason for this behaviour is provided in the next subsection). The condition V(q₁) > V(q₃) (resp. V(q₁) < V(q₃)) can be rewritten as A₁ > A₃ (resp. A₁ < A₃), where

$$\begin{eqnarray} &&{A_1}{\rm{ }}: = - \int_{{q_1}}^{{q_2}} {[{W_ + }(q) - {W_ - }(q)]{\rm{d}}q > 0,} \end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray} &&{A_3}: = \int_{{q_2}}^{{q_3}} {[{W_ + }(q) - {W_ - }(q)]{\rm{d}}q > 0,} \end{eqnarray} \tag{ 48 }$$

Note that A₁ and A₃ correspond to the two areas enclosed by the curve W₊(q) − W₋(q) and the horizontal axis in figure 1. Without mutation (μ_A = μ_B = 0), the condition A₁ > A₃ (resp. A₁ < A₃) is clearly equivalent to the condition q₂ > 1/2 (resp. q₂ < 1/2), where q₂ = (d − b)/(a − c+d − b), because W₊(q) − W₋(q) is a cubic function of q. Thus, the condition V(q₁) > V(q₃) (resp. V(q₁) < V(q₃)) is equivalent to d − b > a − c (resp. d − b < a − c). With a small mutation, although the precise condition is slightly modified, the above discussion still holds approximately.

The observation above shows that in the limit of small mutation, the lifetimes of the metastable states are simply characterized by the relative size of a − c and d − b in the original game (1). The value a − c (> 0) measures the coordination advantage of playing A over playing B when matched with an A individual. Similarly, the value d − b (> 0) measures the coordination advantage of playing B over playing A when matched with a B individual. The strategy with the larger (resp. smaller) coordination advantage is called a 'risk-dominant' (resp. 'risk-dominated') strategy in game theory [4] (e.g. if d − b > a − c, strategy B is risk-dominant and strategy A is risk-dominated). Therefore, our results are summarized as follows: the lifetime of the metastable state of the risk-dominant strategy keeps growing exponentially with the system size, while that of the risk-dominated one saturates after the system size reaches a certain threshold. Hence (when the system is sufficiently large) spatial structure drastically accelerates the selection of the strategy with a larger coordination advantage (= risk-dominant strategy), facilitates its invasion against a risk-dominated strategy and gives a prediction of how natural selection resolves coordination problems in biology.

A similar spatial effect was already found by Ellison for an evolutionary game model (Theorem 3 in [22]). However, a direct comparison with his result is not appropriate because his model assumes one individual per each patch, while our model assumes N_p (≫ 1) individuals per patch.

Here we intuitively discuss the origin of the qualitative differences in the system size dependence of the lifetimes between the well-mixed model and the spatial model. We first point out that during the transition from one metastable state to the other, the system must cross a 'dividing point', which is a marginal unstable steady solution of the system located between two metastable states. Once the system has crossed this dividing point, it evolves with a very high probability along a 'deterministic' path (the trajectory of the expectation value satisfying (9) or (28)), which quickly carries the system to the other metastable state. This view implies that the lifetime of a metastable state is mostly determined by the difficulty of crossing a dividing point from a given metastable state. Then, the system size dependence of the lifetimes can be intuitively explained as follows (again, we restrict our consideration to the case V(q₁) > V(q₃)):

• In the well-mixed model or the spatial model for L < L_c, the unstable state q₂ plays the role of a dividing point. To make a transition from one metastable state to the other, the entire system has to cross the dividing point q₂ (figure 10(a)). As the system size increases, the probability that the system changes from a given metastable state to the dividing point decays exponentially with N (or L), resulting in the exponential dependence of the lifetimes on N (or L).
• In the spatial model for L > L_c, the critical nucleus also plays the role of a dividing point. This indicates that the transition from one metastable state to the other can be caused by 'nucleation', i.e. transition via the critical nucleus. The details of the nucleation process depend on the direction of transition.
  • For a transition from q₃ to q₁ to occur, it is sufficient for the system to change from q₃ to the critical nucleus, which means that the change of only a part of the system, whose length scale is approximately L_c, is sufficient (see figure 10(b) and the spatial form of the critical nucleus in the right panel of figure 7). The difficulty of this change does not depend on L and hence neither does the difficulty of transition from q₃ to q₁. Therefore, the lifetime of q₃ does not grow with L.
  • For a transition from q₁ to q₃ to occur, it is necessary for the system to change from q₁ to the critical nucleus, which means that almost the entire system has to change (figure 10(c)). This change becomes more difficult as L grows and so does the transition from q₁ to q₃. Hence, the lifetime of q₁ grows exponentially with L.

Zoom In Zoom Out Reset image size

It should be noted that the existence of two kinds of transition paths (shown in figure 10) is directly related to two kinds of bounce solution, i.e. the uniform bounce solution and the non-uniform bounce solution, discussed in section 3.3.2. We stress that the essential change of the transition process is described by the change in the bounce solutions.

Note that when V(q₁) < V(q₃), the behaviour of the lifetime of q₁ and q₃ is reversed. This is because the form of the critical nucleus becomes 'upside down' compared with the right panel of figure 7, so that the roles of strategy A (red) and B (blue) in figure 10 should be replaced by each other.

To confirm the intuition described above, we performed a Monte Carlo simulation of the spatial model, specified by the master equation (23) employing the Gillespie algorithm [39]. Figure 11 shows results of the simulation. For a small system (the upper panel, L < L_c), the transition from one metastable state (the red region) to the other state (the blue region) occurs almost uniformly. For a large system (the lower panel, L > L_c), however, the transition proceeds via nucleation; transition occurs first in a small region and then it spreads over the whole system.

Zoom In Zoom Out Reset image size

We point out that the intuition presented above suggests the generality of our results; although we have demonstrated a spatial effect on bistable evolutionary games with a specific update rule, the qualitative behaviour of the lifetimes is expected to be the same for models with other update rules as long as the underlying evolutionary game is bistable. This is because our results rely on the presence of a critical nucleus (which has a characteristic length scale L_c), which in turn originates from the structure of the phase portrait (see figures 1 and 7). It can easily be shown that models with other update rules (e.g. a pairwise comparison process with other forms of p(Δ Π) or the Moran process) also generate phase portraits with the same structure as figures 1 and 7, provided that the payoff matrix elements satisfy (3) and mutation rates are sufficiently small.

The characteristic length L_c in the spatial model is an important quantity that determines the population size dependence of the lifetimes. In this section, we discuss its parameter dependence in detail. The characteristic length L_c is calculated in section 3.2.2 and is given by

$$\begin{equation} L_c = \frac{2\pi}{\sqrt{W^{\prime}_{+}(q_2)- W^{\prime}_{-}(q_2)}}. \end{equation} \tag{ 49 }$$

To observe the parameter dependence more clearly, we derive an approximate expression for L_c: because q₁, q₂ and q₃ are the solutions of the cubic equation W₊(q) − W₊(q) = 0, we obtain an alternative expression for W₊(q) − W₋(q) (see (10)):

$$\begin{equation} \fl W_{+}(q)- W_{-}(q) = - w\left( 1- \frac{\mu_A + \mu_B}{2} \right) (a-b-c+d) (q-q_1)(q-q_2)(q-q_3). \end{equation} \tag{ 50 }$$

By using (49), we obtain

$$\begin{equation} L_c = 2\pi \left/\sqrt{w \left(1-\frac{\mu_A + \mu_B}{2}\right)(a-b-c+d)(q_3-q_2)(q_2-q_1)} \right. . \end{equation} \tag{ 51 }$$

Assuming that μ_A and μ_B are sufficiently small, we arrive at a concise approximate form

$$\begin{equation} L_c \simeq 2\pi \sqrt{\frac{1}{w} \left( \frac{1}{a-c} + \frac{1}{d-b} \right)}, \end{equation} \tag{ 52 }$$

from which we can conclude that L_c decreases with w, a − c and d − b.

Figure 12 shows the w dependence of L_c, where the circles and solid line denote the results obtained from a numerical evaluation of (49) and the approximate expression (52), respectively. One can see that the present approximation works well as w increases.

Zoom In Zoom Out Reset image size

Figure 13 shows the payoff matrix dependence of L_c; the left and right panels are obtained from the numerical evaluation of (49) and the approximate expression (52), respectively. One can see that L_c decreases as a − c and d − b increase and that the approximate expression (52) captures the qualitative feature well.

Zoom In Zoom Out Reset image size

For a comparison with real systems, it is important to judge whether the system under consideration can be regarded as well-mixed or spatial. In this section, we show the criterion in terms of migration rate σ and examine the effect of σ on the lifetimes of the metastable states.

As was pointed out in section 4.1, the system is effectively well-mixed if L < L_c. Because L is defined by $L := \tilde{L}/\sqrt{D} = M \sqrt{\lambda/\sigma}$ (M: the number of patches), the condition L < L_c is equivalent to

$$\begin{equation} M < \sqrt{\frac{\sigma}{\lambda}} L_c. \end{equation} \tag{ 53 }$$

Note that L_c depends only on the game theoretical parameters given in table 1 but not on σ or λ (see (49)). Therefore, in terms of the migration rate σ, the system is effectively well-mixed if

$$\begin{equation} \sigma/ \lambda > M^2 / L_{c}^{2}, \end{equation} \tag{ 54 }$$

i.e. if the ratio between migration rate and update rate is larger than the value determined by M and the game theoretical parameters. On the other hand, when σ is small so that (54) is not satisfied, the system should be described as the spatial model. Then, the transition between the metastable states is affected by the nucleation process described in section 4.2.

Finally, we discuss the σ dependence of the lifetimes of the metastable states. We focus on the exponents of the lifetimes

$$\begin{equation} N_p \sqrt{\frac{\sigma}{\lambda}} \ S_{q_i, {\rm sp}}\left(M \sqrt{\frac{\lambda}{\sigma}}\right) \end{equation} \tag{ 55 }$$

(see (46)). The lifetimes take constant values when σ is large so that (54) is satisfied, i.e. when the system is effectively well-mixed. This is because S_{qi, sp}(L) grows linearly with L when L < L_c and therefore factors depending on σ are cancelled out (see (55)). When σ becomes smaller than the critical value in (54), the lifetimes decrease as σ decreases. Thus, it can be concluded that spatial structure accelerates the transition between metastable states, leading to a drastic reduction of the lifetimes compared to those of the well-mixed model.

When we study coordination games, it is important to determine which metastable state has a longer lifetime, because it gives a criterion for which one of A and B is better off [17]. We briefly discuss the spatial effect on this problem.

For the well-mixed model, the metastable state with the larger value of action S_{qi, wm} has a longer lifetime as can be seen from (45). For the spatial model with a sufficient size (L > L_c), on the other hand, the metastable state with the larger value of potential V(q_i) has a longer lifetime (see section 3.4 and 4.1).

In the limit of small mutation rates (μ_A, μ_B → 0), it can easily be shown that these two conditions coincide: V(q₁) > V(q₃) and S_{q1, wm} > S_{q3, wm} are, respectively, equivalent to the condition d − b > a − c, suggesting that the risk-dominant strategy has a longer lifetime, which is consistent with previous studies [17, 22]. However, in the presence of asymmetric mutation (μ_A, μ_B > 0, μ_A≠μ_B), the conditions are not necessarily equivalent. This observation indicates that the long–lived metastable state may switch from one state to the other as the system size increases. A detailed discussion of this problem is left for a future study.