PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Evolution. Author manuscript; available in PMC 2010 November 22.
Published in final edited form as:
Evolution. 2009 February; 63(2): 418–431.
doi:  10.1111/j.1558-5646.2008.00577.x
PMCID: PMC2989730
NIHMSID: NIHMS249698

DYNAMICS OF HYBRID INCOMPATIBILITY IN GENE NETWORKS IN A CONSTANT ENVIRONMENT

Abstract

After an ancestral population splits into two allopatric populations, different mutations may fix in each. When pairs of mutations are brought together in a hybrid offspring, epistasis may cause reduced fitness. Such pairs are known as Bateson-Dobzhansky-Muller (BDM) incompatibilities. A well-known model of BDM incompatibility due to Orr suggests that the fitness load on hybrids should initially accelerate, and continue to increase as the number of potentially incompatible substitutions increases (the “snowball effect”). In the gene networks model, which violates a key assumption of Orr’s model (independence of fixation probabilities), the snowball effect often does not occur. Instead, we describe three possible dynamics in a constant environment: 1) Stabilizing selection can constrain two allopatric populations to remain near-perfectly compatible. 2) Despite constancy of environment, punctuated evolution may obtain; populations may experience rare adaptations asynchronously, permitting incompatibility. 3) Despite stabilizing selection, developmental system drift may permit genetic change, allowing two populations to drift in and out of compatibility. We reinterpret Orr’s model in terms of genetic distance. We extend Orr’s model to the finite loci case, which can limit incompatibility. Finally, we suggest that neutral evolution of gene regulation in nature, to the point of speciation, is a distinct possibility.

Keywords: hybrid incompatibility, speciation, parallel evolution, neutral networks, drift, genetic distance, gene networks

Darwin (1859) compared attempts to cross two varieties of a single species (in many cases producing hybrid offspring of reduced fertility or viability) and attempts to cross two truly distinct species (producing almost no viable, or fertile, hybrid offspring). He came to the general conclusion that, “… there is no fundamental distinction between species and varieties,” in the sense that the two cases differ only in the degree to which similar mechanisms cause infertility or inviability (with rare exception). Darwin also recognized, but never solved, the problem of how hybrid infertility or inviability might increase in apparent opposition to selection, “inasmuch as the sterility of hybrids could not possibly be of any advantage to them, and therefore could not have been acquired by the continued preservation of successive profitable degrees of sterility.” Dobzhansky (1936) and Muller (1942) used a genetic model to solve the paradox. In this model, pairs of incompatible alleles can arise by two mutations in two allopatric populations. The important feature is that the alleles are deleterious only in combination, and that this combination is never “tested” by natural selection until the alleles are finally brought together in an unfortunate hybrid individual. Pre-existing reproductive isolation between two populations would allow the accumulation of such incompatibilities, whose collective effect could make reproductive isolation permanent even if the original cause of isolation were removed. This model had been previously proposed by Bateson (1909), though apparently neither Dobzhansky nor Muller were aware of this (Orr 1996). Such pairs (or triplets, etc.) of alleles are therefore referred to as Bateson-Dobzhansky-Muller incompatibilities (BDMIs).

Orr (1995) derived a relationship between K, the number of substitutions since the ancestral population split into two allopatric populations, and L, the expected “fitness load” suffered by hybrids in terms of the average deleterious effect [r with macron] of each such potentially incompatible pair of alleles:

L=1w1j=1K(1r¯)j11eK2r¯/2
(1)

For K2[r with macron] [double less-than sign] 1 – the start of the speciation process – Orr suggests that the strength of reproductive isolation should increase nearly as fast as the square of the number of substitutions:

LK2r¯/2
(2)

This initial rapid acceleration of the fitness load on hybrids is known as the “snowball effect” (Orr and Turelli 2001). As speciation continues, L asymptotically approaches one, corresponding to complete reproductive isolation. Orr’s result does not refer to any particular mode of fixation (e.g., under stabilizing selection, directional selection, or neutral drift) but assumes that fixations of alleles are independent and equiprobable events, and that one locus can experience only one substitution in either population, which eliminates the possibility of parallel genetic evolution.

Here we explore some violations of Orr’s assumptions: specifically, we allow the fixations of different alleles to have different probabilities, and those fixation probabilities not to be independent, but instead to depend on what fixations have previously occurred. In addition, we permit multiple substitutions across two populations at the same locus, thereby allowing parallel evolution. We do not seek to falsify Orr’s (1995) results, but rather to extend them.

Our experimental model for this exploration is in silico populations of “gene networks.” This model has previously been used to study evolvability and canalization (Wagner 1996; Siegal and Bergman 2002). Gene networks model the expression dynamics of sets of genes that are up- and down-regulated by one another’s gene products. The model permits complex epistasis and pleiotropy, hence fixation probabilities are not independent.

In this article, we consider the case of a constant environment, identical in the two allopatric populations. We choose this focus because the case of directional selection has already been covered well using a similar model: Johnson and Porter (Johnson and Porter 2000; Porter and Johnson 2002; Johnson and Porter 2007) showed that linear gene regulatory pathways can accumulate BDMIs in the presence of directional selection (a slowly ever-receding phenotypic optimum, which induces the population mean phenotype to “chase” the optimum). That directional selection would accelerate speciation is perhaps intuitive, because it results in substitutions, but theoretical predictions indicate that hybrid incompatibility should also increase due to substitutions occurring by drift alone, even under conditions of stabilizing selection (i.e., Orr’s model does not refer to the cause of the substitutions). Notably, Johnson and Porter did not observe accumulation of hybrid incompatibility under stabilizing selection. “The interaction of genetic drift and mutation, even at very high rates, did not reduce hybrid fitness at all on the time-scales we considered.” (Johnson and Porter 2000) In contrast, using the gene networks model, we demonstrate the growth of hybrid incompatibility (as well as other dynamics), even in a constant environment. We attribute this difference to the multi-locus epistasis and abundant antagonistic pleiotropy present in the gene networks model, but not in the linear regulatory pathways of the Johnson and Porter model. Johnson and Porter (2007) did include a single point of pleiotropic influence; however, much larger neutral networks are apparently required for the effects we describe to emerge. The gene networks model also easily generates hybrid incompatibility under directional selection, but that is not the focus of this article.

