Healthy and General

Eco-evolutionary model on spatial graphs reveals how habitat

14 min read

Eco-evolutionary model on spatial graphs

We establish an individual-based model (IBM) where individuals are structured over a trait space and a graph representing a landscape. For the sake of simplicity, we consider the case of asexual reproduction and haploid genetics29. Individuals die, reproduce, mutate and migrate in a stochastic fashion, which together results in macroscopic properties. The formulation of the stochastic IBM allows an analytical description of the dynamics at the population level, which links emergent properties to the elementary processes that generate them.

The trait space \(\mathcalX\subseteq \mathbbR^d\) is continuous and can be split into a neutral trait space \(\mathcalU\) and an adaptive trait space \(\mathcalS\). We refer to neutral traits \(u\in \mathcalU\) as traits that are not under selection, in contrast to adaptive traits \(s\in \mathcalS\), which experience selection. The graph denoted by G is composed of a set of vertices v1,v2,…,vM that correspond to habitat patches (suitable geographical areas), and a set of edges that constrain the movement of individuals between the habitat patches. We use the original measure of genetic differentiation for quantitative traits QST (standing for Q-statistics) in the case of haploid populations45,46. We denote the neutral trait value of the kth individual on vi as \(u_k^(i)\), the number of individuals on vi as N(i), the mean neutral trait on vi as \(\overlineu^(i)\), and the mean neutral trait in the metapopulation as \(\overlineu\). It follows that we quantify neutral differentiation QST,u as

$$Q_ST,u=\sigma _B,u^2/(\sigma _B,u^2+\sigma _W,u^2)$$


where \(\sigma _B,u^2=\mathbbE[\frac1M\sum _i\left(\overlineu^(i)-\overlineu\right)^2]\) denotes the expected neutral trait variance between the vertices and \(\sigma _W,u^2=\frac1M\mathop\sum \nolimits_i^M\mathbbE\left[\frac1N^(i)\sum _k\left(u_k^(i)-\overlineu^(i)\right)^2\right]\) denotes the average expected neutral trait variance within vertices. We similarly quantify adaptive differentiation QST,s.

Following the Gillespie update rule47, individuals with trait \(x_k\in {{{{{\mathcalX}}}}}\) on vertex vi are randomly selected to give birth at rate b(i)(xk) and die at rate d(N(i)) = N(i)/K, where K is the local carrying capacity. The definition of d therefore captures competition, which is proportional to the number of individuals on a vertex and does not depend on the individuals’ traits (we relax this assumption later on). The offspring resulting from a birth event inherits the parental traits, which can independently be affected by mutations with probability μ. A mutated trait differs from the parental trait by a random change that follows a normal distribution with variance \(\sigma _\mu ^2\) (corresponding to the continuum of alleles model48). The offspring can further migrate to neighbouring vertices by executing a simple random walk on G with probability m. A schematic overview of the two different settings considered is provided in Fig. 1. Under the setting with no selection, individuals are only characterised by neutral traits so that \({{{{{{{\mathcalX}}}}}}}={{{{{{{\mathcalU}}}}}}}\). For individuals on a vertex with trait xk ≡ uk we define b(i)(xk) ≡ b, so that the birth rate is constant. This ensures that neutral traits do not provide any selective advantage. Under the setting with heterogeneous selection, each vertex of the graph vi is labelled by a habitat type with environmental condition Θi that specifies the optimal adaptive trait value on vi. It follows that, for individuals with traits \(x_k=(u_k,s_k)\in {{{{{{{\mathcalU}}}}}}}\times {{\mathcalS}}\) on vi, we define

$$b^(i)(x_k)\equiv b^(i)(s_k)=b(1-p(s_k-\Theta _i)^2)$$


where p is the selection strength41. This ensures that the maximum birth rate on vi is attained for sk = Θi, which results in a differential advantage that acts as an evolutionary stabilising force. In the following we consider two habitat types denoted by I and II with symmetric environmental conditions θI and θII, so that ΘiθI, θII and θII = − θI = θ, where θ can be viewed as the habitat heterogeneity41.

