Dynamical community structure of populations evolving on genotype networks

2.1. Definitions and dynamical equations

Each genotype i in the network is represented by $n_i(t)$ sequences at time $t\text{,}i=1\text{,}\dots \text{,}m$, with m the number of different genotypes (or nodes) in the network. A fixed population size is used, $N={\sum }_i{n}_i(t)$, and we work in the limit of infinitely large populations, $N\to \infty$, so that finite size effects are discarded and, in practice, we consider the fraction of population at each node. The initial distribution of sequences on the network at $t=0$ is $\overrightarrow{n}(0)$. In particular, we will be interested in homogeneous initial populations, that is, in populations which, at time $t=0$ are formed by N identical genotypes. The degree $k_i$ of each genotype i specifies how many neutral neighbors it has.

At each time step, all sequences replicate synchronously. Daughter sequences mutate to one of the $3l$ nearest neighbors with probability $\mu$, and remain equal to their mother sequence with probability $1-\mu$. In our representation $0<\mu ⩽1$. The singular case $\mu =0$ is excluded to avoid trivial dynamics and guarantee evolution towards a unique equilibrium state. With probability $k_i/(3l)$, the mutated sequence exists in the neutral network and it adds to the population of the corresponding neighboring genotype. Otherwise, it falls off the network and disappears. These conditions implicitly represent a peak landscape in phenotypes, that is a phenotype with value 1 and any other possibility with value 0 [6].

The mean-field equations describing the dynamics of the population on the network in matrix form are(1)$\overrightarrow{n}(t+1)=(2-\mu )\mathbf{I}\overrightarrow{n}(t)+\frac{\mu }{3l}\mathbf{C}\overrightarrow{n}(t)\text{,}$where $\mathbf{I}$ is the identity matrix. The transition matrix $\mathbf{M}$ is defined as(2)$\mathbf{M}=(2-\mu )\mathbf{I}+\frac{\mu }{3l}\mathbf{C}\text{.}$

The set of m eigenvalues (all real) of $\mathbf{M}$ is $\{{\lambda }_i\}$, and they are ordered such that ${\lambda }_i⩾{\lambda }_{i+1}$. The corresponding m eigenvectors are $\{{\overrightarrow{u}}_i\}$, and since $\mathbf{M}$ is real and symmetric they can be chosen such that ${\overrightarrow{u}}_i·{\overrightarrow{u}}_j=0\text{,}\forall i\ne j$ and $|{\overrightarrow{u}}_i|=1\text{,}\forall i$. Matrix $\mathbf{M}$ is irreducible by definition (the underlying network is connected) and has positive values in the diagonal. It is therefore primitive, so the Perron–Frobenius theorem guarantees that the largest eigenvalue of $\mathbf{M}$ is positive, ${\lambda }_1>|{\lambda }_i|\text{,}\forall i>1$, and its associated eigenvector is also positive (i.e., ${({\overrightarrow{u}}_1)}_i>0\text{,}\forall i$) in the interval of $\mu$ values used [12].

The dynamics of the system, Eq. (1), can thus be written as(3)$\overrightarrow{n}(t)={\mathbf{M}}^t\overrightarrow{n}(0)={\displaystyle \sum _{i=1}^m}{\lambda }_i^t{\alpha }_i{\overrightarrow{u}}_i\text{,}$where we have defined ${\alpha }_i$ as the projection of the initial condition on the ith eigenvector of $\mathbf{M}$,(4)${\alpha }_i=\overrightarrow{n}(0)·{\overrightarrow{u}}_i\text{.}$

Furthermore, as ${\lambda }_1>|{\lambda }_i|\text{,}\forall i>1$, there exists a unique asymptotic state of the population that is independent of the initial condition $\overrightarrow{n}(0)$ and is proportional to the eigenvector that corresponds to the largest eigenvalue, ${\overrightarrow{u}}_1$:(5)${\displaystyle \underset{t\to \infty }{\lim }}\left.\frac{\overrightarrow{n}(t)}{{\lambda }_1^t{\alpha }_1}\right)={\overrightarrow{u}}_1\text{,}$while the largest eigenvalue ${\lambda }_1$ yields the growth rate of the population at equilibrium (in the absence of rescaling). For convenience, in the following, and without any loss of generality, we normalize the population $\overrightarrow{n}(t)$ such that $|\overrightarrow{n}(t)|=1$ after each generation. With this normalization, $\overrightarrow{n}(t)\to {\overrightarrow{u}}_1$ when $t\to \infty$.