To summarize our results, we find that the gene networks model does not in general obey the predictions of Orr (1995). The hybrid fitness load L does not show an initial increase proportional to K2, and, in most cases, does not asymptotically approach one, as predicted by Equations 1 and 2, above. Instead, depending on the several parameters of the model, under constant selection, we observe three dynamics: 1) hybrid incompatibility between two allopatric populations may not increase at all; or 2) it may increase to large values; or 3) a pair of populations may “drift” in and out of compatibility.

The first of these occurs with certain model parameter values for which a constant environment produces a basin of attraction toward an optimal genotypic state. In this case, a constant environment is both stabilizing and prevents the development of hybrid incompatibility. This affirms an observation of Johnson and Porter, but conflicts with Orr’s predictions.

The second dynamic, accumulation of hybrid incompatibility to high levels, is observed when punctuated equilibria (Bergman and Feldman 2003) occur over an extended period despite a constant environment. During this process, the mean fitness intermittently increases as the mean phenotype jumps toward, but does not reach, a static optimum. The third dynamic, drift of a pair of populations in and out of compatibility, occurs in a constant environment when selection is truly stabilizing. In this case, the mean phenotype and mean fitness do not change, but “developmental system drift” (True and Haag 2001) occurs, permitting hybrid incompatibility to arise. (This is the case that Johnson and Porter did not observe in their model.)

To make the distinction clear: a constant environment does not necessarily imply stabilizing selection. If only deleterious variation is present, we say that constant selection is stabilizing: it removes the deleterious variation, selecting for the current, relatively adapted phenotype. However, in the same constant environment, if adaptive variation occurs rarely, then periods of stasis will be interrupted by periods of adaptation. Pleiotropy and complex epistasis can permit the second dynamic to occur. Furthermore, antagonistic pleiotropy increases the size of the near-neutral networks in the genetic neighborhood of local fitness optima, permitting the third dynamic to occur by developmental system drift.

We also extend the model of Orr to the case that the number of interacting loci,, is finite, in which case the growth of hybrid incompatibilities may be curtailed.

Parallel Evolution

We define parallel evolution as the independent evolution of similar characters in allopatric populations sharing a common ancestral state. Parallel evolution may be genetic or phenotypic; where relevant, the distinction is indicated below.

A central hypothesis of this article is that parallel evolution can prevent the development of hybrid compatibility. Parallel genetic evolution (absence of divergent substitutions in two populations) acts directly to prevent the development of BDMIs. Parallel phenotypic evolution (parallel evolution in the external phenotype) may or may not entail parallel genetic evolution. If it does (see Figure 4, for example), then the growth of hybrid incompatibility may be curtailed. Parallel phenotypic evolution may be caused by directional selection, which can strongly drive two allopatric populations in the same evolutionary direction. Alternatively, even if two populations are evolving neutrally, variational constraints (the topology of the neutral network, in the neutral case) dictate the evolutionary pathways along which they may drift. Furthermore, even when there is no fitness increase in individual phenotypes, there can still be indirect selection for mutational robustness, such that certain parts of the neutral network will tend to be more populated (van Nimwegen et al. 1999), resulting in parallel genetic evolution, and curtailing hybrid incompatibility.

Figure 4
Parallel genetic evolution in two allopatric populations. (Parameters used: a=1, integer target, c=1.0, μ=0.01, N=100, Td =1,000. Truncation selection was applied.)

Methods

The Phenotype and the Genotype

In the gene networks model of Wagner (1996) and Siegal and Bergman (2002), the phenotype of an individual, shown on the right of Figure 1, is an M-dimensional vector representing the amounts of M gene products it produces. The heritable genotype of the individual is represented by an M × M matrix, shown on the left of Figure 1. Element (i, j) of the matrix represents the (positive, negative, or zero) regulatory action of gene product j upon the expression of gene i.

Figure 1
The “developmental process” of the gene networks model. The M × M genotype matrix ||gij|| and an initial vector are provided. The individual’s M-dimensional phenotype vector (elements p1 through pM) is produced by repeated ...

Initially, a random fraction c (the parameter c is fixed for any given experiment) of the matrix elements is drawn from the standard normal distribution; the rest (1-c) are set to zero, which permits some control of the amount of epistasis present in the network.

Segregation by Rows, and Mutation

Upon mating, an offspring is generated, row by row, by copying each new row from the matrix of one parent or the other with equal probability, akin to segregation of chromosomes (rows), with no recombination within chromosomes. Then mutation is applied at rate μ per genome and entails replacement of a random nonzero matrix element with another number from the standard normal distribution. Thus the per (nonzero) element mutation rate will be μ/(cM2).

The Developmental Process

At the beginning of the developmental process, each individual’s phenotype is set to a fixed initial vector, the elements of which are in the closed interval [0,1]. Each iteration of the developmental process, illustrated in Figure 1, consists of two parts: first, the genotype matrix and the current phenotype vector are multiplied; and second, the entries of the resulting vector are normalized. The following sigmoid normalization function is used to map each element x of the product back to the interval [0, 1]:

sig(x,a)=(1+eax)1
(3)

The vector of normalized elements becomes the phenotype vector for the next iteration of the developmental process. The parameter a (fixed for any given experiment) controls the steepness of the sigmoid function, which is plotted for two values of a in Figure 2.

Figure 2
The sigmoid normalization function sig(x, a) from Equation 3. Each element of the product of the genotype matrix and the current phenotype vector is normalized back to the interval [0, 1] with this function, producing the corresponding element of the ...

A high value of a (e.g., 10) creates a tendency for the elements of the phenotype vector to assume near-zero or near-one values. A value near 0.5 at an a of 10 requires a precise balance of regulatory effects. Low values of a (e.g., 1) permit elements of the phenotype vector to more easily assume intermediate values and make it difficult to attain near-zero or near-one values.

We restrict the phenotypic elements to the interval [0, 1], rather than the interval [−1, 1] used by Wagner (1996) and Siegal and Bergman (2002), because we interpret the value of each phenotypic element as the presence (positive value) or absence (zero value) of the corresponding gene product. Positive genotype matrix entries imply up regulation, and negative entries down regulation. Therefore, a mutation is required to reverse the direction of regulation; i.e., there is no critical point about which a slight change in the value of a phenotypic element will suddenly reverse the direction of regulation.

Integer- and Real-valued Phenotypic Targets

