MDSINE: Microbial Dynamical Systems INference Engine for microbiome time-series analyses

Abstract

Predicting dynamics of host-microbial ecosystems is crucial for the rational design of bacteriotherapies. We present MDSINE, a suite of algorithms for inferring dynamical systems models from microbiome time-series data and predicting temporal behaviors. Using simulated data, we demonstrate that MDSINE significantly outperforms the existing inference method. We then show MDSINE's utility on two new gnotobiotic mice datasets, investigating infection with Clostridium difficile and an immune-modulatory probiotic. Using these datasets, we demonstrate new capabilities, including accurate forecasting of microbial dynamics, prediction of stable sub-communities that inhibit pathogen growth, and identification of bacteria most crucial to community integrity in response to perturbations.

Fig. 1

Schematic of the MDSINE software, which provides a comprehensive toolbox for dynamical systems analyses of microbiota time-series data. MDSINE implements a new algorithm, Bayesian Adaptive Penalized Counts Splines (BAPCS), for estimating microbial growth concentrations (trajectories) and their changes over time (gradients) from sequencing data; optionally, gradients can instead be estimated using our previously described first-order difference method. The software implements three new algorithms for dynamical systems inference: maximum likelihood constrained ridge regression (MLCRR), Bayesian adaptive lasso (BAL), and Bayesian variable selection (BVS). Our previously published method [8], the maximum likelihood unconstrained ridge regression algorithm (MLRR), is also implemented in MDSINE for comparison

Fig. 2

New inference algorithms in MDSINE outperform our previously published method on simulated data. Data were simulated to capture key features of real microbiome surveys, including noise and compositionality. Simulations assumed an underlying dynamical systems model with ten species observed over 30 days and an invading species at day 10. The number of time points sampled was varied between 8 and 27 to mimic common experimental designs and sequencing depths of 1000 or 25,000 reads were evaluated. Performance of the four MDSINE inference algorithms, maximum likelihood ridge regression (MLRR), maximum likelihood constrained ridge regression (MLCRR), Bayesian adaptive lasso (BAL), and Bayesian variable selection (BVS), were compared. Algorithm performance was assessed using four different metrics: root mean square error (RMSE) for microbial growth rates (a); RMSE for microbial interaction parameters (b); RMSE for prediction of microbe trajectories on held-out subjects given only initial microbe concentrations for the held-out subject (c); and area under the receiver operator curve (AUC ROC) for the underlying microbial interaction network (d). Lower RMSE values indicate superior performance, whereas higher AUC ROC values indicate superior performance

Fig. 3

Application of MDSINE to an experimental dataset evaluating the dynamics of Clostridium difficile infection in gnotobiotic mice. Germ-free mice were pre-colonized with the GnotoComplex microflora, a mixture of human commensal bacterial type strains chosen to capture the phylogenetic diversity and key physiologic capabilities of a native gut microflora. After the commensal microbiota were allowed to establish for 28 days, mice were infected with C. difficile spores and monitored for an additional 28 days. Throughout the experiment, 26 fecal samples per mouse were collected and interrogated via high-throughput 16S rRNA sequencing to determine abundances of species and 16S rRNA qPCR using universal primers to estimate the total bacterial biomass present. a Predicted directed microbe–microbe interaction network. Edge thickness denotes the magnitude of the evidence favoring presence of the interaction (Bayes factors [35]); only edges with strong evidence (Bayes factor ≥10) are displayed. b Predicted stable combinations of strains for each possible size of sub-community that optimally inhibit C. difficile colonization. Each row depicts the sub-community (combination of commensal strains) of a given size that is predicted to stably colonize the gut in the absence of the pathogen and is predicted to maximally inhibit C. difficile infection at the experimental end point (28 days). CDI median predicted median concentration of C. difficile at 28 days, CDI mad predicted absolute deviation of the median of the C. difficile concentration at 28 days

Fig. 4

Application of MDSINE to an experimental dataset evaluating stability of a probiotic cocktail in gnotobiotic mice. Germ-free mice were inoculated with 13 Clostridia strains from the VE-202 cocktail, a mixture of bacteria previously shown to be Treg inducers [24]. Over the 9-week experiment, mice were alternated between a standard high-fiber diet and a low-fiber dietary perturbation. We collected 56 fecal samples per mouse and these were interrogated using strain-specific qPCR primers to estimate strain concentrations. a Example of inferred growth and interaction parameters and their variability. The top grid displays mean parameter estimates and the bottom grid displays standard deviations of parameter estimates. Strains are ordered by their mean estimated growth rates on the standard diet. Pert. perturbation effect, St. strain. b Example forecasts of microbial concentration trajectories. Forecasts were obtained using a hold-one-subject-out procedure. Briefly, MDSINE was run on all data from all but one of the mice (the held-out subject) and model parameters were inferred. Using the inferred model parameters (including for the perturbation) and the measured concentrations of the microbiota at an initial time point for the held-out mouse, the trajectories of the microbiota for the held-out mouse were then forecast for all the remaining time points; the procedure was repeated for each mouse. Solid lines denote predicted trajectories and circles denote data. c Keystoneness analyses for high-fiber (top) and low-fiber (bottom) diets. Rows represent all possible stable states in which each strain has been removed in turn and the others are present (if that configuration is stable). The grid displays predicted steady-state concentrations of strains (log₁₀ ng strain DNA/μg total fecal DNA), with white entries indicating absent strains. Ky keystoneness, a measure assessing the marginal predicted quantitative effect of removing each strain from the full community, with larger values indicating greater effects on the overall ecosystem