Fig. 1: Graphical representation of the structure of individuals in the eco-evolutionary model.
figure 1

a Setting with no selection, where individuals are characterised by a set of neutral traits \(u\in {{{{{{{\mathcalU}}}}}}}\). The scatter plots represent a projection of the first two components of u for the individuals present on the designated vertices at time t = 1000, obtained from one simulation of the IBM. b Setting with heterogeneous selection. In this setting, individuals are additionally characterised by adaptive traits \(s\in {{{{{{{\mathcalS}}}}}}}\). Blue vertices favour the optimal adaptive trait value θI, while red vertices favour θII. The scatter plots represent a projection of the first component of u and s for the individuals present on the designated vertices at time t = 1000, obtained from one simulation. The majority of individuals are locally well-adapted and have an adaptive trait close to the optimal value, but some maladaptive individuals originating from neighbouring vertices are also present. m = 0.05.

Deterministic approximation of the population dynamics under no selection

The model can be formulated as a measure-valued point process (30 and Supplementary Note). Under this formalism, we demonstrate in the Supplementary Note how the population size and the trait dynamics show a deterministic behaviour when a stabilising force dampens the stochastic fluctuations. This makes it possible to express the dynamics of the macroscopic properties with deterministic differential equations, connecting emergent patterns to the processes that generate them. In particular, in the setting of no selection, competition stabilises the population size fluctuations, and the dynamics can be considered deterministic and expressed as

$$\partial _tN_t^(i)=N_t^(i)\left[b(1-m)-\fracN_t^(i)K\right]+mb\mathop\sum\limits_j\ne i\fraca_i,jd_jN_t^(j)$$


where \(A=(a_i,j)_1\le i,j\le M\) is the adjacency matrix of the graph G and D = (d1,d2,…,dM) is a vector containing the degree of each vertex (number of edges incident to the vertex). The first term on the right-hand side corresponds to logistic growth, which accounts for birth and death events of non-migrating individuals. The second term captures the gains due to migrations, which depend on the graph topology. Assuming that all vertices with the same degree have an equivalent position on the graph, corresponding to a “mean field” approach (see Methods), one can obtain a closed-form solution from Eq. (3) (see Eq. (12)), which shows that the average population size \(\overlineN\) scales with \(\langle \sqrtk\rangle ^2/\langle k\rangle\), where 〈k〉 is the average vertex degree and \(\langle \sqrtk\rangle\) is the average square-rooted vertex degree. The quantity \(\langle \sqrtk\rangle ^2/\langle k\rangle\), denoted as hd, relates to the homogeneity in vertex degree of the graph and can therefore be viewed as a measure negatively associated with heterogeneity in connectivity. Simulations of the IBM illustrate that hd can explain differences in population size for complex graph topologies with varying migration regimes (Fig. 2a for graphs with M = 7 vertices and Supplementary Fig. 1a for M = 9). This analytical result is connected to theoretical work on reaction-diffusion processes49 and highlights that irregular graphs (graphs whose vertices do not have the same degree) result in unbalanced migration fluxes that affect the ecological balance between births and deaths. Highly connected vertices present an oversaturated carrying capacity (N(i) > bK, see Methods), increasing local competition and lowering total population size compared with regular graphs (Fig. 2a). Because populations with small sizes experience more drift (31 and Supplementary Fig. 2), this result indicates that graph topology affects neutral differentiation not only through population isolation, but also by affecting population dynamics.

Fig. 2: Effect of <l> and hd on average population size \(\overlineN\) and neutral differentiation QST,u in the setting with no selection.
figure 2

a Response of \(\overlineN\) to homogeneity in degree \(h_d=\langle \sqrtk\rangle ^2/\langle k\rangle\) for all undirected connected graphs with M = 7 vertices and m = 0.5. b Response of QST,u to average path length <l> for similar simulations obtained with m = 0.01. c Response of QST,u to homogeneity in degree hd for the same data. In a, b, and c, each dot represents average results from 5 replicate simulations of the IBM, the colour scale corresponds to the proportion of the graphs with similar x and y-axis values (graph density), and the blue line corresponds to a linear fit. d Standardized effect of hd and <l> on QST,u, obtained from multivariate regression models independently fitted on similar data obtained for m = 0.01 and m = 0.5. The contributions of <l> and hd to QST,u are alike for low migration regimes. Error bars show 95% confidence intervals. Analogous results on graphs with M = 9 vertices are presented in Supplementary Fig. 1 and all regression details can be found in Supplementary Table 2.

