Mobile microbiome networks across spatiotemporal gradients

Background: The microbiome, a community of microbes that co-reside in biotic or abiotic environments, underlies biogeochemical cycling, plant and animal development, and human health. Increasing evidence shows that much of the role the microbiome plays is executed through complex interactions among microbes. Thus, network reconstruction has been increasingly used as a tool to disentangle internal workings within microbial communities. Results: We developed a general framework for recovering microbial interaction networks from any design of microbial experiment. This framework represents a quasi-dynamic game model (qdGM) derived from the seamless integration of evolutionary game theory and allometric scaling laws. The qdGM can not only characterize how individual microbes act singly, but also reveal how different microbes interact with each other to govern microbial community assembly. Beyond existing approaches that can only identify a single overall microbial network from a number of samples, our framework can track and visualize how interaction networks vary from sample to sample and covert sample-specific (personalized) networks into context-specific networks. More importantly, this framework can reconstruct such mobile microbial networks from steady-state data, facilitating the widespread use of network tools to understand the impact of the microbiome on natural processes. Conclusions: As proof of concept, we used the new framework to analyze human gut microbiota data and interspecific animal-associated microbiota data. Mobile networks reconstructed from each dataset can characterize previously unknown mechanisms that drive the change of microbial interaction architecture and organization along spatiotemporal gradients. This framework provides a tool to generate process-specific microbiome networks that can be readily translated into various biotechnological applications and evolutionary studies. or not informative for biological modelling. We develop a statistical model that can recover sample-specific informative networks from steady-state data. We integrate allometric scaling theory and evolutionary game theory to quantify sample-specific changes of microbial abundance and interactions by deriving a system of quasi-dynamic ordinary differential equations (qdODEs). We integrate functional clustering and variable selection to identify distinct network communities from high- or ultrahigh-dimensional microbiome data and the most significant links a given microbe receive to form sparse networks. We incorporate a maximum likelihood approach for estimating qdODE parameters that constitute the stability of microbial networks. We formulate a statistical platform for comparing and testing how microbial networks evolve over a phylogenetic clade of hosts, how they change across spatiotemporal gradients, and how they establish causal interrelationships with biogeochemical processes. We derive the mathematical properties of qdODE solving and perform extensive computer simulation studies to investigate the statistical behavior of microbial networks. The networks inferred from model can reveal the full information of microbial interactions, interrogate their ecological properties, and serve as an approach for visualizing how microbial communities assemble, stabilize, and evolve in a range of biological and biotechnological contexts.

microbes vary independently and how they interact with each other over natural processes.
Variable selection equips qdGM with a computational capacity to reconstruct highdimensional but sparse microbial networks. Taken together, not relying on temporal or perturbed data, our framework can reconstruct mobile microbial networks that topologically vary among hosts, environments, and times.