An additional model parameter is the “target type.” The elements of the initial phenotype vector, and of the target phenotype vector, may either be “integer-valued” (chosen with equal probability to be zero or one), or “real-valued” (chosen uniformly from the real numbers in the closed interval [0, 1]). The initial phenotype vector in a given experiment will have the same type as the target phenotype vector. A ten-dimensional integer-valued target is plotted as the ten large orange circles in both panels of the first row of Figure 3; the vector represented is (1, 0, 1, 0, 1, 1, 0, 0, 0, 1). A real-valued target is plotted as the ten large orange circles in each panel of the second row. (The other features of this figure are discussed in detail below.)

Figure 3
Phenotypic plots at generation 0 (left column) and generation 100 (right column) for four runs using four combinations of a values (a=1 or a=10) and target types (integer or real). The ten large orange circles in each plot represent the ten-dimensional ...

Because high values of a tend to produce phenotypes with near-one or near-zero elements, experiments with real-valued targets and a high value of a tend to have difficulty producing a phenotype near the target phenotype, but they easily approximate integer-valued targets. In contrast, if a=1, intermediate phenotype values are easily attained, but extreme matrix values are required to produce phenotypes with near-zero or near-one elements.

Viability, Fitness, and Mating

At the beginning of the developmental process, each individual’s phenotype is set to the fixed initial vector uinit. If the individual’s phenotype converges to any vector in less than 100 developmental iterations (rather than entering a cycle, or not converging), then that individual is declared to be viable; otherwise, it is removed from the population. Convergence is determined as follows. Let d be the distance metric between two phenotype vectors ua and ub defined by:

d(ua,ub)=i=1M(uiauib)2/M
(4)

At the end of each iteration, an average phenotypic vector uavg over the last τ=10 iterations is computed. If the average d between the last τ phenotype vectors and uavg is below the threshold value ε=10−4, then convergence has occurred, and the individual is declared viable.

At the beginning of each experiment, an M-dimensional phenotypic “target vector” is randomly chosen; this target defines the best possible phenotype. For all viable offspring, we compute a probability of survival to maturity (“fitness”), f, using the distance between the final, converged phenotype and the target vector:

f=ed(uphen,utarg)σ
(5)

The parameter σ permits control of the severity of selection; we use a value of σ=0.1 here.

At each generation, we repeat the following operations until there are exactly N surviving individuals: random mating with segregation of chromosomes; mutation; viability selection (i.e., whether the phenotype converges); fitness selection by Equation 5.

Population Metrics

We define I, the incompatibility between two allopatric populations, as:

I=12H/(W1+W2),
(6)

where W1 and W2 are the average fitnesses within the two allopatric populations, and H is the average fitness (counting nonviables as having zero fitness) of the hybrid offspring produced by mating pairs with one parent coming from each population. This is intended to be analogous to Orr’s hybrid fitness load L in Equations 1 and 2; thus Orr’s predictions should apply to I.

Define the “self-incompatibility,” Is, within a single population as:

Is=1Woffspring/W,
(7)

where Woffspring is the average fitness of offspring of the current generation (nonviables have zero fitness), and W is the average fitness of the current generation. Typically in our models, Is is positive immediately after the founding of a population, but it quickly decreases to near zero.

Gavrilets (2004) noted two stages in the evolutionary dynamics of a simple two-locus BDM model: the first, driven by strong selection, during which most of the population rapidly moves onto a neutral fitness ridge, and a second in which the dynamics are driven by weak forces like mutation. The first is analogous to the rapid elimination of Is, which we routinely observe. Nei et al. (1983) noted a similar effect in their models.

We measure Ki, the number of substitutions in population i since the population split, as follows. At the moment of population split, each Ki is set to zero. The frequencies of all polymorphisms at each matrix element of the N genotypes in each population are tracked. When the frequency of a particular allele first exceeds 0.95 in population i, Ki is incremented by one. That allele is then marked as fixed unless its frequency subsequently decreases below 0.5. In comparing two of the allopatric populations, i and j, we sum the number of substitutions between them: K = Ki + Kj.

Results

All experimental populations described in this article were created from six founder individuals that were randomly generated and subjected to one round of viability and fitness selection. The populations were thereafter maintained at the constant size of N=100 (Figures 3 and and4)4) or N=250 (all subsequent Figures involving simulation of gene networks) individuals. For Figures 3, ,4,4, ,5,5, ,6,6, and and8,8, animated versions of the figures are available in online supplementary material.

Figure 5
Incompatibility, I, for the ten replicates r1, r2, …, r10 as a function of K2,, for a=1. Left: integer-valued target; right: real-valued target. (c=1.0, μ=0.003, N=250, Td=1,000).
Figure 6
Incompatibility, I, for the ten replicates r1, r2, …, r10 as a function of K2, for a=10. Left: integer-valued target; right: real-valued target. (c=1.0, μ=0.003, N=250, Td =1,000).
Figure 8
Incompatibility, I, for the ten replicates r1, r2, …, r10 as a function of K2, for a=10 and real target, c=0.1 (μ=0.003, N=250, Td=1,000).

A View of the Phenotype for Four Parameter Sets

Figure 3 illustrates four experimental runs (four rows) of a single population at two points in time: generation 0 (left column) and generation 100 (right column). The four runs correspond to the four combinations of two values of the parameter a (a=1, or a=10) and two types of phenotypic target (integer-valued, or real-valued). The left plot in the first row shows this run at generation 0. Since a is low (a=1), the values of the elements of the phenotype vectors are well distributed through the intermediate values between zero and one. The dots are generally colored blue and blue-green, indicating that no phenotypes very near the target have yet evolved. In contrast, the right side of the first row of the figure shows generation 100 of the same run: the elements of the phenotypes are clustered closer to their target values, and the orange and red colors of the dots reflect the evolved closeness of each phenotype vector to the target vector.

The second row of Figure 3 shows a run with a=1 and a real-valued target. At generation 0, the dots representing the phenotypes are more green than blue. Because the elements of the target vector (orange circles) have intermediate rather than extreme values, the initial phenotype vectors are in general closer to the target from the outset. At generation 100, the phenotype vectors again cluster near the target vector, giving the orange-red colors of the dots.

The third row of Figure 3 shows a run with a=10 and an integer-valued target. The initial distribution (left) of the elements of the phenotype vectors is biased toward the extremes of 0 and 1 due to the sharpness of the sigmoid normalization curve (see Figure 2), which makes it difficult to obtain intermediate element values. At generation 100, the elements of the phenotypes are clustered about the elements of the target, as reflected by the dark red color of the dots. The high value of a causes a bias toward the production of phenotypes with extreme values, causing the population to evolve very fit phenotypes in this case.