Nonetheless, the stochasticity of the processes at the individual level can propagate to the population level and substantially affect the macroscopic properties. In particular, neutral differentiation emerges from the stochastic fluctuations of the populations’ neutral trait distribution. These fluctuations complicate an analytical underpinning of the dynamics, and in this case simulations of IBM offer a straightforward approach to evaluate the level of neutral differentiation.

Effect of graph topology on neutral differentiation under no selection

We study a setting with no selection and investigate the effect of the graph topology on neutral differentiation. When migration is limited, individuals’ traits are coherent on each vertex but stochastic drift at the population level generates neutral differentiation between the vertices. Migration attenuates neutral differentiation because it has a correlative effect on local trait distributions. Following21,22,26, we expect that the intensity of the correlative effect depends on the average path length of the graph 〈l〉, defined as the average shortest path between all pairs of vertices50. For a constant number of vertices, 〈l〉 is strictly related to the mean betweenness centrality and quantifies the graph connectivity50. High 〈l〉 implies low connectivity and greater isolation of populations, and hence we expect that graphs with high 〈l〉 are associated with high differentiation levels. We consider various graphs with an identical number of vertices and run simulations of the IBM to obtain the neutral differentiation level QST,u attained after a time long enough to discard transient dynamics (see Methods). We then interpret the discrepancies in QST,u across the simulations by relating them to the underlying graph topologies.

We observe strong differences in QST,u across graphs for varying m, and find that 〈l〉 explains at least 55% of the variation in QST,u across all graphs with M = 7 vertices for (Fig. 2b). Nonetheless, some specific graphs, such as the star graph, present higher levels of QST,u than expected by their average path length. To explain this discrepancy, we explore the effect of homogeneity in vertex degree hd, as we showed in Eq. (12) that it decreases population size, which should in turn increase QST,u by intensifying stochastic drift. We find that hd explains 57% of the variation for low m (Fig. 2c). However, the fit remains similar after correcting for differences in population size (see Supplementary Table 1), indicating that irregular graphs structurally amplify the isolation of populations. Unbalanced migration fluxes lead central vertices to host more individuals than allowed by their carrying capacity. This causes increased competition that results in a higher death rate, so that migrants have a lower probability of further spreading their trait. Highly connected vertices therefore behave as bottlenecks, increasing the isolation of peripheral vertices and consequently amplifying QST,u.

We then evaluate the concurrent effect of 〈l〉 and hd on QST,u with a multivariate regression model that we fit independently for low and high migration regimes (Fig. 2d). The multivariate regression model explains at least 70% of the variation in QST,u for the migration regimes considered and for graphs with M = 7 vertices (see Supplementary Table 2 for details). Moreover, we find that 〈l〉 and hd have akin contributions to neutral differentiation for low m, but the effect of 〈l〉 increases for higher migration regimes while the effect of hd decreases. To ensure that these conclusions can be generalised to larger graphs, we conduct the same analysis on a subset of graphs with M = 9 vertices and find congruent results (Supplementary Fig. 1). In the absence of selection and with competitive interactions, graphs with a high average path length 〈l〉 and low homogeneity in vertex degree hd, or similarly graphs with low connectivity and high heterogeneity in connectivity, show high levels of neutral differentiation.

Deterministic approximation of the population dynamics and adaptation under heterogeneous selection