The qdGM formulation
Networks are snapshots of biological processes and particularly central to the functionality of microbes in ecological communities across time and space. Thus, the reconstruction of dynamic networks is of paramount importance to better understand and reveal general principles that mediate microbe-environment (biotic and abiotic) relationships. Not relying on temporal data, the qdGM can infer dynamic microbial networks. To formulate this theory, it is useful to consider a series of n samples, each representing an independent microbial community, which are drawn to reflect a spectrum of ecological or health processes. In ecology, a habitat is defined as the sum of all resources and environmental conditions in a community all species can jointly exploit to live and grow [20], whereas a niche defined as the subset of those factors that supports a specific species' survival and growth [21]. Being part of the habitat, different niches may have the same or different types and numbers of environmental dimensions (including biotic and abiotic) [22], which establishes the material, energy, or signal basis for various species to interact with each other through cooperation or competition. In this regard, the pattern of how each niche (part) scales with the habitat (whole) determines the type and strength of interspecific interactions. We argue that the allometric scaling of a niche with the habitat as a part-whole relationship can be explained by resource allocation theory [23] that obeys the power law [24].
Under natural conditions, modeling the niche-habitat relationship via its underlying environmental components is impossible given the confounding effects of multiple environmental factors. Since the number of microbes that inhabit the community virtually reflects its carrying capacity to supply essential resources and conditions that maintain these microbes' survival and growth [25,26], we coin the total amount of abundance of all microbes in a sample as the habitat index (HI) of this sample and the abundance level of a specific microbe in the sample as the niche index (NI) of this microbe within this sample. Let Mji denote the abundance level (or NI) of a microbe j (i = 1, .., m) from sample i (i = 1, .., n) and = ∑ =1 denote the HI of sample i. Thus, modeling how Mji scales with Wi can facilitate the characterization of how microbes diversify and interact with each other across samples. The concepts of HI and NI, analogous to the environmental index defined to quantify the overall quality of site in terms of the accumulative growth of all plants [27,28], can reflect the overall adaptation of all microbes or a specific microbe to the microbial assembly environment.
The microbes, co-inhabiting the same environment, promote or inhibit each other in a complex network. This process acts like a game in which each player may choose to compete or cooperate with its opponents to maximize its payoff. Such strategic choices include a rational judgement based on a player's accrued knowledge of the environment affected by other players [29]. Evolutionary game theory models the dynamic change of strategy frequency in the population as a result of evolutionary pressure [19]. In an evolving population, any strategy used by an individual to maximize its payoff would be constrained by strategies of other individuals that also strive to maximize their own payoffs and, ultimately, this process through natural selection, would optimize the structure and organization of the population, making it reach maximum payoff [19]. We integrate such evolutionary game theory and concepts of NI and HI through Lokta-Volterra prey-predator equations [30] to develop a system of ordinary differential equations (ODEs), expressed as where i denotes the ith sample of context c when this context is considered.
Thus, we will have C such context-specific microbial networks. We can test and compare how microbial networks change under different treatments, with environmental signals, and over times. All this is of fundamental importance for our detailed understanding of microbial impacts on natural processes.
An approach for testing the relationship between microbial networks and contexts is to draw the HI-varying curve for each microbe and decompose it into its endogenous and exogenous curve components. By comparing these two types of curves, we can assess how the existence of microbes is influenced by others and how this influence is determined by contexts.

Ecological and social dissection of microbial networks
In the context l-specific network, ( ) represents the intrinsic capacity with which microbe j grows in isolation, and ′ ( ) reflects the pattern and strength with which microbe j is affected by microbe j under context l. If ′ ( ) is positive or negative, this suggests that microbe j is promoted or inhibited by microbe j, respectively. The zero value of ′ ( ) suggests that microbe j is neutral to microbe j. The size of ′ ( ) describes the strength of promotion or inhibition. By comparing ′ ( ) and ′ ( ), we can determine whether and how these two microbes reciprocally trigger impacts on each other. If both are positive or negative and equal in magnitude, this suggests that microbes j and j have established a relationship defined as symmetric mutualism or symmetric antagonism. If ′ ( ) and ′ ( ) are positive or negative but not equal in magnitude, these two microbes have a relationship defined by asymmetric mutualism or asymmetric antagonism, respectively. If ′ ( ) is positive or negative but ′ ( ) is zero, this suggests that microbe j, respectively, exerts commensalism or amensalism toward microbe j. If ′ ( ) is positive but ′ ( ) is negative, this suggests that microbe j offers altruism toward microbe j or microbe j triggers parasitism or predation on microbe j. Taken together, we reconstruct a fully informative network that encapsulate bidirectional, signed, and weighted microbial interactions.
If a microbe has more links than the average of all microbes (i.e., connectivity) within a network, this microbe is regarded as a keystone microbe. A keystone species plays a disproportionately large role in maintaining the structure of an ecological community relative to its abundance [11,31]. Such species affect many other organisms in an ecosystem and help to determine the types and numbers of various other species in the community. The determination of a keystone species is based on the total number of links with it. Yet, a link can be active or outgoing in which case a species affects other species through promotion or inhibition. A link can also be passive or incoming, which means that a species is promoted or inhibited by other species. If a microbe has more outgoing links compared to its incoming links, this microbe is regarded as being "social." If the number of incoming links of a microbe is larger than that of its outgoing links, this microbe is "subordinate." A microbe can be viewed as a "leader," if the number of its outgoing links outperforms the average of all microbes. A "solitary" microbe is one that has fewer links, both outgoing and incoming. The classification of microbes into social, subordinate, leader, and solitary types can help researchers not only better understand the role of each microbe, but also exacerbate the mass or signals that are essential for relevant interaction patterns.

