extfieldsofstudy:
Biology
Computer Science
Medicine

text:
Neutral Models of Microbiome Evolution

There has been an explosion of research on host-associated microbial communities (i.e.,microbiomes). Much of this research has focused on surveys of microbial diversities across a variety of host species, including humans, with a view to understanding how these microbiomes are distributed across space and time, and how they correlate with host health, disease, phenotype, physiology and ecology. Fewer studies have focused on how these microbiomes may have evolved. In this paper, we develop an agent-based framework to study the dynamics of microbiome evolution. Our framework incorporates neutral models of how hosts acquire their microbiomes, and how the environmental microbial community that is available to the hosts is assembled. Most importantly, our framework also incorporates a Wright-Fisher genealogical model of hosts, so that the dynamics of microbiome evolution is studied on an evolutionary timescale. Our results indicate that the extent of parental contribution to microbial availability from one generation to the next significantly impacts the diversity of microbiomes: the greater the parental contribution, the less diverse the microbiomes. In contrast, even when there is only a very small contribution from a constant environmental pool, microbial communities can remain highly diverse. Finally, we show that our models may be used to construct hypotheses about the types of processes that operate to assemble microbiomes over evolutionary time.

Introduction
Microbial communities associated with animals and plants (i.e., microbiomes) are implicated in the day-to-day functioning of their hosts in a variety of ways; microbes provide hosts with access to nutrients [1][2][3][4], they protect against pathogens [5,6], confer drought resistance [7], mediate hosts' social interactions [8], and modulate host behavior [9][10][11]. Not surprisingly, many studies have related human health to the human microbiome. Researchers have proposed that perturbations in human microbiome composition may be associated with a range of disorders, including inflammatory bowel disease [12][13][14], anxiety and depression [15], obesity [16,17], autism [18,19], allergic responses [20] and respiratory ailments [21][22][23]. There is a growing number of large-scale projects underway to characterize and analyze the collection of microbes linked to human health and disease using advanced sequencing methods and technologies, including the Human Microbiome Project [24] and the European Metagenomics of the Human Intestinal Tract [25].
Many microbiome studies focus on descriptions of microbial dynamics within hosts over short periods [26][27][28], or on the composition and diversity of microbial communities amongst hosts with different phenotypes or under different treatment regimes [17,29], but there is increasing recognition that host-microbe dynamics need to be studied within the context of a unified ecological or evolutionary framework. In a recent paper, Costello et al [30] listed four ecological processes that mediate the diversity of human microbiota: environmental selection, whereby the host environment favors the presence and persistence of certain microbial taxa; historical contingency, in which differences in the timing and order of microbial colonization leads to differences in succession and climax communities; random sampling, where stochastic factors influence community assemblages; and dispersal limitation, where the availability of microbial taxa is restricted by the local structure of host communities and environments. Yeoman et al [31] discuss many of these same factors, and frame the challenge of understanding microbial diversity through the lens of evolution, identifying selection induced by competition between individuals of the same or different microbial species, and host influences including phylogeny, as important drivers of host-microbial variation.
In this paper, we take the first steps in developing an explicit framework for modelling the evolutionary and ecological dynamics of microbial communities within a population of hosts. The models we use are related to some of those that metacommunity theorists work with [32,33], specifically, neutral models of ecology and biodiversity [34], but there is one significant difference. Standard spatial models of community assembly, biogeography and biodiversity assume that areas available for colonization remain indefinitely and are static. These models do not apply to the evolution of microbiomes because the "spaces" these microbial communities colonize-that is, the hosts-share an evolutionary history characterized by lineages that persist or go extinct. For this reason, the pairing of host genealogy with microbiome assembly may lead to different patterns of biodiversity than those expected under existing models of microbial community assembly and metacommunity theory.
We develop an agent-based framework, similar to that used recently by Hellweger et al [35], to capture the emergent patterns of microbial diversity over many host generations. In contrast to Hellweger et al, we do not model mutations and speciation in microbial lineages but focus instead on how ecological and evolutionary sampling processes affect the standing variation of microbes in the environment and hosts. The neutral processes that are encapsulated in our framework are based on those that have been proposed by others [30,31] with the added dimensions of an evolutionary timescale and a genealogy of hosts, both of which contribute to shaping the eventual composition and variation of microbial communities within and between hosts. Our models do not assume that hosts differ in their reproductive success as a consequence of their microbial communities; nor do our models assume differences amongst microbes in their propensities to persist in the host or the environment, or in their abilities to disperse. Consequently, whereas our models do not capture the complexity of processes that are likely to mediate the evolution of microbiomes, they serve as minimalist null models against which empirical patterns may be compared.
Our framework is a generic one, and we have not developed it specifically for any one species of host, or indeed, any one systemic compartment within a host. Instead, our models incorporate simple host mechanisms for microbial acquisition, and we explore the effects of parental inheritance of the microbiome, and microbial recruitment from the environment. Our results indicate that parental inheritance tends to reduce microbial diversity and increase homogeneity within hosts, while ongoing environmental acquisition works to maintain microbial heterogeneity within hosts. Interestingly, we observe a non-linear relationship between the degree of parental inheritance and between-host differences in microbiome composition. Finally, we show how these neutral patterns allow us to make predictions about the processes that are important in shaping microbiome diversity when applied to empirical data.

The models
Our framework applies to a population of hosts and an available pool of microbial colonists. As a first step, we assume that hosts do not exert any preferences on the microbial taxa they acquire. Similarly, we assume that microbes do not interfere with host reproductive capacity, or the survivorship and reproductive success of other microbes in the community. As will become clear, the only indirect effects that influence microbial recruitment and persistence from one generation of hosts to the next are competition for space within hosts and the relative abundance of microbial taxa. Simply put, we assume that the ecological and evolutionary processes that operate on hosts and their microbiomes are neutral; in this regard, our framework is analogous to neutral theories in evolutionary biology [36,37], ecology and biodiversity [34]. We expand on this analogy later, but for now, we note that neutral theories provide parsimonious accounts of the types of patterns that can emerge in complex systems, they serve as null models for statistical hypothesis tests, and they provide platforms upon which we may construct more elaborate representations of these same systems [38].
In this framework, hosts reproduce asexually in discrete generations, following a neutral Wright-Fisher process [39,40], where each individual in a succeeding generation chooses a parent randomly from the preceding generation. Hence, with a population of hosts of constant size N, all asexual individuals will share a common ancestor after 2N generations, on average. In our models, asexual reproduction is a computational convenience, and can be replaced with sexual reproduction without changing the essential patterns that we observe.
We model how hosts acquire their microbiomes in three ways (Fig 1). First, under a strict "parental-acquisition" (PA) process, all hosts acquire their microbial communities directly from their parents. Second, with strict "environmental-acquisition" (EA), hosts acquire their microbiomes solely from the environment. Between these two extremes, we also allow a third "mixed-acquisition" (MA x ) process, whereby hosts acquire some percentage, x%, of their microbiomes from their parents and (100-x)% from the environment. MA 0 is exactly equivalent to EA, and MA 100 to PA; as such, EA and PA designate boundary conditions of the ecological processes that mediate microbial acquisition in hosts.
It is worth pausing at this point to clarify what we mean when we say that hosts acquire their microbiomes from their "parents" or their "environments". Our models do not explicitly take account of the life events-illness, infections, changes in environments or diets-of each host within a generation, nor does it consider microbial fluxes within the lifespan of each host. Instead, the microbial composition of each host is essentially measured as an aggregate over the single generation that the host exists. Consequently, when we quantify the percentage of microbes from parents and environment using, say, MA 10 , we mean that over the life of the host, 90% of its microbes come from the environment and 10% from its parent. In our models, it is possible that the parental contribution happened in the first 10% of the host's life, or it may be that over the entire lifespan of the host, there was an ongoing contribution by the parent that amounted to 10% of the microbial composition.
Since we allow hosts to recruit microbes from an "environment", we need to define how the microbial content of this environment is constituted. In simulations, we characterize microbial composition using a distribution of taxa' relative abundances. We propose three processes that determine the composition of the pool of microbes available for recruitment. First, we assume that the environment has a microbial composition that remains fixed over time. For the "fixed environment" (FE), all taxa are present in the environment throughout the simulation, and are available to every generation of hosts. The second process we propose involves a changing environmental microbial profile, whereby the relative abundance of each microbial taxon available to the hosts in a given generation, is an aggregate of their abundances from all hosts of the preceding generation. Under this "pooled-environment" (PE), microbial composition is reflection of what was present in the parents of the current generation of hosts. A third, intermediate, process is a combination of the previous two "environments": the environmental microbial pool available for recruitment contains a percentage, y%, from the parental pool of microbes, and (100-y)% of microbes from the fixed environment. Under this "mixed environment" (ME y ), the proportion of contribution from host microbiomes is given by y. As with our acquisition models, ME 0 and ME 100 are equivalent to the boundaries FE and PE, respectively.
Our framework allows us to combine different host-acquisition processes with different ways of constructing the pool of available microbes in the environment. Conceptually, each of these combinations is a particular neutral model, capturing some of the elements previously discussed in the literature. For instance, PA or MA x incorporate the phylogenetic dependencies that Yeoman et al [31] discuss, and EA x PE is equivalent to what Costello et al [30] call dispersal limitation, whereby the local host community influences microbial composition. It is worth noting that the combinations PA x (FE, ME y, PE)-read as "PA in combination with FE, with ME y or with PE"-will give identical results. This is because, in all cases, the environment contributes nothing to host microbial content (see the first row in Fig 1).
We model the construction of the microbial community in each host by competitive random sampling with replacement. Under this process, each host allows only a fixed and limited number of microbes to populate its microbiome. If microbial acquisition occurs under EA, each host samples randomly from the available pool of taxa according to the relative abundance of each taxon in the environment. In the case of MA x , x% of microbes are selected from the parent and (100-x)% from the environment. If hosts acquire their microbial taxa under PA, then all microbes are inherited from the hosts' parents, although the relative abundance of each taxon fluctuates multinomially. By constructing microbial communities in this way, we allow stochastic factors and indirect competition to modify taxon composition within and between hosts, as proposed by Costello et al [30].
By simulating combinations of PA, MA x and EA against FE, ME y and PE forward in time over many host generations and over a range of conditions, we are able to recover data on the behavior of individual microbial taxa, as well as a variety of summary statistics, including the expected time it takes individual taxa to invade all hosts or go extinct in the host population, and the trajectories of microbial taxonomic richness (measured simply as the total number of microbial taxa) and microbial taxonomic evenness (measuring the similarity in the frequency of each taxon), microbial diversity within hosts (α-diversity), inter-host variation in microbial composition (β-diversity) and the aggregate microbial diversity from all hosts in the population (γ-diversity). Here, we report only on the latter three measures of diversity.