We next consider heterogeneous selection and investigate the response of adaptive differentiation to the spatial distribution of habitat types, denoted as the Θ-spatial distribution. Adaptive differentiation emerges from local adaptation, but migration destabilises adaptation as a result of the influx of maladaptive migrants. We expect that higher connectivity between vertices of similar habitat type increases the level of adaptive differentiation, because it increases the proportion of well-adapted migrants. Local adaptation can be investigated by approximating the stochastic dynamics of the trait distribution with a deterministic partial differential equation (PDE). We demonstrate under mean-field assumption how the deterministic approximation can be reduced to an equivalent two-habitat model. We analyse the reduced model with the theory of adaptive dynamics36,41 and find a critical migration threshold m that determines local adaptation. m depends on a quantity coined the habitat assortativity rΘ, and we demonstrate with numerical simulations that rΘ determines the overall adaptive differentiation level QST,s reached at steady state in the deterministic approximation.

Heterogeneous selection, captured by the dependence of the birth rate on Θi, generates a stabilising force that dampens the stochastic fluctuations of the adaptive trait distribution. The dynamics of the adaptive trait distribution consequently shows a deterministic behavior and we demonstrate in the Supplementary Note and Supplementary Figs. 3 and 4 that the number of individuals on vi with traits \(s\in \Omega \subset {{{{{{{\mathcalS}}}}}}}\) can be approximated by the quantity ∫Ωn(i)(s)ds, where n(i) is a continuous function solution of the PDE

$$\partial _tn_t^(i)(s)= \, n_t^(i)(s)\left[b^(i)(s)(1-m)-\frac1K\int_{{{{{{{{\mathcalS}}}}}}}}n_t^(i)(\bfs)d{{{\bfs}}}\right]\\ +m\mathop\sum\limits_j\ne ib_j(s)\fraca_i,jd_jn_t^(j)(s)+\frac12\mu \sigma _\mu ^2\Delta _s\left[b^(i)(s)n_t^(i)(s)\right]$$


Equation (4) is similar to Eq. (3), except that it incorporates an additional term corresponding to mutation processes and that the birth rate is trait-dependent. We show how Eq. (4) can be reduced to an equivalent two-habitat model under mean-field assumption. The mean-field approach differs slightly from the setting with no selection because vertices are labelled with Θi. Here we assume that vertices with similar habitat types have an equivalent position on the graph (see Supplementary Fig. 5 for a graphical representation), so that all vertices with habitat type I are characterised by the identical adaptive trait distribution that we denote by \(\overlinen^\bfI\), and are associated with the birth rate \(b^\bfI(s)=b(1-p(s-\theta _\bfI)^2)\). Let P(I, II) denote the proportion of edges connecting a vertex vi of type II to a vertex vj of type I, and let P(I) denote the proportion of vertices vi of type I. By further assuming that habitats are homogeneously distributed on the graph so that \(P(\bfI)=P(\bfII)=\frac12\), Eq. (4) transforms into

$$\partial _t\overlinen_t^\bfI(s)= \,\overlinen_t^\bfI(s)\left[b^{{{{\bfI}}}}(s)(1-m)-\frac1K\int_{{{{{{{{\mathcalS}}}}}}}}\overlinen_t^{{{{{{{{\bfI}}}}}}}}({{{{{{{\bfs}}}}}}})d{{{{{{{\bfs}}}}}}}\right]+\frac12\mu \sigma _\mu ^2(\Delta _sb^{{{{{{{{\bfI}}}}}}}}\overlinen_t^{{{{{{{{\bfI}}}}}}}})(s)\\ +\fracm2\,[(1-r_\Theta )b^{\bfII}(s)\overlinen_t^{{{{{{{{\bfII}}}}}}}}(s)+(1+r_\Theta )b^{{{{{{{{\bfI}}}}}}}}(s)\overlinen_t^{{{{{{{{\bfI}}}}}}}}(t)]$$


(see Methods), where we define

$$r_\Theta =2\left(P({{{{{{{\bfI}}}}}}},{{{{{{{\bfI}}}}}}})-P({{{{{{{\bfI}}}}}}},{{{{{{{\bfII}}}}}}})\right)$$