Validating and testing the qdGM
We show how qdGM can be used for network reconstruction in practice. As the proof of concept, we use two examples, human gut microbiota and microbiota residing in bodies of different animal species, although qdGM has a more widespread use in microbial research.
We also perform a simulation study to statistically validate its usefulness and utility.

Networking the gut microbiota
Our model was used to analyze a gut microbial abundance dataset containing eight phyla, 23 classes, 28 orders, 50 families, and 101 genera collected from 127 heathy hosts from a founder, the Hutterites [32]. These hosts include different sexes, different ages, and two treatments, use vs. no use of antibiotic intervention, of whom 93 were measured in winter, 91 in the following summer, and 57 in both. We define a measurement as a sample and, thus, this study contains 184 samples. We first reconstructed 184 sample-specific (i.e., dynamic) networks at the genus level based on the context-agnostic model of equation (1) and then test whether these networks should be classified into season-, sex-, age-, and treatment-dependent microbial networks using the context-specific model of equation (2).
We calculated the log-likelihood ratio of microbial abundance (817.73) under the seasonspecific and season-agnostic model, which is larger than the critical threshold at the 5% Paenibacillus, and Acetivibrio. There are some overlaps in secondary keystones between winter and summer networks.
As observed above, microbial networks display season-specific differences in terms of their organization rather than structure (Fig. 1). This statement is supported by keystoneness as a quantitative descriptor of the organization of microbial networks. Keystoneness contributes considerably to season-dependent difference (Fig. 1B). In the winter network, only a small portion of genera display a measureable capacity to link with other microbes through outgoing keystoneness (primary keystones) and incoming keystoneness (secondary keystones). Yet, keystoneness increases dramatically from winter to summer for a majority of genera, including primary keystones, secondary keystones, and non-keystones. Only a few genera, especially Veillonella, display a decrease from winter to summer High variability was observed in keystoneness between different hosts in both winter and summer, although, for most genera, inter-host variability is much more pronounced in summer than in winter (Fig. 2), suggesting that keystoneness can serve as a major determinant of interpersonal variability in gut-microbial interaction networks. Taken together, our findings suggest that, on average, gut-microbial community assembly responds to seasonal change not through altering the overall architecture of its gut-microbial interactions, rather through quantitatively changing the role of individual microbes, especially those highly social microbes, in microbial network organization.
We also compared sex-, age-, and treatment-specific microbial networks. We did not found qualitative differences of network architecture between different sexes, ages (young, ages = 19-30 years, vs old, ages > 60 years), and treatments (use vs. or no use of antibiotics), but these context-specific network differences are displayed in network organization and can be quantitatively explained especially by keystoneness ( Fig. S1 -S3).
We decomposed the overall abundance of each genus into its underlying endogenous component (determined by the intrinsic capacity of this genus) and exogenous component (affected by ecological interactions between other genera and this genus) and visualized how these two components contribute differently to microbial abundance across samples (Fig. 4).
The strength of promotion and inhibition increases exponentially with HI, despite genusdependent difference in increase slope. Together, the influences of these genera make Actinomyces's overall abundance higher than its endogenous component. However, in summer, the endogenous component of Actinomyces is higher than its overall level because of a considerably high promotion extent by genus Paenibacillus. Also, the types of genera that promote or inhibit Actinomyces change from winter to summer. In winter, genus Bacteroides increases its abundance with HI, a process influenced by numerous genera, of which two act through promotion and six through inhibition. Yet, in summer, its endogenous abundance decreases with HI, but with a positive promotion by Lactonifactor which makes Bacteroides outperform its endogenous component. These decomposition analysis provides a global and detailed picture of genus-genus interactions, their signs, magnitudes, and contextdependent changes, within the gut microbiota.