Relative abundances of microbial taxa in the host population
Microbial diversity within the host population is a function of the proportion of microbes that parents contribute directly to offspring and the proportion they contribute to the environment. Fig 2 illustrates how population-level taxon abundances change under various combinations of these proportions. In our simulations, the distributions of taxon abundances under high levels of parental contributions are skewed, and may be approximated by commonly-applied distributions, including the log-normal distribution and the Dirichlet multinomial (DM) distribution [41] (Fig 3; the DM distribution has the advantage of allowing α-, βand γ-diversities to be simulated-see S2 Fig). The ability to recover skewed abundance distributions is interesting, because we begin our simulations with a uniform distribution of microbial taxa, and we retain this uniform distribution in the fixed environment throughout the evolutionary history of the host population. Consequently, the emergence of dominant and rare taxa is a consequence of repeated parental contributions either directly to the next generation of hosts or indirectly to the environment.  (50) or 50% parental acquisition); (b) distribution of relative taxon abundances with different parental contributions to the environment (ME), when the proportion of direct parental acquisition of microbes is high (i.e., MA(98.4) or 98.4% parental acquisition); (c) distribution of relative taxon abundances with different proportions of microbes directly acquired from parents (PA), when the proportion of parental microbial contribution to the environmental pool is low (i.e., ME (50) or 50% parental contribution); (d) distribution of relative taxon abundances with different proportions of microbes directly acquired from parents (PA), when the proportion of parental microbial contribution to the environmental pool is high (i.e., ME(99.9) or 99.9% parental contribution). In fact, all simulations in which there was complete parental acquisition of microbes (i.e., PA) resulted in the loss of all but one microbial taxon in the host population. Similarly, when the environment was reconstituted each generation exclusively with microbes from the parents (i.e., PE), the same pattern was observed with only a single microbial taxon remaining. These result are consistent with predictions made under neutral models of community ecology [42], and highlight the strong depressive effect of parental transmission, either directly from parent to offspring or via parental contributions to a local pool of microbes, on population-level microbiome diversity.
With EA x FE, microbes are obtained randomly from a fixed environment that persists over the evolutionary history of the hosts; unsurprisingly, the host population retains all microbes found in the environment. Interestingly, when microbes are obtained both from parents and a fixed environment (MA x FE), we still see the persistence of all or almost all microbes in the host population (see Fig 2A and 2B, first column of each bar chart; ME(0) is equivalent to a fixed environment with no microbial contributions from parents). This is true even when the proportion of microbial taxa that an individual host acquires from the fixed environment at each generation is very small, on the order of 0.001. Therefore, a very small contribution from a constant environmental source of microbes is sufficient to retain high levels of microbial diversity in the host population.