as the habitat assortativity of the graph, which ranges from −1 to 1. When rΘ = − 1, all edges connect dissimilar habitat types (disassortative graph), while as rΘ tends towards 1 the graph is composed of two clusters of vertices with identical habitat types (assortative graph). Eq. (5) can be analysed with the theory of adaptive dynamics36,38,41, a mathematical framework that provides analytical insights by assuming a “trait substitution process”. Following this assumption, the mutation term in Eq. (5) is omitted and the phenotypic distribution results in a collection of discrete individual types that are gradually replaced by others until evolutionary stability is reached (see Methods and36,38,41 for details). By applying the theory of adaptive dynamics, we find a critical migration rate m

$$m^\star =\frac1(1-r_\Theta )\frac4p\theta ^2(1+3p\theta ^2)$$


so that when m > m, a single type of individual exists with adaptive trait \(s^* =\left(\theta _{{{{{{{{\bfII}}}}}}}}+\theta _{{{{{{{{\bfI}}}}}}}}\right)/2=0\) in the steady-state (see Methods for the derivation of Eq. (7)). In this case, adaptive differentiation QST,s is nil and the average population size is given by \(\overlineN=bK(1-p\theta )^2\). In contrast, when m = 0 and/or rΘ = 1, all individuals are locally well-adapted with trait Θi on vi, and it follows that the average population size is higher and equal to \(\overlineN=bK\), while adaptive differentiation is maximal and equal to \(Q_ST,s={\rmVar}(\Theta )/\left({{{{{{{\rmVar}}}}}}}(\Theta )+0\right)=1\). When 0 < m < m, the coexistence of two types of individuals on each vertex vi is predicted but the calculation of the trait values is more subtle. To understand the effect of m and rΘ on the local trait distributions and on QST,s, we therefore leave behind the adaptive dynamics framework and numerically solve Eq. (5) by including the mutation term. When 0 < m < m, the local trait distributions are bimodal with peaks corresponding to the two types of individuals predicted by the adaptive dynamics. The highest peak corresponds to the well-adapted individuals, whose adaptation is destabilised by the influx of maladaptive migrants (Fig. 3a). This phenomenon is dampened as rΘ increases, since the proportion of maladaptive migrants is reduced in assortative graphs (Fig. 3b). As a consequence, the habitat assortativity rΘ increases the differentiation QST,s when 0 < m < m (Fig. 3c). The simulations further confirm that the adaptive dynamics prediction given by Eq. (7) is still valid when the continuous accumulation of mutations is considered, so that for m > m the local trait distributions obtained from Eq. (5) are unimodal and QST,s vanishes (Fig. 3a,c). Our analysis of the mean-field deterministic approximation Eq. (5) therefore demonstrates that assortative graphs present high levels of adaptive differentiation QST,s. On the other hand, the analysis shows that QST,s rapidly declines with increasing m on disassortative graphs, until QST,s vanishes when m > m.

Fig. 3: Effect of habitat assortativity rΘ and migration m on the local adaptive trait distribution \(\overlinen^{{{{{{{{\bfI}}}}}}}}\) and on the adaptive differentiation level QST,s under the mean field, deterministic approximation Eq. (5).
figure 3

a Effect of m and rΘ on \(\overlinen^{{{{{{{{\bfI}}}}}}}}\). Migration induces the apparition of maladaptive individuals (centred around θII = 0.5), which destabilise local adaptation by displacing the mean value of the well-adapted individuals (centred around θI = − 0.5). Together with the decrease in local adaptation, migration causes a displacement of the mean value of the local trait distribution (represented by the vertical dashed lines), which decreases local population size and adaptive differentiation QST,s. b Similar data for higher rΘ. Increasing rΘ increases population size and QST,s. c Effect of rΘ on QST,s. The red line indicates the critical migration threshold m predicted by Eq. (7); QST,s vanishes when m > m.

Effect of graph topology on adaptive differentiation under heterogeneous selection

To generalise the conclusions drawn from the mean-field deterministic approximation Eq. (5), we generate different Θ-spatial distributions for varying graph topology, and compare outputs of the IBM simulations with those of Eq. (5) (see Methods for the details of the simulations). For each combination of Θ-spatial distribution and graph, we compute the habitat assortativity rΘ, since rΘ can be generalised from Eq. (6) to any graph topology following the original definition of51 as