It is easy to demonstrate that the eigenvalues ${\lambda }_i$ of the transition matrix $\mathbf{M}$ are related to the eigenvalues ${\gamma }_i$ of the adjacency matrix $\mathbf{C}$ through ${\lambda }_i=(2-\mu )+\frac{\mu }{3l}{\gamma }_i$. Furthermore, the eigenvectors of both matrices are identical [12]. This result implies that the asymptotic state of the population only depends on the topology of the genotype network.

2.2. Time to equilibrium

Eq. (3) describes the dynamics towards equilibrium from an initial condition $\overrightarrow{n}(0)$. The distance $\mathrm{\Delta }(t)$ to the equilibrium state can be written as(6)$\mathrm{\Delta }(t)\equiv \left.\frac{{\mathbf{M}}^t\overrightarrow{n}(0)}{{\lambda }_1^t{\alpha }_1}-{\overrightarrow{u}}_1\right|=\left.{\displaystyle \sum _{i=2}^m}\overrightarrow{{\mathrm{\Delta }}_i}(t)\right|=\left.{\displaystyle \sum _{i=2}^m}\frac{{\alpha }_i}{{\alpha }_1}{\left.\frac{{\lambda }_i}{{\lambda }_1}\right)}^t{\overrightarrow{u}}_i\right|\text{.}$In order to estimate how many generations elapse before equilibrium is reached, we fix a threshold $∊$, and define the time to equilibrium $t_{∊}$ as the number of generations required for $\mathrm{\Delta }(t_{∊})<∊$.

When ${\alpha }_2\ne 0\text{,}{\lambda }_2\ne 0$ and ${\lambda }_2\ne {\lambda }_3\text{,}{t}_{∊}$ can be approximated to first order by(7)$t_{∊}^1\simeq \frac{\ln |{\alpha }_2/{\alpha }_1|-\ln ∊}{\ln |{\lambda }_1/{\lambda }_2|}\text{.}$This approximation turns out to be extremely good in most cases thanks to the exponentially fast suppression of the contributions due to higher-order terms (since ${\lambda }_i⩾{\lambda }_{i+1}\text{,}\forall i$). An evaluation of situations where approximation (7) fails can be found in [12]. However, all cases shown in this work are well approximated by (7) within an error of one generation with respect to the exact time to equilibrium implicitly defined through expression (6).

An explicit relationship between the time to equilibrium and the mutation rate can be obtained by expanding Eq. (7) in powers of $\mu$,(8)$t_{∊}^1=\ln \left.\left.\frac{{\alpha }_2}{∊{\alpha }_1}\right|\right)\left.\frac{a}{\mu }+b-c\mu \right]+O({\mu }^2)\text{,}$where(9)$a=\frac{6l}{({\gamma }_1-{\gamma }_2)}\text{,}b=\frac{{\gamma }_1+{\gamma }_2-6l}{2({\gamma }_1-{\gamma }_2)}\text{,}c=\frac{{\gamma }_1-{\gamma }_2}{72l}\text{.}$

Since $c\ll a$, the dependence of the time to equilibrium with the mutation rate follows $t_{∊}^1\propto {\mu }^{-1}$ [12]. For a fixed topology of the genotype network, the mutation rate $\mu$ sets the rate at which equilibrium is approached. An additional factor, that we will explore in the following, is the effect of the initial condition on $t_{∊}^1$, which is implicit in the projections ${\alpha }_1$ and ${\alpha }_2$, as previously defined.

3.1. Equilibrium properties