Networking the impact of microbiomes on animal evolution
More recently, a new concept of hologenome (host genome plus microbiome) has been proposed to postulate that the holobiont (host plus symbionts) evolves through natural selection because microbial symbionts can be transmitted from parent to offspring via cytoplasmic inheritance, coprophagy, direct contact during and after birth, and the environment [33,34]. This is contrast to the traditional Darwinian theory, stating that hosts are the units acted upon by natural selection [35]. Increasing studies have demonstrated that microbial symbionts contribute to the anatomy, physiology, development, innate and adaptive immunity, and behavior and finally also to genetic variation and the origin and evolution of species [36,37]. Given these findings, it is of great interest and importance to study how microbial interactions co-evolve with hosts.
We analyzed gut bacterial 16S rRNA gene sequence data collected on 265 specimens from 64 species [38]. The bacterial lineages collected were placed into phylogenetic taxa, approximately including 65 families that are common to all animal species. By modeling the total number of bacterial lineages per gut sample as a function of animal mass, a positive association was observed with a slope of 0.053 (P < 0.05), paralleling the known species-area relationship theory derived from island biogeography, consistent with the second law of thermodynamics [39,40]. We modified the system of qdODEs in equation (1) by replacing the habitat index derivative with the animal mass derivative and used these equations to reconstruct specimen-specific microbial networks which are further converted into speciesspecific networks. Figure 4 provides the examples of three animal species, bedbug (small size), hermitcrab (moderate size), and whale (large size). We found that the architecture of species-specific microbial network is predominated by commensalism and antagonism, together accounting for almost all links of microbial networks. The connectivity of microbiome networks is proportional to animal size; network architecture is more complex on whale than on hermitcrab as well as on hermitcrab than on bedbug. Family Oxalobacteraceae is a microbe found to serve as a primary keystone in microbial networks on all animals, regardless of their size. Yet, families Verrucomicrobiaceae and Enterobacteriaceae add to a list of keystones in small-size bedbug mostly by receiving incoming links, whereas moderatesize hermitcrab and large-size whale have family Helicobacteraceae as a secondary keystone through outgoing links although, for the latter, many tertiary keystones are identified.
Network properties, including the number of links for each bacterial family, changes dramatically among animal species, which may also be a determinant of animal adaptation to ecological environments.

Computer simulation
Results from simulation studies suggest that the model can capture a general trend of how each microbe interacts with others to form a network when the data contain a modest sample size (n = 50) and modest measurement error variance. The estimated networks approach the true networks when sample size increases and/or the error variance decreases. It appears that a sample size of 250 can reasonably well estimate a 100-node network from our model for error variance of intermediate size. Figure S4 illustrates the comparison of estimated and true abundance curves as a function of habitat index for two microbes chosen from a pool of 100 microbes. A sample size of 50 and modest error variance can estimate and decompose the overall abundance curve of a microbe into its endogenous and exogenous components, but with a fairly large bias. To obtain an acceptable curve estimate, a sample size of 250 is needed with modest error variance. The bias of estimation due to unpredictable large error variance can be compensated with a larger sample size.