$$r_\Theta =\frac{{{{{\rmCov}}}}(\Theta _\times ,\Theta _\wedge )}{\sigma _{\Theta _\times }\sigma _{\Theta _\wedge }}$$


where Θ× and Θ denote the sets of habitats found at the toe and tip of each directed vertex of graph V, and 〈Θ×〉, 〈Θ〉 and \(\sigma _{{{{\Theta }}}_\times },\sigma _{{{{\Theta }}}_\wedge }\) denote their respective means and standard deviations (see Supplementary Note). The mean-field deterministic approximation Eq. (5) is in very good agreement with the IBM simulations for general graph ensembles at low migration regimes, and captures the response of \(\overlineN\) and QST,s to rΘ (Fig. 4). Nonetheless, under high migration regimes, higher levels of QST,s are observed in the stochastic simulations compared with the mean field deterministic approximation (Supplementary Fig. 6). We hypothesize that this reinforcement is generated by stochastic drift, which must become the main driver of differentiation when local adaptation is lost for m > m, and perform a multivariate regression analysis to investigate the additional effect of 〈l〉 and hd on QST,s. As expected, the analysis highlights that the effect of 〈l〉 and hd are substantial and complement the effect of rθ for high m (Fig. 5c for graphs with M = 7 vertices and Supplementary Fig. 7a for M = 9), further explaining the discrepancies observed (see Supplementary Table 3).

Fig. 4: Effect of habitat heterogeneity rΘ on QST,s and average population size \(\overlineN\) for general graph ensembles.
figure 4

a Effect of rΘ on QST,s for all undirected connected graphs with M = 7 vertices and varying rΘ, for m = 0.1. b Effect of rΘ on average population size \(\overlineN\) for the same simulations. In a and b, each dot represents average results from 5 replicate simulations of the IBM, the colour scale corresponds to the proportion of the graphs with similar x and y axis values (graph density), and the blue lines correspond to results obtained from the mean-field approximation Eq. (5). Insights from Eq. (5) are congruent with the IBM simulations for complex habitat connectivity patterns at low m. Similar results with m = 0.5 are presented in Supplementary Fig. 6.

Fig. 5: Effect of rΘ, 〈l〉, and hd on QST,s and QST,u in the setting with heterogeneous selection.
figure 5

a Comparison of the response of QST,u to migration with the response of QST,u in the setting with no selection for the complete graph. The dashed vertical blue line corresponds to the critical migration regime m predicted by Eq. (7). Heterogeneous selection increases QST,u when m < m, but local adaptation is lost when m > m, and in this case QST,u reaches similar levels as QST,u in the setting with no selection. b Response of QST,u to rΘ and migration for the path graph. rΘ correlates positively with QST,u for high m, but correlates negatively for low m. In a, b, each plain dot represents average results from 5 replicate simulations, the bars represent one standard deviation, and each fade dot represents a single replicate value. c, d Standardized effect of hd, 〈l〉, and rΘ on QST,s, and QST,u obtained from a multivariate regression model independently fitted for low and high migration regimes on average results from 5 replicate simulations of the IBM on all undirected connected graphs with M = 7 vertices and varying rΘ (see Methods). The ambivalence of the effect of rΘ on QST,u found for the path graph holds for general graph ensembles and adds up to that of 〈l〉 and hd. Error bars show 95% confidence intervals. Analogous results on graphs with M = 9 vertices are presented in Supplementary Fig. 7 and all regression details can be found in Supplementary Table 3.

We extend our analyses to realistic landscapes with a continuum of habitat types by running simulations on graphs obtained from real spatial habitat datasets and by considering mean annual temperature as a proxy for habitat type (see Supplementary Fig. 8 and Supplementary Table 4). We also consider simulations accounting for trait-dependent competition to test whether our results hold under more complex ecological processes (see Supplementary Note for the implementation details and Supplementary Table 5 for the results). The simulations are congruent and show that the effects of rΘ, hd, and 〈l〉 are similar under these alternative settings, underlining the robustness of these metrics and the generality of our conclusions. Taken together, these results indicate that under sufficiently strong selection and sufficiently high habitat heterogeneity, adaptive differentiation QST,s is mainly driven by habitat assortativity rΘ. Nonetheless, local adaptation is lost in disassortative graphs when m > m, such that 〈l〉 and hd become complementary determinants of QST,s for high migration regimes.