In all the results that will be presented, we work in the limit of infinite population size and thus obtain the dynamical properties of the system by numerically solving the corresponding equations. To evaluate equilibrium properties, as discussed, it is enough to consider the adjacency matrix $\mathbf{C}$. The first eigenvector of the transition matrix $\mathbf{M}$ yields the population of each node at equilibrium, according to (5). Two important properties of genotypes are their degree ($k_i$, number of neutral neighbors) and the fraction of population they accumulate at equilibrium, ${({\overrightarrow{u}}_1)}_i$. We will refer to them as node degree and node population throughout the paper. When these two quantities are represented as a function of each other, we observe the appearance of three disjoint clusters of nodes, Fig. 2(A). From now on, these clusters will be communities 1, 2, and 3 (C1, C2, and C3 for short), ranked in decreasing order with respect to the maximally populated genotype (or node). The separation in three communities is also observed when we represent a histogram of the number of nodes as a function of their population, Fig. 2(B). We obtain a distribution with three well-defined regions of abundance, in agreement with the clusters in Fig. 2(A). Finally, a third indicator of the division in communities appears in a rank-ordering representation of nodes according to their population, Fig. 2(C). There are visible jumps in this curve that coincide with changes of community.

1. Download : Download full-size image

Fig. 2. Equilibrium composition of a neutrally evolving population of RNA sequences: dynamical community indicators. (A) Node degree versus node population for all $22\text{,}434$ sequences in the network associated to phenotype ..((….))..; three clusters corresponding to three communities identified through the dynamics are observed. (B) Histogram of the population of nodes. (C) Rank plot of the node population. The population in node (or genotype) i is obtained as the ith element of the first eigenvector of the transition matrix $\mathbf{M}$.

3.2. Time to equilibrium

The previous analysis of equilibrium properties reveals that, regarding dynamical properties, genotype networks seem to organize in communities. In order to further investigate this possibility and the hierarchical organization of the network, we continue our analysis with the study of non-equilibrium properties, and their dependence on quantities characterizing the nodes of the network. A first question is whether the initial condition affects in a significant and meaningful way the time to equilibrium. From a biological viewpoint, this quantity is of interest as well if we consider how an evolving population might find and fix a new phenotype. Previous phenotypes correspond to the exterior of the network we are investigating, and the new phenotype is associated to the current network. The probability to enter the network through a particular node depends on the number of outgoing links it has, that is its outward degree $3l-k_i$. Also, the larger the network the more likely that the phenotype it represents is localized. Nodes that are more connected in the network have lower probability of being the first node visited by an external population. This fact has implications that have been discussed elsewhere [19]. In addition, it is also common that populations enter a new phenotype through a single node, affected by a sort of genomic bottleneck caused by the difficulty to find and fix new phenotypes, and also usually by evolving at not-that-high mutation rates (specifically, in the limit $\mu N\ll 1$) [20].

The time to equilibrium varies approximately twofold depending on the initial distribution of the population. We have compared $t_{∊}^1$, Eq. (7), in several different situations, in particular when all the initial population is concentrated in a single node of communities 1, 2, or 3, or when it spreads uniformly over the whole network. Further, we have also considered whether C1, which is the subset of the network accumulating most of the sequence population (see Fig. 2), could be a good representative of the dynamics we observe over the whole network. To this end, we have simulated the evolution towards equilibrium in a population spreading only on C1 (with 7460 nodes) and compared with the previous cases. All these results are summarized in Fig. 3(A). The slowest dynamics corresponds to initial conditions starting in C3, pointing out at its relative isolation within the network, while an initial condition uniformly distribution over the whole network leads to the fastest dynamics. The comparison between the subnetwork C1 isolated from the rest (sC1) and the whole network reveals that dynamics is faster in the former, though both are quite comparable. Differences arise from a variation in the two largest eigenvalues, as shown in Fig. 3(B).

1. Download : Download full-size image