The fourth row of Figure 3 has a=10 and a real-valued target. The initial distribution of phenotypes is similar to that of the third row with elements biased toward extreme values (0 or 1). However, in this case, the target has intermediate values, so the initial phenotypes are far from it, as reflected by the blue shades of the dots. By generation 100, the population is still unable to evolve a very fit phenotype: the colors of the dots are an intermediate green, and, in most cases, the phenotypic elements are not clustered near the target elements.

We will use these four combinations of a and the target type as our primary means of manipulating the evolutionary challenge faced by our populations. Later, we will also discuss the effect of varying the c parameter.

Parallel Evolution in the Genotype

Each experiment below starts with a single population. At generation Td (“time of divergence”), the population is cloned into two populations, between which there is no subsequent migration. We typically wait to split the population until it has become well adapted (attained near-optimal fitness), except when this is not possible due to very slow adaptation, as described below. This is also sufficient time to eliminate the initial Is that is due to the founding of the population with randomly generated individuals.

In order to generate a visually striking example of parallel genetic evolution, we increase the selective intensity by applying truncation selection. Specifically, in addition to our normal selective regime, the bottom 25% of individuals by fitness rank were not permitted to mate. Genotypes, rather than phenotypes, are plotted in Figure 4. At generation 0 (first column), the experiment starts with a single population. The ten × ten grid of colored squares represents the genotype matrix of one individual randomly selected from the population of 100 individuals. The darkest reds represent high-valued matrix elements (+3 and higher), and the darkest blues represent low-valued matrix elements (−3 and lower). Zero is represented by an intermediate green color. At generation 1,000 (second column), the population is cloned, forming two identical populations. The genotype of one randomly selected individual is displayed from the original Population 1 (top), and from the cloned Population 2 (bottom), respectively. There is no further migration between the two populations for the rest of the experiment.

By generation 10,000 (third column), it is visually apparent that parallel genetic evolution has occurred. The target vector for this run was the integer-valued target (i.e., the orange circles) shown in both panels in the first row of Figure 3. The patterns of red and blue rows of both of the genotype matrices in the last column of Figure 4 correspond exactly to the pattern of high and low elements of the target vector, i.e., (1, 0, 1, 0, 1, 1, 0, 0, 0, 1): high values in row i of the genotype matrix up-regulate gene product i, and low values do not. The standard fitness model (without truncation selection), used in the rest of our results, also produces parallel evolution, but not to such a visually striking degree as in Figure 4.

The Biological Interpretation of the Parameter ‘a

We interpret the parameter a as specifying the sensitivity of the regulatory response to the intermediate phenotype vectors that are computed during development. Do the genes act additively, at least in the middle of the range of their minimum and maximum combined effect (low a)? Or do they behave more like threshold characters, for which there is a sudden qualitative change when a threshold is passed (high a)?

When a=1, a broad basin of attraction is produced around local phenotypic optima; the population can quickly “hill climb” the fitness gradient, which is nonzero nearly everywhere, toward local optima, and to the global optimum. At the same time, a=1 makes it difficult to produce extreme (near zero or one) phenotypic values. There is, therefore, a distinct fitness peak with a very small neutral network at the top. For the metaphorical “fitness landscape”, a=1 corresponds to broad-based, sloping hills that are peaked, rather than flat, on top.

When a=10, it is much easier to produce extreme phenotypic element values, since for production of such extremes it is not necessary that multiple regulatory pathways act in the same direction. Therefore much larger neutral networks are created about the optima. At the same time, the basins of attraction about optima become much narrower, since for most of the input range of the sigmoid function, the gradient is near zero. For the metaphorical “fitness landscape”, there is a low plateau containing mesas: narrow elevated areas with steep sides and a flat top.

The parameter a could also be placed under evolutionary control, an analysis not pursued here.

I vs.K2

The extent of hybrid incompatibility was measured by our metric I, for the four combinations of a and target type illustrated in Figure 3. We performed ten replicate runs of eight parameter sets created using: two values of a, the steepness of the sigmoid curve (1, and 10); two target types (integer-valued, and real-valued); and two values of c (0.1, and 1.0). Following Johnson and Porter (2000), we used a population of N = 250 and a genomic mutation rate of μ = 0.003. The initial population was split into two allopatric populations of size 250 at time Td = 1,000.

All ten replicates of each parameter set were started with the same, fixed random seed. Thus, for each repetition of the same parameter set, each detail of the simulation (the phenotypic targets, the founder individuals, and all other individuals in the population at each generation) was identical, up to the moment of population split. At the moment of split, a new random seed was set, depending on the repetition number. This scheme permits us to eliminate variation between runs due to random events occurring before the population split. We further replicated all of the above (the ten replicates of eight parameter sets) with ten different initial random seeds (yielding a total of 800 runs), the choice of initial random seed had no significant effect. One set of 80 runs is shown below, but all sets are available in the online supplementary material.

We chose a Td of 1,000, since this delay allows all combinations of a and the target type (with the exception of a=10 and real-valued target) to attain high fitness before the population is split; it also allows the self-incompatibility Is to be largely eliminated. All experiments ran for Td + 30,000 generations.

The a=1 Case: Parallel Evolution Limits the Increase of Incompatibility

When a=1, there is a single evolutionary dynamic. Whether the target is real- or integer-valued, the population can quickly “hill climb” to the global optimum where there is little neutral drift. Instead strong stabilizing selection keeps the population near the optimum genotype, and pairs of allopatric populations are prevented from developing hybrid incompatibility.

On the left of Figure 5, we plot I against K2 for the ten runs for a=1, an integer-valued target. Clearly I does not increase as the square of K, as would be predicted by Orr. The range of the I values is very small: by the end of each experiment (corresponding to 31,000 generations, or K2 ~ 20,000–35,000 in the figure), I is less than 0.001, indicating that the populations are extremely compatible despite many generations of allopatric evolution. Speciation, which would require values of I approaching 1.0, does not occur. We expect strong parallel evolution for this parameter set, since, for a=1, extreme matrix values are required to produce near-zero- and near-one-valued phenotypic elements to match the integer-valued target elements. The average fitness (see online supplementary material) hovers around 0.98 to 0.99.