Effect of habitat assortativity on neutral differentiation under heterogeneous selection

We finally consider a setting with heterogeneous selection where individuals carry both neutral and adaptive traits. With distinct habitat types, selection promotes neutral differentiation by reducing the birth rate of maladaptive migrants, reinforcing the isolation of local populations. We have shown above that adaptive differentiation QST,s is driven by habitat assortativity rΘ, so we expect rΘ, together with the topological metrics found in the setting with no selection, to influence the level of neutral differentiation QST,u. We first investigate how the response of QST,u to migration compares between the setting with no selection and the setting with heterogeneous selection for graphs with an identical topology. We then examine how the response compares between graphs with an identical topology but different rΘ. We finally consider simulations on different graphs with varying rΘ to assess the concurrent effect of 〈l〉, hd, and rΘ on QST,u.

Migration has a fitness cost because maladaptive migrants present lower fitness. Under an equivalent migration regime, migrants therefore have a lower probability of reproduction, increasing the populations’ isolation compared with a setting without selection. Simulations with varying m on the complete graph confirm that selection in heterogeneous habitats reinforces QST,u compared with a setting without selection (Fig. 5a). Nonetheless, previous results show that adaptive differentiation QST,s vanishes on a disassortative graph when m > m, implying that individuals become equally fit in all habitats. In this case, the isolation effect of heterogeneous selection is lost and QST,u reaches a similar level as in the setting with no selection for m > m (Fig. 5a), although QST,u is slightly higher in the setting with heterogeneous selection due to lower population size (\(\overlineN=bK(1-p\theta )\) vs. \(\overlineN=bK\), see section above and Methods). This suggests that rΘ reinforces QST,u, as assortative graphs sustain higher levels of adaptive differentiation (Figs. 3 and 4). Simulations on the path graph with varying Θ-spatial distribution support this conclusion for high migration regimes, but show the opposite relationship under low migration regimes, where the habitat assortativity rΘ decreases QST,u (Fig. 5b). Assortative graphs are composed of large clusters of vertices with similar habitats, within which migrants can circulate without fitness losses. Local neutral trait distributions become more correlated within these clusters, resulting in a decline in QST,u for assortative graphs compared with disassortative graphs. Figure 5b therefore highlights the ambivalent effect of rΘ on QST,u. rΘ reinforces QST,u by favouring adaptive differentiation, but also decreases QST,u by decreasing population isolation within clusters of vertices with the same habitat type.

We compare the effect of rΘ on QST,u to the effect of the topology metrics 〈l〉 and hd found in the setting with no selection using multivariate regression analysis on simulation results obtained for different graphs with varying Θ-spatial distribution (Fig. 5d for graphs with M = 7 vertices and Supplementary Fig. 7b for M = 9). The multivariate model explains the discrepancies in QST,u across the simulations for low and high migration regimes (see Supplementary Table 3 for details), and we find that rΘ, 〈l〉, and hd contribute similarly to neutral differentiation. Hence, the effects of rΘ and the topology metrics 〈l〉 and hd add up under heterogeneous selection. A change in sign of the standardized effect of rΘ on QST,s for low and high migration regimes verifies that the ambivalent effect of rΘ on QST,u found on the path graph holds for general graph ensembles. Simulations with trait-dependent competition and simulations on realistic graphs with a continuum of habitat types equally confirm the ambivalent effect of rΘ and further support the complementary effect of 〈l〉 and hd on QST,u (see Supplementary Fig. 8). 〈l〉 and hd therefore drive neutral differentiation with and without heterogeneous selection. rΘ becomes an additional determinant of neutral differentiation under heterogeneous selection. In contrast to the non-ambivalent, positive effect of habitat assortativity on adaptive differentiation, rΘ can amplify or depress neutral differentiation depending on the migration regime considered.