Discussion
Although microbial community assembly exists as a highly complex, heterogeneous, and dynamic system, its underlying inner workings may still obey some certain fundamental principles that guide how microbes interact with each other to achieve maximum inclusive fitness [20]. No single theory can effectively illustrate the overall picture of microbial communications and interactions and reveal the detailed mechanisms underlying these interactions. We propose a quasi-dynamic game model (qdGM) by combining elements of biogeographic ecology theory and evolutionary game theory through statistical reasoning, aimed to unveil a basic rule for gut-microbial diversity, structure, and function. One main merit of qdGM lies in its capacity to pack steady-state abundance data into dynamic and mechanistic networks that encode and quantify all possible patterns of microbe-microbe interactions at play through bidirectional, signed, and weighted links. In the past, such fully informative networks only could be reconstructed from dense temporal data, although this type of data may be unavailable and less informative in gut microbiota studies.
In a population, the maximum payoff achieved by an organism may be quickly counteracted by the strategies of other organisms which also strive to maximize their payoffs [37]. This process will ultimately make the whole population achieve an optimal (equilibrium) payoff [41], which leads us to establish a statistical foundation for reconstructing a stable microbial network through maximizing the likelihood of the abundance data of all microbes [42]. The qdGM provides a way for quantitatively inferring such microbial networks, but not relying on temporal data. If the two microbes are mutualistic, this implies that they produce certain designer microbial consortia study more efficient and predictive.
Existing approaches can only identify a single overall microbial network from a large number of samples, failing to characterize network variability. The most remarkable advantage of qdGM is that it can reconstruct sample-specific networks (i.e., mobile networks) and convert them into context-specific networks through which to test how microbial interactions alter from treatment to treatment and vary along spatiotemporal gradients. Reanalyzing Davenport et al.'s [32] Hutterites data, we found that the overall structure of gut-microbiome networks is predominated by commensalistic or amensalistic interactions. This phenomenon was also observed in animal species-specific microbial networks. The overwhelming existence of commensalism or amensalism is consistent with widely identified cyclic synergism or antagonism (e.g., cyclic dominance) in nature, guided by the rock-paper-scissors game [44].
In winter's gut microbiota, we identified a few leaders, such as genera Prevotella, Veillonella, and Victivallis, which manipulate many other microbes and play a pivotal role in mediating the organization and function of microbial networks. Over a half of microbes in the gut are solitary because they only receive one or two incoming links from other microbes.
Although gut microbial networks do not make qualitative structural changes in response to seasonal change, sex difference, aging, and antibiotic treatment in their structure, they respond to these environmental or developmental changes in a quantitative way through The application of the qdGM to analyzing microbial data from different taxa of animals show its promising utility to reveal the evolution and ecology of holobionts, a new concept that has received increasing attention in modern biology [36,37]. The new model can not only identify which microbes drive animal evolution and ecology, but also reveal how different microbes work together to shape evolutionary processes. Based on data reanalysis of microbial abundance from 64 animals, we found that the role of most microbes as a mediator of networks change from animal to animal. As a primary keystone, family Oxalobacteraceae affects microbial networks across all animals of different sizes, but its role is played differently, depending animal size; by both incoming and outgoing links in small-and moderate-size animals but only by outgoing in large-size animals. We found a clear trend, i.e., microbial networks have greater connectivity for large than small animals. Taken together, the qdGM is armed with a capacity to identify key microbial indicators that may have played an important role in mediating animal energy expenditure and gain and, further, animals' adaptation and evolution over environmental change.
Results from simulation studies show that the qdGM has favorable statistical properties, and is expected to produce reasonably accurate microbial network with a modest sample size. As with all network inference methods, however, we cannot expect the networks inferred from the qdGM to cover all issues typical of microbial communities. The key step of the network inference is to solve a system of ODEs. Currently, we implemented the widely used fourthorder Runge-Kutta algorithm to estimate equation parameters. More robust algorithms, especially those with a capacity to handle high-or even ultrahigh-dimensional systems of differential equations and tackle errors due to experimental uncertainty or mis-assignment of sequencing reads into operational taxonomic units through stochastic differential equation modeling [11], are needed. The qdGM focuses on a set of homogeneous samples, but it can be readily extended to assemble microbial data from different sources with various perturbations [45]. With these and other modifications, the qdGM could potentially provide a powerful approach for inferring biologically meaningful microbial networks that contribute substantially to translational medicine and evolutionary research.