On the right of Figure 5, a similar plot is shown for the case of a=1 and real-valued target. Significant parallel evolution is still expected for this parameter set (because even a real-valued target has some extreme values, driving some parallel evolution); and incompatibility remains nearly as low as for the first parameter set: less than one percent by the end of each run. The average fitness (see online supplementary material) quickly attains 0.98 to 0.99 and remains at about that level on average.

The a=10 Case: Two Different Mechanisms Permit Hybrid Incompatibility

For a=10, two distinct evolutionary dynamics occur: 1) For integer targets, the variational bias towards extreme output values allows the population to quickly find the tops of the metaphorical “mesas” described above; the population then drifts about these large neutral networks. 2) In the case of real targets, variation when a=10 is so biased toward producing extreme phenotypic values that the population has great difficulty finding local optima; a dynamic of punctuated equilibria ensues, as occasional adaptation occurs, followed by long periods of stasis. Near-optimal fitness is never attained during the experiments. Under this dynamic, two allopatric populations may evolve asynchronously (one ahead of the other), or they may optimize the elements of the phenotype in differing order; either type of asynchrony allows the temporary development of hybrid incompatibility, which may later be reduced as the second population “catches up” to the first.

True and Haag (2001) describe a process they name “Developmental System Drift” (DSD) which means change in the developmental system (implying genetic change) even without phenotypic change. True and Haag state that since DSD happens by chance alone, rather than by selection, “drift” is an appropriate term. DSD does occur in the gene networks model. We suggest, however, that evolution of the developmental system can also be caused by indirect selection for robustness (van Nimwegen et al. 1999). Such evolution can be considered “neutral” in the sense that all (viable) phenotypes produced may be of identical fitness (probability of survival to adulthood) although, depending on the topology of the neutral network, certain genotypes may experience greater mutational load (producing a higher proportion of inviable mutants), resulting in an indirect form of selection.

Furthermore, evolution of the developmental system is also important during episodes of adaptation, not only during movement of the population on a neutral ridge. Consider a population split in two with the two allopatric populations subject to the same directional selection. If we observe identical fitness increase and parallel phenotypic evolution, this does not imply that the populations have also experienced parallel genetic evolution. If the two populations do now differ genetically, even if this difference is cryptic (not expressed in any phenotypic difference), then again there is the potential for hybrid incompatibility to arise. We have observed this effect in gene networks simulations with time-varying phenotypic targets (data not shown). Hence, we prefer to use the term “Developmental System Evolution” (DSE) rather than True and Haag’s DSD. We may then speak of “DSE under neutrality” or “DSE under phenotypic constancy” to indicate constancy of fitness, or phenotype, respectively, or “divergent DSE despite parallel phenotypic evolution” to indicate that two populations whose phenotypes have evolved in parallel are nonetheless genetically distinct.

On the left side of Figure 6, I is plotted against K2 for the case of a=10 and integer-valued target. Here I does increase to significant values (e.g., 20–40%) during several runs. However, increase in I often seems to be temporary, as I wanders to a peak, and then falls to low values (less than 5% incompatibility).

Fitness (see online supplementary material) is already near optimal (> 0.997) before the population splits, and remains extremely high. However, due to DSE, there is a neutral network about the global optimum, and drift about this neutral network allows the development of hybrid incompatibility, which can then fade away as pairs of populations drift back into compatibility. Hence, this mechanism of generation of incompatibility is potentially reversible; whether it reverses depends on the size and topology of the neutral network at the “top of the mesa”. Gavrilets (2004, p. 174) suggests that hybrid incompatibility should increase with genetic distance. (This relationship has been found in Drosophila (Coyne and Orr 1989; Coyne and Orr 1997), frogs (Sasa et al. 1998), birds (Lijtmaer et al. 2003), and Lepidoptera (Presgraves 2002), although it was not detected in darters (Mendelson 2003) and was found in only one of three genera of angiosperms (Moyle et al. 2004).) If genetic distance were to increase over time, as would be expected for very large neutral networks, then hybrid incompatibility should increase over time by drift. In contrast, we believe that the dynamic of “drifting” in and out of compatibility, as illustrated on the left side of Figure 6, is best explained in terms of random walks on neutral fitness ridges that are relatively confined, and where it cannot be assumed that genetic distance will tend to increase by drift. Rather, genetic distance may wax and wane as two populations randomly move relative to one another about a confined set of neutral ridges; hybrid incompatibility (as a roughly monotonic function of genetic distance) may therefore also wax and wane.

On the right side of Figure 6, I is plotted against K2 for a=10 with a real-valued target vector. Here, I does increase, in a volatile manner, to significant levels: 30–80%, depending on the run. Plotting the average fitness, W, against K2 (see online supplementary material), this volatility is clearly due to punctuated equilibrium effects. The average fitnesses of all populations show clear, persistent plateaus, punctuated when an adaptive variant is discovered by sudden jumps up to a higher plateau, which then persists for a time. In this case, despite a constant environment, the selective pressures are not constant: long periods of stasis due to stabilizing selection persist when adaptive variation is absent, but a new adaptive variant may sweep through its population, potentially causing a cascade of fixations that can contribute to incompatibility.