Fig. 3. Dependence of time to equilibrium on the initial condition. (A) Time to equilibrium as a function of the mutation rate $\mu$ when the evolutionary process takes place only in isolated community 1 (sC1, green dashed curve) or in the whole network (remaining curves) for a sample of different initial conditions. Time to equilibrium is shorter in the former case up to a factor 2. Different situations are detailed in the main text and in the legend. We observe a significant variation among nodes of C2, two examples are given. (B) Comparison between the first and second eigenvalues ${\lambda }_1$ and ${\lambda }_2$ of the whole network and the isolated C1 as a function of $\mu$. (C) Time to equilibrium versus node degree when the population initially occupies a single node i, as a function of the degree $k_i$ of that node. (D) Time to equilibrium versus equilibrium population ${({\overrightarrow{u}}_1)}_i$. C2 is clearly split in two subcommunities. The coloring of each node in (C) and (D) follows the division presented in Fig. 2(A). Parameters are $∊=0.0001$ in (A-D) and $\mu =0.1$ in (C) and (D). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Next, we have represented the time to equilibrium $t_{∊}^1$ as a function of the degree $k_i$ of each node i, fixing the mutation rate $\mu =0.1$ and repeating the calculation above for all nodes in the whole network. The relation between these two quantities offers little discriminatory power regarding community structure, since only two major communities can be resolved, see Fig. 3(C). On the contrary, a plot of the time to equilibrium as a function of the population at equilibrium (now relating two dynamical quantities) reveals a rich structure in communities and hints at a possible hierarchical organization, Fig. 3(D). C2 is now clearly separated into two independent communities (labeled C2A and C2B) and further subdivisions of C3 and C2A can be hypothesized.

Fig. 4 compares the overall properties of the four communities C1, C2A, C2B, and C3 clearly detected in Fig. 3(D). It is interesting that three of these quantities appear correlated (the fraction of the total population in each community, the average — per node — population, and the average degree), while two others are weakly dependent on each other (size of the community and the average time to equilibrium) and apparently uncorrelated to the former three.

1. Download : Download full-size image

Fig. 4. Overall properties of the communities detected in Fig. 3(D). While C1 is formed by only 34% of the nodes (A), it attracts 57% of the population at equilibrium (B) and reaches the largest average population per node at the equilibrium (C). In fact, its largest average degree per node (which correlates with higher folding stability) makes C1 more robust under mutations and an attractor for the population (D). The time to equilibrium depends in a complex manner on the location of nodes within communities, as well as on the particular topology of the network (E). C3 behaves quite differently, since it is a relatively large community (A) but appears as isolated (E) and demonstrates little attracting power (B,C), due among others to its low average degree (D).

3.3. Topological communities

Our results up to now indicate that populations evolving on genotype networks organize in communities that are revealed through dynamical quantities, notably the time to achieve equilibrium when a population enters the network from a single node and the final population that same node (or genotype) succeeds at attracting at equilibrium. Now we wish to compare these dynamical indicators of community structure with two other methods to detect communities based solely on network topology [21], [22].

The modularity Q of a given network is quantified as the fraction of edges that fall within the groups specified by a particular division in communities minus the expected fraction should those edges be distributed at random. Therefore, communities are defined as sets of nodes sharing more links among them than with nodes outside the community. Formally, modularity is defined as(10)$Q={\displaystyle \sum _{v=1}^c}\left.e_{\mathit{vv}}-a_v^2\right)\text{,}$where $e_{\mathit{vv}}$ is the fraction of edges with both end vertices in the same community v and $a_i$ is the fraction of ends of edges that are attached to vertices in community v (with origin in nodes outside v). The sum runs over the different c communities in the partition tested. Optimal divisions into communities are obtained by maximizing the value of Q [23], [24].

We have implemented two different methods to detect topological communities in our genotype network. First, we have used the stochastic block model inference method in [25].³ It is based on a nested generative model that uses a hierarchical characterization of the entire network at different scales and allows to perform a correct statistical inference and a proper detection of its modular structure. This method permits to establish a priori the number of communities in the network, so we have fixed it to 4 in order to compare the result with the dynamical communities in Fig. 5(A). The obtained modularity is $Q=0.66$, which is not optimal (since we fix the number of communities) but takes a value higher than in the dynamical case (with $Q=0.46$). However, the stochastic block model method mixes nodes from all four communities, as shown by the distribution of different colors in Fig. 5(B). Second, we use a Potts model approach [26] which optimizes Q and freely selects the corresponding number of communities. In this case we obtain a division in 12 communities and an optimal modularity $Q=0.74$. As with the previous topological method, we observe a high mixing of all 12 communities in Fig. 5(C) in comparison to the dynamical result shown in Fig. 5(A).