Conclusions
Networks are central to the functionality of complex systems. Existing approaches can only reconstruct a single aggregate microbial network from a number of samples, failing to characterize population heterogeneity. Furthermore, identifying the full information of networks that encapsulate bidirectional, signed, and weighted microbe-microbe interactiomes critically relies on the availability of temporal or perturbed data, which are difficult to collect or not informative for biological modelling. We develop a statistical model that can recover sample-specific informative networks from steady-state data. We integrate allometric scaling theory and evolutionary game theory to quantify sample-specific changes of microbial abundance and interactions by deriving a system of quasi-dynamic ordinary differential equations (qdODEs). We integrate functional clustering and variable selection to identify distinct network communities from high-or ultrahigh-dimensional microbiome data and the most significant links a given microbe receive to form sparse networks. We incorporate a maximum likelihood approach for estimating qdODE parameters that constitute the stability of microbial networks. We formulate a statistical platform for comparing and testing how microbial networks evolve over a phylogenetic clade of hosts, how they change across spatiotemporal gradients, and how they establish causal interrelationships with biogeochemical processes. We derive the mathematical properties of qdODE solving and perform extensive computer simulation studies to investigate the statistical behavior of microbial networks. The networks inferred from model can reveal the full information of microbial interactions, interrogate their ecological properties, and serve as an approach for visualizing how microbial communities assemble, stabilize, and evolve in a range of biological and biotechnological contexts.
Based on the structure of qdODEs in equation (1), this microbe's habitat index-varying abundance level can be described by a multiple regression model, expressed as In equation (4), () and ′ () are the habitat index-varying independent and dependent abundance of microbe j, whose derivatives are () and ′ () of equation (1), respectively, and ( ) is the residual error of microbe j at sample i, obeying a multivariate normal distribution with mean vector 0 and sample-dependent covariance matrix for microbe j. We assume that the residual errors of microbial abundance are independent among samples so that is structured as = 2 where 2 is the residual variance of microbe j at the same sample and In is the identity matrix. In equation (

Variable selection
In statistics, it is extremely difficult to model and estimate the effects of all other microbes on a focal microbe when the number of microbes is large. Indeed, this rarely happens in biology.
For these reasons, we need to choose a small set of the most significant microbes that affect a focal microbe. To do so, we need to choose appropriate functions to fit () and ′ () for equation (4) where the observed abundance level of microbe j is treated as a response and the observed abundance levels of all other microbes treated as predictors. Although their biological meanings can be hardly validated, non-parametric approaches, such as B-spline or Legendre Orthogonal Polynomials (LOP), can be statistically effective and efficient in fitting such functions through their linearization. As a first step toward network reconstruction, we incorporate the LOP into model (5) for the following variable selection step.

Group LASSO:
Considering a focal microbe j as a response, we use nonparametric independent parameters aj to fit its (). As predictors, (m -1) microbes contribute to microbe j's dependent component through unknown nonparametric unknown nonparametric dependent parameters j = (j1,…,j(j-1),i(j+1),…,jm). Thus, we have m -1 groups of dependent parameters that reflects the influence of other microbes on the focal microbe. We implemented group LASSO [46] to select those nonzero groups. The group LASSO estimators of dependent parameters, denoted as ̇ = (j1, …, j(j-1), j(j+1), …, ̇) , where dj ≪ m is the number of the most significant microbes that interact with microbe j. It can be obtained by minimizing the following penalized weighted least-square criterion,

Adaptive group LASSO:
In group LASSO, penalty parameters of each group are treated equally, without considering the relative importance of different groups. It has been recognized from traditional linear regression analysis that the over-penalization of parameters for some predictors may reduce the efficiency of parameter estimation and the continuity of variable selection [46]. To overcome this limit of group LASSO, Wang and Leng [47] integrated it with adaptive LASSO to create adaptive group LASSO. This integrative approach selects significant groups by weighted penalty parameters. Weight wj|j is obtained as ‖ | ′ ‖ 2 −1 if ‖ | ′ ‖ 2 > 0 and , otherwise. The adaptive group LASSO estimators of dependent parameters are obtained by minimizing the penalized weighted least-square criterion as follow: where 2j is a penalty parameter determined by BIC or extended BIC.

Likelihood formulation
Through variable selection, we change the number of incoming links for a microbe j from m to dj. A number of approaches, including non-linear least-squares and maximum likelihood, can be implemented to solve the shrunk qdODEs from equation (1). Since we argue that microbial community assembly tends to reach its maximum overall payoff guided by evolutionary game theory, a maximum likelihood approach that is founded on the most probable existence of all microbes is chosen for our qdODE solving. Let = ( ; ) denote the parameters that explain the regression model. The likelihood function of given the abundance data is written as where the number of microbes that regulate microbe j has been reduced to dj. Now, we strive to find a biologically relevant approach for fitting () and ′ () and a statistically robust approach for modeling .  (Fig. S5A). The abundance of genus Bacteroides follows a reverse pattern of change, i.e., increasing exponentially with niche index ( = 3.12; P < 0.001) (Fig. S5B). We found the random distribution of residuals independent from the power equation fitting of these two genera across niche index (Fig. S6).
In reanalyzing Sherrill-Mix et al.'s [38] data, a similar phenomenon was also observed ( and an LOP-based nonparametric approach to fit ′ (). We will model the covariance where is the sample-dependent covariance matrix of microbe j, and 1 2 is the sampledependent covariance matrix between microbes j1 and j2. 1