This incompatibility is also potentially reversible in an interesting way: suppose that two populations optimize the several elements of the phenotype one at a time (slowly, via punctuated equilibria), and in differing orders. Further suppose that after a very long time, they ultimately optimize all phenotypic elements. The intermediate state could include many possible phenotypes (e.g., there are many permutations of five out of ten “correct,” i.e., near-optimal elements, and a correspondingly large number of (potentially incompatible) genotypes that could produce these phenotypes. In contrast, there is only one way to get all ten elements “correct”, and fewer genetic alternatives to produce this phenotype. Therefore one could expect an initial increase in incompatibility, followed by a decline of incompatibility, despite allopatric conditions. (See also (Orr 2005) regarding the probability of parallel evolution.)

This situation (continuing rare adaptation, even over extended periods) may be common in biological systems in a constant environment. For example, in populations of E. coli that were maintained in a constant environment for more than 20,000 generations, adaptation continued throughout the entire period (Cooper and Lenski 2000). Parallel genetic and phenotypic evolution among the several populations was directly observed (Cooper et al. 2003; Philippe et al. 2007).

Significant incompatibilities will be frequent in such situations of evolutionary asynchronicity. Compatibility may, or may not, be attained by parallel evolution in the very long term. This will depend on the structure of the network of possible evolutionary pathways, in particular whether the pathways ultimately converge, or lead to distinct (and potentially incompatible) endpoints (see Discussion).

Extending Orr’s Model to the Case of Finite Loci

Here we present an extension of Orr’s (1995) analysis, which considers an infinite number of loci. Our simulations obviously involve a finite number of loci (specifically, cM2, or between 10 and 100 non-zero loci for the experiments shown here). In addition, Orr does not allow more than one substitution at a locus, whereas we do. We now extend Orr’s (1995) model to the case of finite loci with multiple substitutions permitted per locus and present an algorithmic solution (rather than a closed form solution) for how L, the hybrid load, should increase with K, the number of substitutions.

Let D be the total number of diverged loci between two populations. The number of potentially incompatible pairs is D(D − 1)/2. In the case that the number of loci,, is finite, D can increase to at most, and it will grow more slowly than K (the number of substitutions) due to the chance of repeated substitutions at the same locus. Modifying Equation 1, we derive L = 1 − [r with macron]D(D − 1)/2as the hybrid load when D loci have diverged. We numerically compute D, and hence L, for the case of finite, as follows. We create an array of integers, all initially set to zero. When the first random substitution occurs at a particular locus, the corresponding array element is set permanently to one: the locus has irreversibly diverged. At any point, we can sum the elements of the array to obtain the total number of diverged loci, D. (This approach neglects the possibility of back mutation, as did Orr’s original model. If back mutation, or alternatively, parallel evolution, were permitted, then D could decrease,) We also introduce an additional complication by partitioning the loci into “incompatibility groups” such that only pairs within a group can be incompatible. Let G be the number of groups, assumed to be of equal size. For example, for the case of Λ =1,000 and G=10, we take 10 groups of 100 loci each, combining their fitness penalties multiplicatively. In Figure 7, we plot L versus K for [r with macron] = 0.001; on the left is G=1 (a single group), and on the right is G=10. Ten replicate runs were done for each of Λ=1,000,000, Λ =1,000, Λ =100, Λ =30, and Λ =10; hence there are ten (stochastically overlapping) lines plotted for each of the first five colors in each plot. For comparison, Orr’s infinitely many loci case, Equation 1, is shown as f(K); and our correction for loci grouped into G “incompatibility groups” in the limit of Λ → ∞, is shown as g(K) = 1 − (e−(K/G)2[r with macron]/2)G For G=1, f(K) = g(K). For both values of G, the plots for Λ =1,000,000 (red) are very near the infinitely many loci case, the black dashed line representing g(K).

Figure 7
Growth of hybrid load L with the number of substitutions K in the case of finite Λg, the number of loci. Left: one “incompatibility group”; right: ten groups. r = 0.001.

On the left of Figure 7 the plots for Λ = 100 (blue) approach complete speciation (L = 1.0) more slowly than Orr’s infinitely many loci case (f(K), black line), due to the possibility of multiple substitutions at the same locus. More dramatically, the plots for Λ = 30 (pink) and Λ = 10 (cyan) on the left of the Figure do not approach 1.0 at all, but plateau around values of 0.35, and 0.04, respectively; L cannot continue to grow after all loci have diverged. We compute that L will asymptotically approach L(Λ,G,[r with macron]) = 1 − ((1 − [r with macron])Λ/G(Λ/G−1)/2)G as K approaches infinity and all loci diverge. Thus, for example, L(Λ = 30,G = 1,[r with macron] = 0.0001) ≈ 0.3529, matching the pink lines on the left of the figure. On the right side of the figure, for example, L(Λ = 100,G = 10,[r with macron] = 0.0001) ≈ 0.3625 matches the blue lines.

Clearly, finitely many loci can strongly limit the growth of incompatibility, as on the left of Figure 7; this is also reflected in our gene networks simulations, below. In addition, grouping of loci in “incompatibility groups” can limit incompatibility, or substantially delay its onset, even for high values of, as on the right of Figure 7. Additional volatility, and non-monotonicity of L, can be generated if other distributions are used for r (see (Orr 1995)); however, to model these cases, the individual values of r drawn for each incompatible pair must be tracked (not shown).

Reducing the Number of Interacting Loci Limits Incompatibility

Returning to our gene networks simulations, we indeed find that when, the number of loci, is reduced, incompatibility is restricted. Reduction of is effected in our simulations by reducing c. (Recall that a fraction 1-c of the matrix elements of the founding individuals are set to zero, and that mutation does not alter zero-valued matrix elements.) In Figure 8, we plot I against K2 for the case of c=0.1 with a=10 and a real-valued target. For this case, parallel genetic evolution did not occur, since the final element values between the two populations varied significantly (see online supplementary material). Despite the lack of parallel evolution, however, all pairs of populations maintain incompatibilities of less than one percent at the end of the run. Compare Figure 8 with the right side of Figure 6, for which c=1.0 and incompatibilities of 30–80% were reached. We may conclude that in this case incompatibility was limited not by parallel evolution, but because the number of interacting loci was small. Therefore, in general, either mechanism may prevent incompatibility from increasing.

Discussion

Another model of incompatibility

Nei et al. (1983) performed computer simulations of hybrid incompatibility in one- and two-locus multi-allele models. In the one-locus model, multiple alleles (labeled A−2, A−1, A0, A1, A2, etc.) are permitted at the locus. Homozygotes are viable, as are heterozygotes whose two alleles are only one step apart (i.e., an Ai Ai−1 heterozygote or Ai Ai+1 heterozygote are viable), but any other heterozygote is inviable. This is not the standard formulation of BDM incompatibility, which considers incompatibility due to epistasis between (at least) two loci. Nei et al.’s two-locus case considers two independent loci of the same description, with no epistatic interaction between them. In both the one- and two-locus cases, the low fitness of heterozygotes more than one step apart greatly retards the fixation of novel alleles, producing a fundamentally different dynamic from the independent fixations of Orr’s model. Nei et al. also observed an apparent delay, then a rapid acceleration of hybrid incompatibility, qualitatively similar to Orr’s “snowball.” However, in the Nei model, this dynamic is due to a delay and then rapid fixation of a novel allele, rather than to a growing number of potential incompatibilities, as in the Orr case.

Dynamics of incompatibility on networks of evolutionary pathways

The “holey landscapes” metaphor of Gavrilets (2004) describes high-dimensional genotype spaces as pervaded by “fitness ridges,” nearly-neutral networks on which populations may drift, constrained to remain on a ridge by stabilizing selection alone. Gavrilets’ original formulation did not include sexual recombination: the spatial structure of the landscapes – the allowable stepwise movements of genotypes on them – was defined in terms of mutational transitions only. Therefore, whereas an asexual population would remain confined to a “connected component” of genotypes by the action of mutation, introducing sexual recombination would break this confinement. Gavrilets extended this model to sexual populations by neglecting polymorphism (Gavrilets and Gravner 1997). In this case, a population can be represented by its single extant genotype, and transitions from one state to the next represent allele fixations in the population. This assumes only one mutation in the entire population for the duration of each fixation event, to ensure that the population remains confined to its current connected component.

We define the network of “evolutionary pathways” that population replicates may traverse as the union of the neutral networks spanning genotype (monomorphic population state) space at various fitness levels, together with the rare adaptive transitions connecting them. The evolutionary outcome depends stochastically upon the topology of this network. In this context, we may reinterpret Orr’s (1995) result as follows. Although it is true that two populations proceeding along divergent pathways will tend to have increasing genetic distance between them (and accelerating hybrid incompatibility, by Equation 2), genetic distance nevertheless may decrease for two reasons: 1) under neutrality, a population need not move in a consistent direction along a pathway; and 2) parallel genetic evolution may occur, with two pathways converging to an identical genotype by different routes (e.g., the same fixations but in different order), or asynchronously. A set of pathways diverging from the same point may in general have different probabilities of being evolutionarily “chosen” by selection or drift. Orr (2005) discusses how this increases the probability of parallel evolution.