1. Download : Download full-size image

Fig. 5. Network communities obtained with different methodologies. The three panels repeat the values plotted in Fig. 3(D) and color the nodes according to the community they are assigned to with each method. (A) The functional dependence between the time to equilibrium and node population retrieves four communities (for comparison, this is a repetition of Fig. 2(D) with four colors). This partition yields a modularity value $Q=0.46$. (B) Communities obtained with a stochastic block model inference method fixing the number of communities to 4. Here, $Q=0.66$. (C) Community structure detected through a modularity optimization method. In this case, the number of communities is not fixed, so modularity reaches a higher value, $Q=0.74$. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

3.4. Dynamical communities and genotype composition

Methods to optimize topological communities yield divisions that differ substantially from those obtained through population dynamics. In order to disentangle the reasons that concentrate populations in certain regions (as C1) while others, though significantly large are comparatively depleted in population (as C3), we have investigated the existence of compositional differences among the four communities in Fig. 5(A). To this end, we need here to reconsider the system that we began with, RNA sequences, and analyse differences and similarities between the composition, in terms of nucleotides, of sequences in each community.

As we illustrated in Fig. 1, not all positions along the RNA sequence admit mutations with the same probability. In the phenotype of the network that we have analysed, corresponding to the secondary structure ..((….)).., there are four positions that are invariant in all network genotypes. These are those forming pairs, which occupy positions 3, 4, 9, and 10, with nucleotides C, C, G, and G, respectively. In order for the phenotype to be maintained, there are some restrictions on the nucleotides occupying positions 2 and 11: they have to be occupied by nucleotides unable to pair; otherwise, a third pair would be stably formed, leading to a different phenotype. This condition excludes six combinations (A-U, U-A, G-U, U-G, G-C, and C-G) for positions 2 and 11. The remaining positions (1, 5, 6, 7, 8, and 12) can be occupied by any possible nucleotide, since they do not participate in structural changes — therefore the phenotype is maintained.

A compositional analysis of genotypes in the network reveals that it is precisely the limitations in the possible pairs at positions 2 and 11 that separates nodes into major communities. Fig. 6 illustrates the relationship between communities and sequence composition. Thanks to this analysis we have realized that three of the communities can be further divided into subcommunities attending to their composition, and that non-trivial dynamical relationships appear. We have calculated the Hamming distance, defined as the number of differences in the composition of two sequences, between all possible pairs in the network. It turns out that some communities are two mutational steps away, thus requiring an intermediate group of genotypes for a population to move from one group to another. The most remarkable example relates communities 2A and 2B, which can only be linked through nodes at communities 1 or 3. Further, if going through C3 the jump can only take place through sequences in community 3.3 going to 2A.2, since all other possibilities are again two mutations away, and thus an intermediate (sub) community is again required. Also, a similar decoupling occurs among the subcommunities of 2A: 2A.2 is two mutations away from 2A.1, but intermediate sequences reside in different communities. Therefore, if staying in 2A the only possible path has at least three mutations: 2A.2 $\to$ 2A.4 $\to$ 2A.3 $\to$ 2A.1.

1. Download : Download full-size image

Fig. 6. Genotype composition and community structure. Sequences in each of the four dynamical communities belong to different compositional groups, as specified. Colors as in Fig. 5(A). Compositional differences and Hamming distances indicate that a further community structure can be identified, in agreement with the existence of additional divisions, first inferred from Fig. 3(D) and here indicated as numbered circles inside major communities. Links (black between major communities, white between subcommunities) join compositional groups that are at a Hamming distance of 1.

The process of division in communities attending to sequence composition could be in principle iterated to unveil a complete hierarchical structure in genotype networks. However, not all divisions are equally meaningful. For instance, the fact that most unpaired nucleotides inside loops of sizes 3 and 4 can take any value without disturbing the secondary structure makes divisions at this level uninteresting regarding population dynamics. From what has been discussed in this section, compositional restrictions in sequence positions that affect RNA secondary structure generates communities that are dynamically detected by evolving populations: these communities are de facto separated by valleys, or bottlenecks, in the space of genomes, a situation that difficults the exchange of population among groups.