is structured as
Under mean-covariance structure modeling by equations (9) and (10), model parameters = ( ; ) become model parameters ′ = [Θ , Θ ′ ( = 1 , … , , ′ = 1, … , − 1, + 1, … , ); 2 , 1 2 ( 1 ≠ 2 = 1, … , )], whose optimal solution can be obtained, by maximizing the likelihood (8), as Intuitively, this maximization implies an optimal topological structure and organization by which microbes interact with each other to maximize the overall abundance of microbial community assemble as a whole. The convex optimization formulation under equation (11) ensures the stability and sparsity of the network reconstructed from qdODEs. Since no constraints are given on the number of outgoing links, the resulting network can be highdimensional and omnidirectional.
If n samples include C contexts, we need to test whether if the data can be better modeled according to different contexts. In this case, we can still formulate a likelihood like equation ), (12) and the covariance matrix modelled similarly as above.
By plugging in the MLEs of mean vectors (9) and (12) into likelihood (9), we obtain the likelihood values L1 (assuming that there are no differences among C contexts) and L2 (assuming that there are difference among C contexts), respectively. We further estimate the log-likelihood ratio, as a statistic used to test if n samples should be sorted into C contexts. By reshuffling n samples randomly into C groups, we calculate the LR value. If this permutation procedure is repeated 1000 times, we obtain the 95th percentile from 1000 LR values and use it as a critical threshold.

Quantifying the keystoneness of microbes
Since its proposal, the keystone species concept has been widely used to define entire ecosystems [31]. The qdODE model allows us to expand the definition of this concept.
Indeed, every species, as long as it links to any another, may play a bridging role in network organization. Based on the maximum likelihood estimates of (•) and ′ (•) under equation (11), we derive the parameter of keystoneness to assess the degree with which a microbe potentially behaves as the keystone. A similar strategy for keystoneness assessment was used by Berry and Widder [49], but our definition is quantitative and extendible to any cases. We define the keystoneness of a microbe as the sum of ratios of its dependent abundance to independent abundance over all microbes linked with it. Because of different behaviors of active and passive influences, outgoing links may contribute to keystoneness differently from incoming links. For a given microbe j, we define its outgoing keystoneness and incoming keystoneness, respectively, as where Dj and dj are the numbers of microbe j's outgoing links and incoming links, respectively. If the outgoing keystoneness of a microbe is larger than its incoming keystoneness, this microbe is regarded as a socially aggressive microbe. If a microbe has a smaller value of keystoneness, including outgoing and incoming, than the average of all microbes, it is regarded as a solitary microbe. As can be seen, keystoneness provide a quantitative classification of all microbes co-existing within the gut.

Computer simulation
The statistical properties of our qdODE network inference are empirically investigated through computer simulation studies. We simulated the abundance data of 100 microbes, yj(Wi) (j = 1, …, 100), for n samples, which vary with HI, Wi (i = 1, …, n), based on equation (5). The residual error of microbe j at sample i is assumed to obey a multivariate normal distribution with mean vector 0 and structured covariance matrix. We design different scenarios by changing sample size and variance, (1) n = 50, 2 = 1.5, (2) n = 50, 2 = 0.7, (3) n = 250, 2 = 1.5, and (4) n = 250, 2 = 0.7. We analyze the data simulated under each of these scenarios. By comparing the estimated and true HI-varying abundance curves, we assess the estimation precision of our new model.

Declarations
Ethics approval and consent to participate: Not applicable.

Consent for publication: Not applicable.
Availability of data and material: Fecal microbial abundance data https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0140301#ack The computer code is available at https://github.com/JiangLibo/microbialnetworks or may be requested from the corresponding author.