Because genetic distance may decrease, we replace K, the number of substitutions (which can never decrease), with D, the genetic distance measured by the number of divergent alleles (which may decrease in certain circumstances). Equation 1 then becomes:

L1j=1D(1r¯)j11eD2r¯/2
(8)

In addition, [r with macron] should not be assumed to be constant between different pairs of pathways. Orr’s result, reinterpreted via Equation 8, does not predict an inevitable increase of hybrid load, although, if two populations do diverge genetically, increasing (indeed, initially accelerating) hybrid load is expected. However, the rate of this increase will depend on the particular [r with macron] between the two pathways along which the two populations are moving.

The evolution of gene expression

A classic paper by King and Wilson (1975) is widely cited for popularizing the notion that much phenotype divergence between species may be due to changes in the regulation of gene expression, as distinct from changes in protein-coding DNA. Although the role of structural mutations clearly remains important in evolution (Hoekstra and Coyne 2007), it has been shown in yeast (Borneman et al. 2007) and in primates (Gilad et al. 2006b) that gene expression levels diverge more strongly between closely related species than their nucleotide sequences. For example, between humans and chimpanzees, it was found that 10% of genes differed in expression in at least one brain region (Khaitovich et al. 2004).

The role of natural selection in the evolution of gene regulation is hotly debated. Many researchers favor the idea that the evolution of gene expression is primarily driven by selection (Nuzhdin et al. 2004; Gilad et al. 2006a; Gilad et al. 2006b; Haygood et al. 2007; Mustonen and Lassig 2007; Hutter et al. 2008), with some studies advocating stabilizing selection or directional selection as more “important”. However, a significant minority urge that the role of neutral evolution must not be neglected (Lynch 2007; Fay and Wittkopp 2008). It is in the context of this debate that we present our results on the simulated evolution of gene networks in a constant environment.

We have focused on constant environments in this article; however, our gene networks also show rapid growth of hybrid incompatibility in the presence of 1) divergent selection (selection toward different optima in the two populations), and 2) directional, nondivergent selection (selection toward the same optimum in both populations). That incompatibility would evolve under divergent selection is intuitive: if two populations evolve toward two different optima, it is not surprising that a hybrid might be less fit in both environments. That incompatibility increases under directional, nondivergent selection is less intuitive; however, it can be understood if one considers that there are multiple genetic means to the same end: this is “divergent DSE despite parallel phenotypic evolution”, as we described above, and that selected phenotypic change may drive a faster pace of genetic change than would proceed by drift.

Although speciation is rapid under both divergent and nondivergent directional selection, we have here focused on constant environments, demonstrating several dynamics, including strong growth of incompatibility, and waxing and waning incompatibility. Our conclusion is that the role of neutral evolution in the growth of incompatibility in biological gene regulatory networks, and hence in speciation, should not be dismissed out of hand.

Conclusions

We have extended Orr’s model in two ways to account for the complex dynamics of incompatibility in the gene networks model: 1) by considering genetic distance, (which may decrease by back mutation or parallel genetic evolution) rather than number of substitutions (which does not decrease), and 2) by considering a finite, rather than an infinite, number of interacting loci.

Importantly, we have shown that directional selection is not required to produce hybrid incompatibility in a model inspired by biological gene networks. Rather, the growth of hybrid incompatibility may proceed by any of several dynamics in a constant environment, including asynchronous parallel evolution, and developmental system evolution with no change in phenotype. The important point is that many processes can cause genetic change, some of which can increase hybrid incompatibility. Divergent directional selection clearly will promote phenotypic divergence between two populations, which necessarily implies genetic divergence between them. Nondivergent directional selection will promote phenotypic change, but not necessarily phenotypic divergence; hence it may or may not result in genetic divergence. Likewise, neutral phenotypic change may or may not result in genetic divergence. Our models suggest that two mechanisms, 1) phenotypically neutral changes in gene regulation, and 2) asynchronous parallel evolution due to rarity of adaptive mutation, may also lead to speciation in biological populations.

Supplementary Material

supplementary material

Acknowledgments

The authors thank Michael Desai, Laurent Lehmann, Joanna Masel, and Sohini Ramachandran for important discussions and suggestions. We are also indebted to Adam Porter (in particular for the suggestion that incompatibility could be limited in the finite loci case) and an anonymous reviewer.