Measures of diversity
Microbial diversities are frequently measured in three ways: α-diversity, β-diversity, and γdiversity. Our simulations indicate that all three measures depend on the percentage of parental contribution to offspring microbiomes and the composition of the environmental microbial pool (see S1-S3 Tables for simulation means and standard deviations). Under our neutral model, in which the absence of host sub-population structure means that all hosts sample their microbes from the same environment, αand γ-diversities remain high, and β-diversity remains low, for a large part of the range of direct or indirect parental contributions (i.e., to offspring or to the environment, respectively). Nonetheless, at high values of parental contributions, there are discernible differences in diversities, and we have also focused our simulations in these areas (Fig 4; see S4-S6 Tables for simulation means and standard deviations). In general, α-diversity (average diversity within hosts) and γ-diversity (overall diversity within the entire population of hosts) increase as we increase the fixed environmental contribution because a fixed environment helps maintain a uniform distribution of taxon abundances and delays the loss of microbial taxa during evolution. Conversely, when hosts acquire increasing proportions of their microbiomes from their parents directly, or indirectly from a pooled environment, the variation of taxon abundance increases and taxon richness tends to decrease, thus lowering both αand γ-diversities (Fig 4A and 4B; S1, S2, S4 and S5 Tables). Inter-host variation in microbial composition, or β-diversity, also depends on the degree of parental inheritance, and the ratio of fixed-to-pooled environmental components (Fig 4C; S3 and S6 Tables).
Under the combination of PA x (FE, ME or PE), β-diversity tends to zero, because all hosts descend from a single common ancestor and, as noted above, only a single microbial taxon remains in all hosts. When we have the parental microbiome as the only source of microbiomes in the next generation, ultimately, all lineages will have acquired their microbiomes from the most recent common ancestor (MRCA) of the population of hosts. Additionally, from one generation to the next, stochastic sampling of microbes over evolutionary time will result in the loss of all but one microbial taxon. Interestingly, with a high percentage of environmental acquisition, β-diversity is also relatively low, because all hosts acquire a large proportion of their microbial taxa from the same environmental pool, and consequently, will tend to acquire the same set of taxa. As noted above, the highest β-diversity occurs in a relatively narrow range of values of pure parental acquisition (between 87-99% of direct parental transmission; S6 Table).
If we focus on the relationship between β-diversity and the environmental pool, we see that its behavior is similar to that of αand γ-diversities: it decreases as we increase the pooled environmental contribution to offspring microbiome. This is because a pooled environment, with contributions from the parental generation, tends to give rise to a non-uniform distribution of microbial taxa. As the degree of parental contribution increases, the environmental community will be dominated by few highly abundant species which are likely shared by most or all hosts within the population, accounting for high between-host similiarity in microbial composition (Fig 4C; S3 and S6 Tables).
It is important to note that our simulations have not been performed with inference or prediction in mind: the number of hosts, the number of microbes, and the number of taxa in our simulations are not necessarily equivalent to those of real-world microbial communities and their hosts, nor have we necessarily chosen the appropriate diversity indices or taxonomic resolution to optimize prediction/inference. Nonetheless, it is helpful to examine how the simulated values of diversity compare to empirical observations, and what these comparisons might tell us about the evolutionary processes that are acting on microbiomes. As an example, we used genus-level taxonomic data from the NIH Human Microbiome Project (HMP) [43], specifically, a table of relative abundance found in different compartments of the human body (http://www.hmpdacc.org/HMSMCP/; see S1 Dataset). Several large samples from the anterior nares, vaginal posterior fornix, stool, buccal mucosa, tongue dorsum and supragingival plaque were chosen to calculate α-, βand γ-diversities on genus level ( Table 1).
Values of αand γ-diversity obtained from all sampled sites of the human microbiome are low, in comparison to most of the values we obtained in our simulations. In fact, human microbiome diversities are generally lower than those of other non-human primates [44,45]. If we compare the empirical diversities to those obtained in our simulations (S1-S6 Tables), we would have to posit very high parental contributions, both direct and indirect (>90%), to account for the αand γ-diversities across all human body sites. In contrast, values of β-diversity appear to provide a little more discrimination amongst body sites: the site with the lowest β-diversity is the vaginal posterior fornix, and its value is consistent with a very low degree of direct parental contribution in our simulations (approximately between 0-15%). The β-diversities at other sites appear to suggest higher levels of parental contribution (again, >90%). In the next section, we discuss the implications of these results, as they relate to human microbiome evolution and how the neutral model may be used to construct hypotheses about relevant evolutionary and ecological processes.

Discussion
In this paper, we introduce a simple and flexible framework to model the evolution of microbiomes within a population of hosts, which takes account of different modes of microbiome The αand βdiversities are mean values and standard errors averaged over all samples or over all pairs of samples. 2 The γdiversity is a single value because it is calculated with a large sample that is pooled from all the samples within the same body site; hence no replicates or ranges are available.
doi:10.1371/journal.pcbi.1004365.t001 acquisition and environmental microbial composition. Under our neutral model, microbiome composition is affected by sampling effects. Stochastic changes in microbial abundances may affect the persistence of microbial taxa in the microbiome over one or a few generations (i.e., ecological drift), or over many generations. The latter may occur because host lineages die out; when this happens, changes in microbial abundance across the whole population of hosts are essentially equivalent to changes in allele frequencies (i.e., genetic drift). The constitution of the microbial community in the environment also plays a considerable role in determining the ultimate fate of microbial taxa within host microbiomes. With a fixed environment, when there is a constant pool of the same microbial taxa from one generation of hosts to the next, microbial taxa never go extinct from the host population as long as hosts obtain some fraction of their microbiome from the environment. This is true, even when that fraction that the environment contributes to each host's microbiome is very small (e.g., 0.1%). In contrast, when the environmental composition of microbes reflects the microbial content of the hosts in previous generations (i.e., PE, the "pooled" environment in our model), microbial diversity of the environment shrinks, as does the diversity of host microbiomes. Therefore, the extent to which parents contribute to the microbiomes of their offspring (either directly or through their contributions to the pooled environment) plays a crucial role in shaping microbiome diversity and constitution. In our simulations, values of αand γdiversities are at their lowest when parental contribution to the microbiome is high. Inevitably, microbial taxa are lost from the population as host lineages are lost.
Thus, under our neutral model, it is possible to recover skewed microbial abundance distributions reminiscent of those obtained with real data, despite a fixed environmental component that remains uniform and constant throughout our simulations. Increasing skewness-essentially, decreasing eveness-is obtained as we increase the degree of parental inheritance. Of course, we don't claim any deep insight here: no one should be surprised that we are able to recover skewed abundance distributions with our models, because there is a large body of literature on the mechanisms-both neutral or otherwise-that may lead to the emergence of skewed abundance distributions (see [46] for an excellent synthesis). Our results reinforce what others [34,47,48] have found, by adding yet another neutral mechanism to account for the emergence of skewed abundance distributions.
Our framework includes sampling effects on an undivided host population, which evolves under a Wright-Fisher process. Consequently, our models have some points of similarity with those that have been developed in population genetics. For instance, Orive et al [49] analyze the evolutionary dynamics of endosymbionts using a discrete-time Moran population genetic model. In their model, endosymbionts are acquired either vertically, passed on from parent to offspring, or horizontally from the environment. This corresponds to our MA x FE model and, in agreement with our results, Orive et al find that increasing the environmental contribution of endosymbionts to host cells results in greater diversity within cells and less diversity between cells.
Our models do not include any mutational process or speciation acting on the microbes, as time moves forward. In reality, of course, microbes acquire mutations in their genomes at a rapid rate, but the measures of diversity we use in our analyses capture differences in taxonomic composition, not genetic diversity. In our models, it is implied that no cladogenetic events have occurred over the course of the simulations.
The models presented here provide an opportunity to construct hypotheses, and make qualitative predictions, about the patterns of diversity we can expect to find in different biological situations. For example, the effects of "pure" pooled versus fixed environments on microbiomes can be found in a comparison of social and solitary bees. Social bees exhibit behaviors that are likely to result in the transmission of microbes from a microbial pool within the colony [50]. In contrast, solitary bees acquire their microbiomes from the environment, through feeding or burrowing. Our model would predict that social bees would have lower α-diversity and lower taxonomic richness than solitary bees. This is consistent with the results obtained by Martinson et al [51] who surveyed the microbiomes of eusocial bee species Apis spp. and Bombus spp. and non-social bees (11 species) and wasps (3 species): they found depauperate microbiomes in social bees compared to non-social bees.
Whereas it is reassuring to obtain empirical corroboration for our models, arguably neutral models are most useful when real-world observations run counter to predicted outcomes. Falsification of neutral models provides a justification for augmenting these models to include additional processes that account for the phenomena under study. In this regard, our analysis of the Human Microbiome Project data is instructive. As we have noted above, our simulations should not be used for inference, and we should be cautious about reading too much into the comparisons between empirical and simulated patterns of diversity. Nonetheless, at least for some sites, i.e., the stool, the tongue dorsum, the supragingival plaque, the anterior nares and the buccal mucosa, the low empirical values of αand γdiversities appear to point consistently to a high level of parental inheritance when compared against values obtained in our simulations. There is evidence that the human microbiomes at various sites are seeded at birth by the mother, particularly if this birth is through the vaginal tract [52]. There is also reasonably strong evidence that families share microbes to a greater extent than unrelated individuals in a population [53], and at least in some human populations, mothers share more microbes in common with their offspring than with unrelated children [54]. It is not clear, based on the studies that have been done to date [55], whether the values of direct or indirect parental contribution we obtain when we compare empirical and simulated diversities are significantly higher than would be obtained in real populations, but we expect that the intuition of mosts microbial ecologists is that percentages of direct and pooled parental contributions > 90% are likely to be too high. Putting to one side the caveats about inference, we accept that while this intuition does not constitute evidence against the neutral model, it is likely to engender scepticism about the model's correctness. If it is, in fact, true that direct or indirect parental contributions to the next generation's microbiomes are not as high as our simulations suggest, how do we account for the apparent depression in αand γ-diversities, and elevation of β-diversities at these sites? One hypothesis that explains low αand γ-diversities, and high β-diversities, and does not require the action of non-neutral processes, is the existence of local host subpopulations. The existence of subpopulations of hosts, with limited immigration and sharing of microbes between subpopulations, is likely to give the appearance of high parental contribution from one generation to the next. Certainly, this is a plausible explanation for patterns of microbiome diversity in the oral cavity (i.e., the buccal mucosa, tongue dorsum and supragingival plaque) and stool samples, because of the likely influence of familial [53] or cultural dietary preferences/practices [56] or lifestyles [57] on these microbiomes. A similarly explanation may account for patterns of microbiome diversity of the anterior nares.
The vaginal posterior fornix presents an interesting contrast to the other body sites because the αand γ-diversities suggest high parental contributions (although they cannot distinguish between direct or indirect contributions), whereas β-diversity suggests a low direct parental contribution. This inconsistency may again cause us to reject the neutral model in favor of an alternative explanation, but in this case, subpopulation structure may play a minor role relative to selection for a vaginal microbial community that is common amongst hosts. Such a selective filter is likely a consequence of a complex suite of factors including host immune defences, hormonal cycling, pregnancy, and the presence of apparently beneficial microbial species (e.g., Lactobacillus spp.) [58]. This hypothesis explains both the high level of αand γ-diversity (i.e., a few abundant species with many rare species), and the low β-diversity.
For the human microbiome, neutral models have the potential to help identify additional processes that may account for patterns of diversity. As noted, of the two processes identified above-host subpopulations and selective filters-the former still remains part of an underlying neutral process, and a plausible extension to the neutral framework presented here. Rejection of a simple neutral model therefore allows us to identify incremental additions that may increase explanatory power.
Another example of empirical data that appears to contradict the expectations of our models is the comparison of microbiome diversities in high microbial abundance (HMA) and low microbial abundance (LMA) sponges [59]. HMA sponges have large numbers of associated microbes, in contrast to LMA sponges. Additionally, researchers have shown that microbial diversity in LMA sponges is lower than that of HMA sponges [60,61]. Based on our results, we would predict that there is a greater degree of vertical transmission in LMA sponges, but it turns out that this is not the case: Schmitt et al [62] have found that "vertical transmission, as a mechanism to obtain bacteria, seems to occur mainly in HMA sponges". Giles et al [60] propose two possible reasons to account for the low diversity in LMA sponges. First, there may be selective filters that permit only certain microbial taxa to colonize the sponges; second, the initial colonization event is stochastic, but serves to constrain or exclude successive colonizations. As with the human microbiome data, the sponge example is important because it does not rely a priori on non-neutral processes to account for the low diversity in LMA sponges; instead, selection (or other ecological and/or evolutionary processes) is invoked only after it is shown that vertical transmission in LMA sponges is unlikely, thus indicating that our neutral models are an inadequate explanation for the observed data.
Riffing on the theme that "Essentially, all models are wrong, but some models are useful" [63], Hubbell, writing about models in community ecology, says "Probably no ecologist in the world with even a modicum of field experience would seriously question the existence of niche differences among competing species on the same trophic level" [64]. But, he continues, "[Neutral theory] begins with the simplest possible hypothesis one can think of . . . and then adds complexity back into the theory only as absolutely required to obtain satisfactory agreement with the data". We agree with Hubbell: to paraphrase, given what we know about the interplay between hosts, their microbial communities, and the environment, we would hesitate to put money on the table and bet that many microbiomes have evolved under the simplest neutral models that we have constructed here. But we would be equally hesitant betting in favor of the null hypotheses evaluated in statistical tests of significance. The value of these hypotheses resides not in their rightness or wrongness but in their ability to protect against overconfidence in our favorite, more complex model. Whereas it is true that biological processes are frequently complex, Occam's Razor dictates that we construct as simple explanations (or models) as possible. In this way, we remain vigilant against the addition of unnecessary and unjustifiable complexity. Much as we do with statistical hypothesis tests, we accept stronger alternative explanations only when we are sufficiently confident that our neutral hypotheses are unlikely. This is not to say that neutral models only serve as strawmen; in molecular evolution, for instance, neutral models are frequently effective at explaining molecular variation [65]. And even in cases when the assumption of neutrality is questionable, the use of neutral models of substitution applied in molecular phylogenetics does not appear to jeopardize the accuracy of tree reconstruction [66]. Consequently, without taking account of the evolutionary processes of mutation, speciation, selection or recombination, or the ecological processes that operate in the context of spatial, environmental, and temporal heterogeneity, what we have developed is a framework on which we can begin to evaluate empirical patterns of diversity, and where necessary, add more elaborate ecological and evolutionary scenarios. We believe that even this simple framework, devoid as it is of all the embellishments afforded by evolution and ecology, can serve a useful purpose: it is a suitable staging ground on which we can construct null models of microbiome diversity in populations of hosts and it allows us to make strong, testable predictions.

Materials and Methods
Simulated host populations consisted of a fixed number of virtual host individuals (N = 500). Each host was allocated a virtual microbiome with a limited capacity or "slots" of microbes (n = 1000). The environmental pool consisted of 150 microbial taxa. Large number of hosts (N = 2000), microbes (n = 100000) per host and microbial taxa (m = 500) were also simulated with our neutral model, and similar patterns of diversity were observed (S1 Fig). The microbiomes of the initial generation of host individuals were seeded randomly, with bacteria sampled from a uniform distribution of taxon abundances. We used an initial uniform distribution of taxa because we wanted to ascertain whether the equilibrium distribution of abundances obtained at the conclusion of our simulations would recover patterns seen in natural microbiomes.
For each subsequent generation, the microbiome of each individual host was simulated by populating each of the available "slots" in the individual's microbiome by sampling microbial taxa with replacement (multinomial choice) from either the environment (with probability given (1-x)) or from the microbiome of a "parent" host individual (selected with uniform random probability from the population of the previous generation).
When sampling from a parental/environmental microbial community, the probability that the new host microbiome will acquire a particular microbial taxon is given by the relative abundance of that taxon within the community (see below for details on how environmental microbial taxon abundances were calculated).
The probability, x, that a particular "slot" in a new individual host's microbiome was occupied by a microbial taxon sampled from a randomly selected parent was varied across simulations. Two sets of simulations were performed: (1) x and y varied linearly, between 0 and 1, with increments of 0.1 (see S1-S3 Tables for means and standard deviations of diversities); and (2) with values of x, y 2 (0.0, 0.5 [0,1,2. . .10] ) (see S4-S6 Tables for diversities). When x = 0.0, a host's microbiome was sampled directly from the environment, i.e., the probability of a microbial taxon being selected was equal to the relative frequency of the microbial taxon in the environment. When x = 0.5 0 = 1, a host's microbiome was sampled multinomially from the microbiome of its "parent" who was selected with uniform random probability from the previous generation. Similarly, when y = 0.0, there is no contribution from the previous host generation to the environmental microbial resource, whereas when y = 0.5 0 , all microbes are replaced each generation by the pool of microbes resident in the hosts of the previous generation.
Three different models were used to model the relative abundance of taxa in the environment. First, for FE ("fixed environment"), the abundances of the microbial taxa in the environment were fixed to the initial uniform distribution and did not vary over the course of the simulation. Second, under PE ("pooled environment") the abundances of the microbial taxa in the environment were composed of the pooled microbiomes of all the hosts of the previous generation, i.e., by summing over the abundances of the respective microbial taxa in hosts' microbiomes, and renormalizing to relative abundances. Third, under ME y ("mixed environment"), the abundances of the microbial taxa in the environment were calculated by combining the fixed environment and pooled environment with y% derived from the pooled environmental component.
A particular simulation regime consisted of a distinct combination of pooled/fixed environmental ratios and environmental factors. A total of ten replicates were run under each simulation regime. The plots of γ-diversity were inspected, and simulations were halted when these stabilized. The number of generations for each simulation varied between 10 4 and 10 6 generations.
After each generation was simulated, the diversities of microbiomes are measured by scaled Shannon-Wiener index (α-diversity and γ-diversity) or Bray-Curtis dissimilarity index (β-diversity). The scaled Shannon-Wiener index is calculated as À where R represents the total number of taxa and p i represents the relative abundance of ith taxon within the community. The calculation of the Bray-Curtis index is given by the formula: , where n represents the total number of hosts in the population, m represents the total number of microbial taxa within the host population, and p ik and p jk represents the relative abundance of k th taxon within the community of host i and host j.
Empirical data from the human microbiome were obtained from the website of the NIH Human Microbiome Project (http://www.hmpdacc.org/HMSMCP/). The original community profiling data is a table of relative abundances for each of 690 samples and 718 taxa of bacteria and archaea (from 2 kingdoms to 397 species). As is described on the HMP website, all the samples were collected from 16 body sites from 103 healthy humans and processed with Whole Genome Shotgun sequencing. Specific information of each sample is available on the website (http://www.hmpdacc.org/HMIWGS/all/). We selected data of microbial communities associated with anterior nares, buccal mucosa, supragingival plaque, stool, tongue dorsum and vaginal posterior fornix because of their large numbers of samples, and removed the replicated samples from the same human subject for each body site (see S1 Dataset). The HMP data provides relative abundances for different taxonomic levels. We calculated diversities for all taxonomic levels ranging from species to kingdom, and found that the values of diversities for genus, family and order were similar. Consequently, we chose to use genus-level diversities. The relative abundances of genera for each site and all samples is given in S1 Dataset.
Fitting simulated results into log-normal and Dirichlet-multinomial distributions was performed in R with methods "fitdistr" of package MASS and "dirmult" of package dirmult.

(DOCX)
S2 Table. γ-diversity under different combinations of acquisition and environment models with linear scales for MA(X) and ME(Y).

(DOCX)
S3 Table. β-diversity under different combinations of acquisition and environment models with linear scales for MA(X) and ME(Y). (DOCX) S4 Table. α-diversity under different combinations of acquisition and environment models with power scales for MA(X) and ME(Y).
S5 Table. γ-diversity under different combinations of acquisition and environment models with power scales for MA(X) and ME(Y).

(DOCX)
S6 Table. β-diversity under different combinations of acquisition and environment models with power scales for MA(X) and ME(Y).


extfieldsofstudy: Mathematics Medicine

text:
Electrical Brain Responses to an Auditory Illusion and the Impact of Musical Expertise

The presentation of two sinusoidal tones, one to each ear, with a slight frequency mismatch yields an auditory illusion of a beating frequency equal to the frequency difference between the two tones; this is known as binaural beat (BB). The effect of brief BB stimulation on scalp EEG is not conclusively demonstrated. Further, no studies have examined the impact of musical training associated with BB stimulation, yet musicians' brains are often associated with enhanced auditory processing. In this study, we analysed EEG brain responses from two groups, musicians and non-musicians, when stimulated by short presentation (1 min) of binaural beats with beat frequency varying from 1 Hz to 48 Hz. We focused our analysis on alpha and gamma band EEG signals, and they were analysed in terms of spectral power, and functional connectivity as measured by two phase synchrony based measures, phase locking value and phase lag index. Finally, these measures were used to characterize the degree of centrality, segregation and integration of the functional brain network. We found that beat frequencies belonging to alpha band produced the most significant steady-state responses across groups. Further, processing of low frequency (delta, theta, alpha) binaural beats had significant impact on cortical network patterns in the alpha band oscillations. Altogether these results provide a neurophysiological account of cortical responses to BB stimulation at varying frequencies, and demonstrate a modulation of cortico-cortical connectivity in musicians' brains, and further suggest a kind of neuronal entrainment of a linear and nonlinear relationship to the beating frequencies.

Introduction
When two sinusoidal tones with a slight frequency mismatch (i.e. 200 and 210 Hz) are separately presented to each ear, the listener perceives a beating frequency equal to the frequency mismatch between the two tones (i.e. 10 Hz); this is termed binaural beating [1,2]. The two tones combined wax and wane as the two frequencies come in and out of phase with one another, and this phase interference produces an amplitude-modulated standing wave, the binaural beat (BB). This beat is not a physical property of the presented sound but a subjectively perceived auditory illusion, which has its origin in the brainstem's superior olivary nucleus [3,4]. BBs are perceived for low frequency mismatch (< 100 Hz) and lower carrier frequencies (<1 kHz) [1]; phase locking in most mammalian brains becomes less precise at carrier frequencies above 1 kHz and disappears completely at frequencies above 5 kHz [5,6].
BB stimulation is assumed to drive and evoke neuronal responses observable on scalp EEG as steady-state responses (SSRs) over a temporal window [7][8][9][10]. For example, BB at 40 Hz elicited a SSR over fronto-central brain regions [8]. Karino et al. [7] reported SSRs evoked by BBs of 4 Hz and 6.66 Hz that were localized in the temporal, parietal and frontal regions. The authors further suggested the possible cognitive encoding of BBs at a cortical level as the phases of the SSRs showed large variability. Pratt et al. [11] studied event-related-potentials (ERP) of two BBs at 3 Hz and 6 Hz, with two base frequencies, 250 Hz and 2000 Hz, and reported distinct event-related-potential (ERP) components such as P50, N100, and P200 following tone-onsets, with larger brain responses to the lower base frequency and to the lower BB frequency. Draganova et al. [10] showed a 40 Hz auditory steady-state response (ASSR) in the auditory cortices during 40 Hz BB stimulation. All these studies suggest that brain activity corresponding to this illusory auditory beats can be identified from, and systematically affect, EEG scalp recordings. However, we have noticed three big gaps as follows in the neuroscientific literature on BB.
First, most studies have looked only at a limited number of frequency mismatch. The largescale neural oscillatory responses to BB frequencies varying systematically over a relative wide range have not been characterized. We believe this is essential as using a range of BBs would allow us to study the effect of cross-frequency coupling in large scale brain responses. The SSRs primarily represent the interaction between neuronal assemblies and external stimulation with identical frequency, but the cross-frequency coupling suggests a more generalized interaction between neuronal oscillations and stimulating frequency. Indeed cross-frequency coupling at the neuronal level has been demonstrated to provide a general mechanism of functional integration at multiple spatiotemporal scales [12].
Second, no study, to our knowledge, has investigated the cortical network patterns against BB stimulation, yet an increasing body of evidence suggests that functional co-operation among near and distant brain areas underlie almost any sensory and cognitive processing [13,14]. Understanding information processing at the cortical network level provides novel insight into the possible propagation of the auditory stimulation and the induction of the illusory beat.
Finally, little is known about the individual differences in the brain responses to BB stimulation. A recent attempt in this direction was made by Goodin and colleagues [15] who investigated EEG spectral power in thirty three participants exposed to BB in either theta (7 Hz) or beta (16 Hz) frequency for a period of 2 minute and measured the personality traits by standard Big Five Factor model; no significant relationships were found between BB-related EEG responses and Big Five personality traits. Instead of looking for personality dependent differences, a more pragmatic and useful approach would be to investigate training related differences. In this context, musicians are an important group to study as the musician's brain is a robust model of neuroplasticity [16,17], and long-term musical training leads to 'auditory fitness', i.e. an enhanced processing of auditory information [18].
The current study was designed to address these issues by recording electrical brain responses as measured by high-density EEG signals obtained from healthy human adults, both musicians and non-musicians, when stimulated by a range of BBs from low frequency delta band to high frequency gamma band. We analysed the EEG signals in terms of neuronal oscillatory activity as measured by spectral power, and phase synchronization as measured by mean phase coherence [19] and phase lag index [20]. Further, we also characterized the underlying network pattern by measures based on graph theory approach [21].
In this study, we strategically focused our analysis on the neuronal oscillations in the alpha and gamma frequency bands of EEG signals for the following reasons. Firstly, large scale brain oscillations in the alpha band (8)(9)(10)(11)(12) are the most studied and understood brain rhythm [22] and have been found ubiquitously in most cortical and subcortical areas [23]. Further, when humans are visually stimulated by flickering stimuli, steady state oscillations in the alpha band and its higher harmonics exhibited most resonant effect [24]. So it would be important to investigate whether alpha oscillations show similar resonance effect in auditory domain. This is in clear alignment with the growing interest on the importance of studying alpha rhythms in audition (see for a review, [25]). Furthermore, alpha oscillations are crucially associated with the perception of various types of illusions, i.e. examples of perceptual accounts that cannot be explained by physical properties of the stimuli, from visual (e.g., wagon wheel illusion [26]), auditory (e.g., tinnitus [27]), to multisensory (e.g., double flash illusion [28]). As alpha oscillations reflect the degree of excitability of sensory cortices, we expected that it would be associated with the perception of the binaural beat, an auditory illusion. Second, there has been increasing evidence that the degree of neuronal synchronization in the gamma band (30 Hz and above) is significantly larger in trained musicians as compared to non-musicians during music related task processing [29][30][31]. Further, enhanced neuronal synchronization in the gamma band has also been observed in other expert groups like long-term meditators [32], proficient bilinguals [33], professional artists [34]. Considering the widespread role of gamma band in providing a general framework for cortical computation [35] and its possible modulation with neuroplasticity, we expected that the gamma band synchrony would be differentially involved between the musicians and non-musicians during the perception of illusory binaural beat. Therefore, we had two specific predictions: (i) across both groups, EEG alpha band power would be largest when stimulated by BBs belonging to alpha band, and (ii) as compared to non-musicians, musicians would produce a larger network response against binaural beat in the gamma band.

Participants
Sixteen musicians (6 males, age 25.5 ± 3.18 years) and sixteen non-musicians (7 males, age 26.1 ± 3.82 years) participated in the study. All but one non-musician were right-handed. No participants reported any auditory deficits or neurological disorders. Our musicians reported an average of 18.4 (± 4.3) years of active engagement with a musical instrument and 14 (± 3.40) years of formal training on a solo instrument. They practiced their principal instrument an average of 21.59 (± 11.42) hours per week with 81.3% reported being professional musicians. All participants were kept blind to the real purposes of the experiment. They gave their written informed consent before the start of the experiment. The experimental protocol was approved by the Ethics Committee of the Department of Psychology at Goldsmiths, and the experiment was conducted in accordance with the Declaration of Helsinki.

Stimuli
Auditory stimuli were presented in 34 blocks of 2 min and 20 s each. Each block consisted of a silent period of 20 s followed by two auditory conditions of 1 min each: the non-binaural beat (NB) and the binaural beat (BB). For both auditory conditions (NB and BB), the base frequency (f b ) was 200 Hz. This selection was made after recent studies on BB stimulation using f b between 100-300 Hz [7,11,36]. For the NB condition, both ears received auditory pure tones of 200 Hz across all blocks. For the BB condition, the left ear received an auditory pure tone of 200 Hz and the right ear an auditory pure tone of 200+f bb Hz (f bb = beat frequency) with f bb varying with a step size of 1 Hz in the range of 1 Hz to 20 Hz and with a step size of 2 Hz from 20 Hz to 48 Hz. Based on the values of the beat frequency (f bb ) lying within the standard EEG frequency bands [37], we formed five groups of BB stimuli as follows: delta-BB (1 Hz f bb 4 Hz), theta-BB (5 Hz f bb 8 Hz), alpha-BB (9 Hz f bb 12 Hz), beta-BB (13 Hz f bb 30 Hz), and gamma-BB (32 Hz f bb 48 Hz). Following this grouping, we had four blocks each for delta-BB, theta-BB and alpha-BB, eight blocks for beta-BB and finally fourteen blocks for gamma-BB. Auditory stimuli were presented through a Philips in-ear headphone with rubber caps that reduced external noise, ensured a precise fit, and decreased effective surface for bone conduction [38,39]. The volume of the auditory stimuli was self-adjusted a priori by each participant and was kept at constant level throughout the experiment. During the entire period of auditory stimulation, participants watched a silent movie with subtitles on and were also asked to ignore the background auditory stimuli. These altogether ensure that any effect elicited by BBs would be largely implicit, i.e. not requiring overt attention. At the end of the experiment, participants rated on a 5-point scale their responses concerning the pleasantness of the background sounds, the interestingness of the movie, and their overall alertness, and finally completed a self-report questionnaire regarding personal details including musical expertise.
Recording and pre-processing EEG signals were recorded from sixty-four active electrodes that were placed according to the extended 10-20 system, and amplified by a Biosemi ActiveTwo amplifier. Four additional electrodes were placed around the eyes to record vertical and horizontal eye movements. The sampling frequency was 512 Hz. Biosemi system has two electrodes-active CMS (common mode sense) electrode and passive DRL (driven right leg) electrode-that together form a feedback loop representing the online reference (see http://www.biosemi.com/faq/cms&drl.htm for details on the Biosemi referencing and grounding procedures). The EEG data were later referenced offline to the average of the left and right earlobe electrodes. The data were further down-sampled offline to 256 Hz to make the file size manageable. A 0.5 Hz high-pass filter was applied to remove drifts and low-frequency artefacts. Trials containing large artefacts (< 2%) were eliminated by visual inspection from subsequent analysis. Eye-blink artefacts were corrected using Independent Component Analysis as implemented in the EEGLAB package [40]. For each auditory condition (NB and BB) of 1 min long, epochs were extracted by excluding the first and last 500 ms to eliminate transient brain responses due to the sudden onset or offset of stimulus, resulting in epochs of 59 s each.

Data analysis
EEG power analysis. The spectral content of the EEG signals was estimated by the method of multitapers [41] as implemented in the Matlab function 'pmtm'. Multitapers are sets of functions that reduce the leakage between neighbouring frequencies, thereby making it suitable for estimating the oscillatory content of neuronal signals [42]. We used a time-bandwidth parameter of 4, a 1024 point FFT (4 s window) with a frequency resolution of 0.25 Hz. We divided the broadband EEG power spectrum into five standard EEG frequency bands [37]: delta-EEG (1-4 Hz), theta-EEG (5-8 Hz), alpha-EEG (9-12 Hz), beta-EEG (13-30 Hz), and gamma-EEG (32)(33)(34)(35)(36)(37)(38)(39)(40)(41)(42)(43)(44)(45)(46)(47)(48); in this article, we focused our primary analysis on alpha-EEG and gamma-EEG frequency bands. The BB-related spectral power was normalized by the NB-related spectral power, and this normalization was done separately for alpha-and gamma-EEG bands, electrode, and at participant level. The normalized spectral power was expressed in dB. The auditory SSRs was estimated by averaging the EEG spectral power over alpha (gamma) frequency band when participants were stimulated by binaural beats belonging to the alpha-BB (gamma-BB). All statistical analyses were initially performed at the global level after averaging across electrodes. The statistical significance level was set at P < 0.05.
Phase synchrony analysis. For the assessment of phase synchronization (PS), it is first necessary to estimate the phases of the EEG signal, x i (t), at each electrode-i, a procedure that consists of two steps. First, the raw EEG data are band-pass filtered (x i,α (t)) in the frequency band of interest (say, the alpha band, α) using a finite impulsive response (FIR), zero-phase distortion filter. Then, this real valued filtered data x i,α (t) is converted into a complex-valued one, where j = ffiffiffiffiffiffi ffi À1 p is the imaginary unit and x H i;a ðtÞ is the Hilbert Transform of the filtered data: p.v. stands for Cauchy principal value. The phase ofx i;a ðtÞ, θ i α (t) is then defined as: The cyclic relative phase (i.e., restricted to the interval [0,2π)) between EEG electrodes i and k (i, k = 1,. . ., 64) is finally obtained as: Assessing the degree of PS between two electrode regions comes down to estimate whether the distribution of Eq (4), which is a circular variable, is different to what would expected for two phase-independent signals. This estimation can be performed in different ways [43,44]. An index, widely used in M/EEG analysis, is the mean phase coherence [19], also termed as Phase Locking Value (PLV); it is a measure of how homogeneously the relative phase spreads over the unit circle: where <> stands for average value and | | for absolute value. This index ranges between 0 (no phase synchrony) and 1 (perfectly synchronized in phase). Though Eq (5) is a powerful, widely used PS indicator, it does not discriminate between zero phase lag and non-zero constant phase lag, yet two EEG signals can be phase-synchronized in both ways. Further, volume conduction effects (i.e. a single neural source affecting two or more EEG electrodes) primarily contribute to zero phase lag synchrony. Therefore, another index, phase lag index, PLI [20], has been developed to deal with this issue, by taking into account that true interaction between neural sources, as opposed to volume conduction effects, occurs with some delay [45], which in turns gives rise to a distribution of Eq (4) asymmetric around 0 (or π). Thus, PLI is defined as: Although, from the point of view of functional connectivity analysis, it may seem enough, then, to look at Eq (6), zero lag synchronization can also be achieved if two neural sources are indirectly connected through a third one, which acts as a dynamical relay [46]. This is actually not a volume conduction effect, but rather an indirect connection between, e.g., two cortical sources through the thalamic relay [47], which would be overlooked should we only focus on PLI. Thus, both indices are complementary, with PLI estimating only direct true connections, whereas PLV is also sensitive, if any, to zero lag, indirect connections [48].
Both indices were calculated using the recently released HERMES toolbox for Matlab [49].

Characterization of the functional brain networks
Once the degree of PS between any two electrodes was assessed by means of PLV or PLI, we had the corresponding interdependence matrix: where a stands for either PLV or PLI, 0 a ij = a ji 1 (i,j = 1,..,n, i 6 ¼ j), n = 64, and we take a ii = 0. As it is well-known, in the context of multivariate EEG analysis this matrix can be regarded as the weighted adjacency matrix of a complex network, where the electrodes are considered as nodes and the PS indices measure the strength of the links between them [21]. From this adjacency matrix, it is possible to estimate different measures that provide information about network's structure and function.

Thresholding
The first step in this direction consists in determining which of the values of the PS indices should be considered significant. There are different approaches for this largely unsolved question [50,51]. A straightforward approach consists of estimating the significance of each individual link a ij at a given level of statistical significance by using, e.g., bivariate surrogate data [52,53]. The problem with this approach, however, is that one may end up with adjacency matrixes with different proportions (densities) of significant links across groups/ conditions, which may bias the values of the network measures [54], thereby giving rise to spurious statistical differences between them. This would be also the problem if we were to take a fixed threshold value T, so that a ij was set to zero if a ij < T. Instead, we used here the so-called fixed density approach [50], in which values of a ij are rank ordered for each adjacency matrix, and the (1-k) Ã N LINKS lowest values were set to zero (0<k<1 is the fixed density, N LINKS = 64 Ã (64-1)/2 = 2048 is the total number of possible links). In line with some recent studies [55,56], we calculated all network measures (see below) for different values of k to check that the robustness of the results against this parameter.

Network measures
We described the structure of the weighted undirected functional brain networks resulting from the above procedure by calculating three commonly used measures (the strength, the clustering and the efficiency), which characterizes network centrality, segregation and integration, respectively [21].
The average strength (S) of the network is defined as: where S i , the strength of each individual node, is the sum of weights of links connected to it, and provides a simple centrality measure that estimates the importance of the node within the network. The clustering coefficient (C) quantifies the tendency of network elements to form local clusters. The presence of network clusters indicates segregated functional dependencies in the brain. In its weighted version [57], C is defined as: where C i is the clustering coefficient of node i (C i = 0 for k i < 2), k i is the degree of node i and t i is the sum of triangle intensities around node i, defined by the geometric mean of triangles around it: Finally, we measured integration in the network by means of the local efficiency [58], which can be understood as a measure of how well each sub-graph exchanges information when the index node is eliminated [59]. It is inversely related to the shortest average path lengths between nodes in the network that contains only neighbours of node i.
All the network measures were calculated using the Brain Connectivity Toolbox for Matlab [21].

Statistical tests for the network measures
The comparison of brain networks can be carried out at two different levels: global and local [60]. Indeed, the three network measures described above are all calculated locally (at each node) and then averaged to get a global estimation of the corresponding feature. Correspondingly, we perform the statistical tests at these two levels.
At the global level, we used a mixed 2x2x5 factorial ANOVA, with Group (musicians vs. non-musicians) as between-subjects factor and Beat (non-binaural vs binaural beat) and Range (delta-BB, theta-BB, alpha-BB, beta-BB and gamma-BB) as within-subjects factors. When any of the factors (or their interaction) was significant at the P<0.05 level, differences were further analysed using appropriate post-hoc comparisons. The test was applied to the global measures for both types of BBs, alpha-BB and gamma-BB. Local comparisons at the sensor level were carried out to elucidate the topography of the global differences, in those cases were the corresponding omnibus test was significant. In this case, a paired (for Beat and Range) or an unpaired (Group) t-test was applied to the three local measures at each electrode, and then corrected for multiple comparison using a Type-I False Discovery Rate (FDR) at the q<0.05 level [61].

Behavioural analysis
No statistically significant differences were found between the two groups (musicians vs. nonmusicians) in terms of their reported pleasantness of the background auditory stimuli (Mann-Whitney U-test: z = 1.25, P > .05, two-tailed), interestingness of the movie (z = 1.75, P > .05), and alertness over the entire duration of the experiment (z = .14, P > 0.5).

Power analysis
First, we calculated normalized SSRs for both alpha-BB and gamma-BB averaged across all electrodes and subsequently applied a 2x2 mixed factorial ANOVA with Group (musicians vs non-musicians) as a between-subjects factor and BeatFrequency (alpha vs gamma) as a withinsubjects factor. We found a significant effect of BeatFrequency: SSR for alpha-BB was larger than SSR for gamma-BB (F(1,32) = 5.03, P = .03). We did not find any effect of Group or Group Ã Beat interaction (P > .6). Next, we investigated normalized SSR individually, as any value systematically larger than zero would suggest a significant frequency following response for that frequency band specific BB stimulation. Two separate one-sample t-tests were conducted and only the alpha-BB revealed a significant effect in its SSR (alpha: t(31) = -3.15, P = 0.004; gamma: t(31) = -1.54, P = .13).
Next we analysed the normalized alpha-EEG power at the global level (i.e. averaged across all electrodes) separately by a 2x5 mixed ANOVA with factors, Group (two levels) and BB (five levels: delta, theta, alpha, beta and gamma). We observed a significant effect of BB (F(4,120) = 4.39, P = .002); the largest increase of normalized alpha-EEG power was observed for alpha-BB followed delta-BB (Fig 1(A)) and only these two showed a significant non-zero response (alpha-BB: t(31) = 3.15, P = .004, delta-BB: t(31) = 2.67, P = .01). The scalp maps of these two effects are shown in Fig 1(B).
No robust differences were observed in the normalized alpha-EEG power between the two groups except for gamma-BB for which the normalized alpha-EEG power in the musicians was higher than non-musicians (t(30) = 2.13, P = .04).
For normalized gamma-EEG power at global level, we did not observe a significant effect of BB (P > .8) but we did find a marginal effect of Group (F(1,30) = 2.4, P = .048) as musicians showed overall higher normalized gamma-EEG power than non-musicians.
Interestingly, when we extended our power analysis to other EEG frequency bands (delta-EEG, theta-EEG and beta-EEG, see Methods), we did not observe significant effect (at the level of Bonferroni corrected P) on the normalized SSRs ( Figure A in S1 File) nor we found any significant effect of cross-frequency responses at other EEG frequency bands ( Figure B in S1 File). The robust results were only found at the alpha-EEG band.

Network analysis
At the global level, differences in network measures were only found for the alpha-EEG band for both phase locking value (PLV) and phase lag index (PLI), as shown in Table 1 for PLV,  Table 2 for PLI, and Figs 2 and 3. In the case of the PLV, the pattern of statistical differences is complex. Thus, the main effect of Beat was significant for all the densities for the network strength (S) and the network efficient (E) indices, but was not present for the network clustering (C). In turn, the main effect of Range was also significant across densities for S, but depended on k for both C and E. The interaction between factors was, however, robust across all the densities analysed for both indices, although the profile was different for the PLV and the PLI. For the PLV (Fig 2), there was a significant Beat-Range interaction (cfr. S: F(4,120) = 5.73, P < 0.001; C: F(4,120) = 5.98, P < 0.001; E: F(4,120) = 4.53; P < 0.01 for k = 0.8), so that BB vs NB differences were only significant in the delta-BB range, with BB > NB for the three network measures. For the PLI (Fig 3), there were no main effects for any index. There was, however, an interaction among the three factors (Beat-Range-Group), with BB < NB for musicians in the theta-and alpha-BB ranges (cfr, S: F(4,120) = 3.9, P < 0.01 for k = 0.8).
The topologies of the differences at the electrode level are shown in Fig 4 for PLV and Fig 5 for PLI. In both cases, they correspond to the density k = 0.8, but are robust across the other values of k analysed. For the PLV, the topology is strikingly similar for the three measures, with differences spreading over many electrodes, the greatest ones showing clear right hemisphere lateralization and antero-frontal and centro-parietal localization, with both C (not shown) and E showing also a left parieto-occipital patch. For the PLI, however, the pattern is different for S (Fig 5) than for the other two measures (Fig 6 for C). Thus, differences in S were mainly concentrated along the midline electrodes (both theta-and alpha-BB range), with an additional right temporo-occipital patch in the alpha band; for C and E, although differences showed a certain trend for right hemisphere lateralization (more posterior for the theta than for the alpha-BB range), they were widespread over many electrodes.
Strength (S), clustering (C) and global efficiency (E) for the PLV index as a function of the density of network links, k. The results are for alpha-EEG band.
Asterisks stand for differences for the main effects (B)eat (Nonbinaural vs. Binaural) and (R)ange (low frequency delta-, δ-, to high frequency, γ-BB), as well as the interaction between them (B-R) *: p < 0.05 **: p < 0.01 ***: p < 0.001. In the B-R column, Greek letters indicate the BB range for which BB-NB differences were significant.
doi:10.1371/journal.pone.0129486.t001  Greek letters indicate the BB range for which BB-NB differences were significant. Note that the differences and the bands were robust across all the density links considered.

Discussion
The current study investigated the relationships between BB stimulation and large-scale neuronal oscillations and synchrony, and the potential influence of musical training on these relationships. We investigated the alpha and gamma band oscillations in the EEG signals against a range of BB stimulations; both local (in terms of spectral power), long-distance (in terms of phase synchronization) neuronal synchronization were investigated in addition to characterizing the underlying brain network in these two bands of neuronal oscillations. Finally, we also investigated the impact of musical expertise on the brain responses to BB stimulations.

Power analysis
We have identified two types of oscillatory responses: (i) narrow band or frequency following responses suggesting a kind of neuronal entrainment or a linear relationship to the beat frequency, and (ii) broadband or cross frequency responses suggesting a kind of nonlinear relationship.
As predicted, we found significant entrainment in the EEG-alpha band across participants; the alpha band EEG power was significantly enhanced during alpha-BB stimulation. Early  Table 2). Asterisks stand for differences between BB and NB, ***: p<0.001. research literature indicates that presentation of a visual rhythmic stimulus (i.e. a flickering light) at alpha band frequencies lead to a robust enhancement of EEG power in alpha and harmonic or sub-harmonic frequencies of alpha [62,63], known as photic driving. However, the evidence for the similar driving response caused by auditory rhythmic stimuli has been mixed [15,[64][65][66][67]. The reported absence of BB-related entrainment could be due to several methodological issues, such as using a higher f b of 400-500 Hz instead of 100-200 Hz, single BB stimulation instead of a range of stimulation or sparse spatial sampling of brain regions. By systematically varying the BB frequency and using high density EEG recording followed by carefully constructed statistical analysis, we here demonstrated the role of alpha band in auditory driving as well.
Alpha, being the dominant neuronal oscillations in the human brain, is classically considered as the 'idling rhythm' [68,69]. However, a large body of recent evidence suggest that alpha band oscillations play an active and important role in both task-specific and ongoing information processing (see [70] for a review). Our alpha effect was largest over posterior parieto-occipital electrode regions, roughly overlying visual cortices. Previous research suggests that parieto-occipital alpha power increases with anticipatory attentional allocation mechanisms [71] and distractor suppression [72]. Altogether, our results are also aligned with behavioural studies suggesting that alpha band rhythmic stimulation could enhance attentional [73] and memory [74,75] performance, and reduce the sensation of unpleasant stimuli [76].
In terms of cross frequency responses, alpha power was significantly enhanced by delta-BB stimulations, with the effect being predominant over parietal electrode regions extending to right frontal regions. By using a variety of auditory stimuli with and without rhythmic component, it was shown that listening to rhythmic stimuli, as compared to non-rhythmic ones, was associated with higher EEG alpha power, and further, the auditory envelope was synchronized with the alpha oscillations [65]; this suggests some kind of frequency selectivity of brain oscillations against auditory or visual stimulation [77], something that has been put forward as one of the possible models of an alpha generator [78]. The spatial locations of this effect, the frontoparietal regions, are aligned with a previous study in which auditory SSRs were evoked by BBs of 4 Hz and 6.66 Hz [7]. A more recent fMRI study has found that acoustic exposure to a more isometric rhythmic sequence was associated with activations in the parieto-frontal brain regions, including dorsal premotor cortex and supplementary motor cortex [79]. An isometric sensation can also be perceived while listening to BBs due to the interaural phase differences (IPD). Bengtsson et al. [79] suggested that activation in these brain areas might reflect the temporal sensory prediction.
Further, the BBs in lower frequencies (e.g. in delta, theta and alpha frequency bands) often elicit a subjective sensation of a moving 'image' that swaps from one side to another side (laterally) within the head [80]. This internal sensation of a moving object from place to place in the head may partly explain the engagement of the fronto-parietal brain network. In fact, an ERP study showed that the sensation of a moving sound as caused by the interaural time differences (ITD) was associated with right inferior parietal and bilateral inferior frontal activities [81]. Nevertheless, further examination is needed to establish the link between consciously perceived sensation of motion induced by BB stimulation and such brain activity patterns. Similar evidences of cross-frequency function of BBs were also found by Kasprzak [36], who examined the effect of BB in 10 Hz. Results revealed higher activation of the EEG 10 Hz strength signal. Interestingly, averaging amplitudes of spectral density into frequency bands lead to a significant decrease for alpha and beta bands and increase for theta frequency bands. He suggested that this phenomenon of alpha suppression in parallel with a 10 Hz follow-up effect might be a natural reflex of the nervous system to be tuned to a strong external stimulus. Cross-frequency results were also reported by Atwater [82] who found a reduction of occipital alpha when stimulated by a complex pattern of BBs mainly in delta band. These results may also give an explanation to our findings underling the non-linear dynamic function of the brain under various BB stimulations. Moreover, encoding of BBs may also be a process of a more complicated mechanism involving both a cognitive process, as Karino et al. [7] suggested, and an encoding of the IPD in the ascending auditory pathway [3,4,10].
Comparisons between musicians and non-musicians revealed higher alpha power over the frontal and parietal regions for gamma-BB stimulations (figure not shown). Previous fMRI studies revealed that frontal and parietal lobes show higher activation bilaterally in musicians as compared to non-musicians during passive listening to piano melodies [83]. Similarly, Meister et al. [84] reported bilateral activations of the parietal-frontal network during music performance and musical imagery. An additional explanation why alpha power increases in musicians across gamma-BB stimulation may be that stimuli higher than 30 Hz (gamma-BB stimulation) falls into the audible range (20 Hz-20 kHz) of the humans' auditory system and is effortlessly recognizable as a bitonal instead of binaural sensation. The frequency intervals, which resulted from the mismatched frequency between both stimuli, are easier to be detected and distinguished as two different tones, especially for well-trained musicians [85]. sensation of a ± minor 3 rd musical interval of the well-tempered tuning system (minor 3 rd ; ratio: f6/5 (f-reference = 200 Hz) =~240 Hz, scale pitch names: near to G 3 and B [ 3 ). Therefore this bitonal sensation could elicit higher frequency selectivity for musicians who are exposed much longer to musical sounds, and therefore could lead to higher alpha power across gamma-BB stimulation. Musicians also showed marginally higher overall gamma power across BB stimulations. This is aligned with previous findings that musicians present larger gamma band oscillations [86] and synchronizations [29] while listening to musical stimuli.

Network analysis
To our knowledge, this is the first study on characterizing the cortical network pattern when stimulated by binaural beats. Our second prediction of finding a greater network in the gamma band oscillations for musicians as compared to non-musicians was not supported. Instead we have found that the differences between the groups and also the effects of binaural beat were found in the alpha band EEG only.
Before we interpret our findings of network analysis, we would like to include some remarks on analysis methodology. It is not exaggerated to state that the network analysis based on graph theoretic methods has been at the forefront of the current neuroimaging research and it has become customary in multivariate EEG studies [87][88][89]. Despite the extensive literature in the subject, however, there are still many important open questions in the field [50]. Among them, we can mention the choice between weighted and unweighted networks, choosing a threshold to consider an edge significant, the effect of network size, which network index best characterizes between group (or tasks) differences, how to deal with the multiple comparison problem and even which functional connectivity index is the right one for a given problem [51,54,90]. Here, we had three popular network indices characterizing different network properties (strength, clustering, efficiency [21]). Moreover, in line with recent EEG studies [55,56], we have used weighted networks, which provide additional information on edge weights as compared to Boolean networks, and have checked that the statistical differences are consistent for a wide range of degree densities (see Tables 1 and 2), both at the global and at the local level, at which we corrected for multiple comparisons by means of the type I FDR algorithm. One important feature of our study, however, is that we analysed the functional networks obtained from two different PS indices, whereas the usual approach entails the use of either the PLV or the PLI. At first sight, it may seem a little odd to use them jointly, since PLI was developed as a kind of refined version of the PLV, robust against to volume conduction due to its insensitiveness to zero lag correlations [20,44]. Yet we think (and the results seem to agree with us) that there are good reasons to include both the PLI and the PLV in the analysis (see also [49] for a recent example), because, as commented in the Methods section, both indices provide complementary rather than redundant information. Indeed, whereas PLI is in principle the right choice for M/EEG graph connectivity analysis as it focuses on lagged (direct) correlations [91,92], it may miss true, zero lag cortico-cortical connectivity through a common relay (i.e., the thalamus [93][94][95]), which could be picked up by the PLV. Thus, by analysing both functional brain networks (the "true" one derived from PLI and the one derived from PLV, which may be affected by volume conduction effects, but includes all type of connections), we can be sure to gather most of the information available from our data in the form of functional connectivity patterns. Indeed, recent results comparing the application of PLV to PS indices robust to volume conduction, clearly suggest [96,97] that PLV may provide additional information not present in other indices, whereas the use of advanced pre-processing methods such as the Laplacian to work with current source densities instead of voltages does not solve the volume conduction problem either [98].
In concrete terms, our results showed that both indices agreed in showing that statistical differences are only present in the alpha frequency band and for the lower range of binaural beats (delta to alpha), with a slight right hemisphere dominance. But the network of direct corticocortical connections (as assessed by the PLI) presented differences only for the group of musicians (Table 2, Figs 3, 5 and 6), with a decrease in overall connectivity, clustering and efficiency during the processing of the binaural beat. Yet this decrease ran in parallel to a (group-unspecific) increase of all network indices for the PLV network (Table 1, Figs 2 and 4), although in this case an alternative explanation might be that during BB there might be an enhanced volume conduction effects due to increased deep source activation. Taken together, these results indicate that processing of low range (delta to alpha) binaural beats significantly increased zero lag (i.e., mediated by deeper relay neural sources) phase synchrony between cortical areas (and possibly also the overall oscillatory activity of these deeper sources) in the alpha band for all subjects regardless of the level of training level in sound processing, but a concomitant reduction in direct cortico-cortical alpha band connectivity is indeed related to the level of expertise in sound processing.

Limitations and future extension of the research
We would like to mention here a few practical remarks. First, sound intensities of this study were not kept constant across participants, instead were self-adjusted by participants at individual level of comfort. Higher intensities may generate possible cross-hearing (bone or fluid conduction) which may affect the presented data [38]. Second, the mismatch frequency was always presented to the right ear, thereby limiting any interpretation on lateralization effect [99]. However, this stands for all of our participants, so any group related effects are unlikely to be affected by this presentation bias. Third, our musicians group was heterogeneous in terms of type of instruments played and years of musical experience, which could explain the marginality of the reported group related effects. Future studies could be focused on different subcategories of musicians such as those with absolute vs. relative pitch [100,101], years of experience [102], type of instrument played [16,103], or even gender [104,105]. Fourth, in addition to musical expertise, the two groups may also differ across other dimensions, like personality. Previous research does suggest that personality traits may partly explain individual differences in the entrainment induced by BB stimulation [106,107], though a recent study did not find any significant relationship between neuronal entrainment and personality variables against a short (2 min) BB stimulation [15]. Finally, our results revealed the neural correlates of BB processing, but due to the lack of behavioural measures related to the reported BB illusion, the specificity and potential implications of the reported effects are less than explicit. On the other hand, these effects could be considered as markers of implicit processing of BB stimuli.

Conclusions
In summary, the present study revealed that binaural beat stimulation can modulate the strength of neuronal oscillations and synchrony obtained noninvasively from the scalp EEG. Various regions of the brain, such as frontal, temporal and parietal lobes, seem to be involved in the procedure of binaural auditory mismatched frequencies stimulations. Alpha brain oscillations appear to be the most prominently entrained to the perceived binaural beat illusion, and further, alpha band network appears to be significantly modulated by low frequency binaural beats.
Supporting Information S1 File. BB-related SSRs and cross-frequency responses. Figure A: Normalized SSRs in each frequency band across all participants, error bars indicate 99% confidence intervals. A significant effect was observed during alpha-BB stimulation. Figure B: Cross-frequency responses against BB stimulations, error bars indicate 99% confidence intervals. A significant enhancement of alpha-EEG power during delta-BB stimulation was found. (DOCX)

Author Contributions
Conceived and designed the experiments: JB CII. Performed the experiments: CII. Analyzed the data: EP JB JPL CII. Contributed reagents/materials/analysis tools: JPL EP JB. Wrote the paper: JB CII EP. Overall research supervision: JB.