Literature Cited

  • Bateson W. Heredity and variation in modern lights. In: Seward AC, editor. Darwin and modern science: essays in commemoration of the centenary of the birth of Charles Darwin and of the fiftieth anniversary of the publication of the Origin of species. Cambridge: Cambridge University Press; 1909. pp. 85–101.
  • Bergman A, Feldman WM. On the population genetics of punctuation. In: Crutchfield JP, editor. Evolutionary dynamics: Exploring the interplay of selection, accident, neutrality, and function. Oxford University Press; 2003. pp. 81–100.
  • Borneman AR, Gianoulis TA, Zhang ZDD, Yu HY, Rozowsky J, Seringhaus MR, Wang LY, et al. Divergence of transcription factor binding sites across related yeast species. Science. 2007;317:815. [PubMed]
  • Cooper TF, Rozen DE, Lenski RE. Parallel changes in gene expression after 20,000 generations of evolution in Escherichia coli. Proceedings of the National Academy of Sciences of the United States of America. 2003;100:1072. [PubMed]
  • Cooper VS, Lenski RE. The population genetics of ecological specialization in evolving Escherichia coli populations. Nature. 2000;407:736. [PubMed]
  • Coyne JA, Orr HA. Patterns of speciation in Drosophila. Evolution. 1989;43:362.
  • Coyne JA, Orr HA. Patterns of speciation in Drosophila revisited. Evolution. 1997;51:295.
  • Darwin C. The Origin of Species by Means of Natural Selection. London: John Murray; 1859.
  • Dobzhansky T. Studies on hybrid sterility. II. Localization of sterility factors in Drosophila pseudoobscura hybrids. Genetics. 1936;21:113–135. [PubMed]
  • Fay JC, Wittkopp PJ. Evaluating the role of natural selection in the evolution of gene regulation. Heredity. 2008;100:191. [PubMed]
  • Gavrilets S. Fitness landscapes and the origin of species. Princeton, N. J: Princeton University Press; 2004.
  • Gavrilets S, Gravner J. Percolation on the fitness hypercube and the evolution of reproductive isolation. Journal of Theoretical Biology. 1997;184:51. [PubMed]
  • Gilad Y, Oshlack A, Rifkin SA. Natural selection on gene expression. Trends in genetics. 2006a;22:456. [PubMed]
  • Gilad Y, Oshlack A, Smyth GK, Speed TP, White KP. Expression profiling in primates reveals a rapid evolution of human transcription factors. Nature. 2006b;440:242. [PubMed]
  • Haygood R, Fedrigo O, Hanson B, Yokoyama KD, Awray G. Promoter regions of many neural- and nutrition-related genes have experienced positive selection during human evolution. Nature genetics. 2007;39:1140. [PubMed]
  • Hoekstra HE, Coyne JA. The locus of evolution: Evo devo and the genetics of adaptation. Evolution. 2007;61:995. [PubMed]
  • Hutter S, Saminadin-Peter SS, Stephan W, Parsch J. Gene expression variation in African and European populations of Drosophila melanogaster. Genome biology. 2008;9:R12. [PMC free article] [PubMed]
  • Johnson NA, Porter AH. Rapid speciation via parallel, directional selection on regulatory genetic pathways. Journal of Theoretical Biology. 2000;205:527. [PubMed]
  • Johnson NA, Porter AH. Evolution of branched regulatory genetic pathways: directional selection on pleiotropic loci accelerates developmental system drift. Genetica. 2007;129:57. [PubMed]
  • Khaitovich P, Muetzel B, She X, Lachmann M, Hellmann I, Dietzsch J, Steigele S, et al. Regional Patterns of Gene Expression in Human and Chimpanzee Brains. Genome Res. 2004;14:1462–1473. [PubMed]
  • King MC, Wilson AC. Evolution at two levels in humans and chimpanzees. Science. 1975;188:107. [PubMed]
  • Lijtmaer DA, Mahler B, Tubaro PL. Hybridization and postzygotic isolation patterns in pigeons and doves. Evolution. 2003;57:1411. [PubMed]
  • Lynch M. The frailty of adaptive hypotheses for the origins of organismal complexity. Proceedings of the National Academy of Sciences of the United States of America. 2007;104:8597. [PubMed]
  • Mendelson TC. Sexual isolation evolves faster than hybrid inviability in a diverse and sexually dimorphic genus of fish (Percidae: Etheostoma) Evolution. 2003;57:317. [PubMed]
  • Moyle LC, Olson MS, Tiffin P. Patterns of reproductive isolation in three angiosperm genera. Evolution. 2004;58:1195. [PubMed]
  • Muller HJ. Isolating mechanisms, evolution, and temperature. Biological Symposia. 1942;6:71–125.
  • Mustonen V, Lassig M. Adaptations to fluctuating selection in Drosophila. Proceedings of the National Academy of Sciences of the United States of America. 2007;104:2277. [PubMed]
  • Nei M, Maruyama T, Wu CI. Models of evolution of reproductive isolation. Genetics. 1983;103:557. [PubMed]
  • Nuzhdin SV, Wayne ML, Harmon KL, McIntyre LM. Common Pattern of Evolution of Gene Expression Level and Protein Sequence in Drosophila. Mol Biol Evol. 2004;21:1308–1317. [PubMed]
  • Orr HA. The population genetics of speciation: The evolution of hybrid incompatibilities. Genetics. 1995;139:1805–1813. [PubMed]
  • Orr HA. Dobzhansky, Bateson, and the genetics of speciation. Genetics. 1996;144:1331. [PubMed]
  • Orr HA. The probability of parallel evolution. Evolution. 2005;59:216. [PubMed]
  • Orr HA, Turelli M. The evolution of postzygotic isolation: Accumulating Dobzhansky-Muller incompatibilities. Evolution. 2001;55:1085. [PubMed]
  • Philippe N, Crozat E, Lenski RE, Schneider D. Evolution of global regulatory networks during a long-term experiment with Escherichia coli. Bioessays. 2007;29:846–860. [PubMed]
  • Porter AH, Johnson NA. Speciation despite gene flow when developmental pathways evolve. Evolution. 2002;56:2103. [PubMed]
  • Presgraves DC. Patterns of postzygotic isolation in Lepidoptera. Evolution. 2002;56:1168. [PubMed]
  • Sasa MM, Chippindale PT, Johnson NA. Patterns of postzygotic isolation in frogs. Evolution. 1998;52:1811.
  • Siegal ML, Bergman A. Waddington’s canalization revisited: Developmental stability and evolution. Proc Natl Acad Sci USA. 2002;99:10528–10532. [PubMed]
  • True JR, Haag ES. Developmental system drift and flexibility in evolutionary trajectories. Evolution & development. 2001;3:109. [PubMed]
  • van Nimwegen E, Crutchfield JP, Huynen M. Neutral evolution of mutational robustness. Proc Natl Acad Sci USA. 1999;96:9716–9720. [PubMed]
  • Wagner A. Does evolutionary plasticity evolve? Evolution. 1996;50:1008–1023.