As illustrated in Figure 2, ARZIMM can be considered as a two-part model which comprises a logistic component and an autoregressive component. To address zero inflation, we consider the zero-inflated mixture model because it assumes both sampling zeros (due to the low sequencing depth) and structural zeros (being truly absent) exist in the data. Specifically, the logistic component models the structure zeros of taxa in the samples, and the autoregressive component models the non-structure-zero abundances of the taxa under the assumption that the changes in abundances from time t − 1 to time t depend only on the observed abundances at time   t − 1 and other time-independent covariates, and the observed abundances before time t − 1 have no direct effect. Since the goal of ARZIMM is to characterize microbial interactions and community stability during a short period after a strong external perturbation like the antibiotic usage in our MIME study, we assume there are no other time-dependent factors exist to affect the microbial stability.

Let Y i m t   denote the observed absolute abundance of bacterial taxon m   ( m = 1 , ... , M ) for subject i   at time t   ( i = 1 ,   2 , ... , n ,   t = 1 , ... , T i ) , and we model Y i m t with a conditional mixture distribution as follow: Y i m t | ν i ( t − 1 ) ∼ { 0 p i m F ( y i m t | ν i ( t − 1 ) ; θ i t m ,   ϕ m ) 1 − p i m (1) where ν i ( t − 1 ) represents all information that is known at time ( t − 1 ) for individual i , including the observed absolute abundance Y i m ( t − 1 ) and later defined coviariates W i and Z i . The parameter p i m represents the probability of the observation Y i m t being structural zero and is assumed time independent. Furthermore, F is assumed to be an exponential dispersion family distribution with the canonical parameter θ i m t and the dispersion parameter ϕ m . Both Poisson and negative binomial (NB) distributions can be used as to model absolute abundance. Below we illustrate the detailed modelling using Poisson model.

The mixture probability parameters p i = ( p i 1 ,…,   p i M ) ′ are modeled by the logistic regression: log ⁡ i t ( p i ) = A W i + a i (2) where W i = ( 1 , w i 1 , … , w i l ) ′ consists of intercept and l time independent covariates for individual i , the parameter A = ( A 1 , … , A M ) ′ is an M × ( l + 1 ) matrix whose elements A m j is the effect of covariate j on the zero proportion of taxon m . a i = ( a i 1 , … ,   a i M ) ′ is an M × 1 vector of random intercepts to model the within-subject heterogeneity of being zero for individual i and has the joint multivariate normal distribution N ( 0 , Σ a ) .

The canonical parameters for Poisson distribution is θ i m t = log E ( Y i m t ) .     We introduce the auto-regressive model by relating   θ i t = ( θ i 1 t ,   … ,   θ i M t ) ′ to the i t h individual’s observed log-transformed absolute abundance vector at time t − 1 :   Y ˜ i ( t − 1 ) = ( log ( Y i 1 ( t − 1 ) + 1 ) , … , log ( Y i M ( t − 1 ) + 1 ) ) ′ (where the pseudo count 1 is added to avoid the undefined logarithm when the absolute abudance is zero), and Z i = ( 1 ,   Z i 1 , … , Z i q ) ′ ,   the intercept and q time-independent covariates of individual i by θ i t | Y ˜ i ( t − 1 ) = B Y ˜ i ( t − 1 ) + C Z i + η i (3) where B is an M × M matrix whose element B m j gives the effect of the abundance of taxon j on the growth rate of taxon m , C is an M × ( q + 1 ) matrix whose element C m j gives the effect of covariate j on taxon m , and η i = ( η i 1 , ... , η i M ) ′   is time-independent random intercepts. Note that, as an autoregressive model, η i is correlated with the fixed effect Y ˜ i ( t − 1 ) and this dependency can be tracked all the way back to the initial observation Y ˜ i 0 . Because the standard random effects model has assumption that the random effects are independent to the other covariates in the model, in order to derive the random effect type maximum likelihood (ML) estimators, we use the Chamberlain type projections (Chamberlain, 1982) to get around this correlation. Specifically, we project η i onto the time 0 observations Y ˜ i 0   by:     η i = Π Y ˜ i 0 + b i (4) where Π is an M × M matrix with diag( Π ) = ( π 1 ,   … ,   π M ) ′ and off-diagonal components being zero. The components of Π represent how much variation in η i   is due to the dependence on subject i ’s initial value Y ˜ i 0 . b i = ( b i 1 , ... , b i M ) ′ is an M × 1 vector, representing the independent subject-specific random effect and follows a joint multivariate normal distribution N ( 0 , Σ b ) .

In the model, our primary interest is to estimate matrix B , which measures the strengths of interactions between taxa. For a microbial community with a given number of species, its stability or dynamics status depends on the changes in the species’ population growth rates due to perturbation, which immediately cause the changes in the population growth rates of other species via species-species interactions (Ives et al., 2000). Interaction between species can be viewed as a filter that amplifies the variability in species’ population growth rates caused by perturbation.

Note that we choose Poisson distribution because of its nice stationary distribution property in the autoregressive model which is crucial for our following stability investigation. To deal with the over-dispersion of microbiome data, we implemented the quasi-Poisson model (Ver Hoef and Boveng, 2007) in the simulation and real data analysis.

To define the joint likelihood of the longitudinal microbial absolute abundance data Y i t , we assume that the vector of time independent random effects c i = ( a i ′ , b i ′ ) ′ underlies both the zero and autoregressive generative processes and these random effects account for the within-subject group heterogeneity in the multivariate logistic component and the multivariate autoregressive component. Denote   D = ( B , C ) = ( D 1 , ... , D M ) ′ , ϕ = ( ϕ 1 , ... , ϕ M ) ′ , and 1 [ ⋅ ] as the indication function that when [ · ]   meets, 1 [ ⋅ ]   = 1 , otherwise, 1 [ ⋅ ]   = 0 . Formally, we have the joint likelihood function as: ℒ ( D , A , Π , ϕ , σ ) = ∏ i = 1 n ∫ { [ ∏ t = 1 t i ∏ m = 1 M f y ( y i t m | θ i t m ( b i m ) , ϕ m , p i m ( a i m ) ) ] g ( c i | Σ ( σ ) ) } d c i (5) where f y is the conditional probability density function and given as f y ( y i t m | θ i t m ( b i m ) , ϕ m , p i m ( a i m ) ) = [ p i m + ( 1 − p i m ) f ( y i t m = 0 | θ i t m , ϕ m ) ] 1 [ y i t m = 0 ] × [ ( 1 − p i m ) f ( y i t m | θ i t m , ϕ m ) ] 1 [ y i t m ≠ 0 ] (6)

The function g ( c i | Σ ( σ ) ) is the joint distribution of c i , and Σ ( σ ) = [ Σ a Σ a b Σ a b Σ b ] represents the corresponding 2 M × 2 M covariance matrix, where σ accounts for all unique non-zero elements of Σ . For the model and computational simplicity, we assume C o v ( a i , b i ) = Σ a b = 0 , i.e. a i and b i are independent.

Assuming that the true underlying fixed effects A and D are sparse, we advocate a Lasso-type approach, which adds an ℓ 1 -penalty for the fixed-effects to the likelihood function. Thus, we consider the following objective function:   Q = − 2 ⁡ log ℒ + ∑ m = 1 M [ μ 1 m | | D m | | 1 + μ 2 m | | A m | | 1 ] . (7)

Maximization of the penalized log-likelihood function corresponding to Eq. 7 with respect to ( D , A , Π , ϕ , σ ) is a computationally challenging task. This is mainly because both integrals with respect to the random effects and the zero-inflated structure do not have analytical solutions. Following the conventional methods, we propose to implement a Laplace approximation on the integral of random effects in Eq. 7 and use the Expectation-Maximization (EM) algorithm to calculate the expectation and compute parameters iteratively, in which the label of zero is treated as “missing data”. The tuning parameters are selected using Bayesian information criterion (BIC).

The existence of a stationary distribution has been investigated for the log-linear Poisson auto-regression model based on the perturbation technique (Fokianos and Tjøstheim, 2011). Here, we prove the existence of a stationary distribution of a zero-inflated Poisson mixed-effect auto-regression model in Theorem 1 utilizing the theory of Markov chains which has been proposed to prove the existence of a stationary distribution of a general class of time series count models (Douc et al., 2013). The detailed proof is provided in the Supplementary Material Section S3.

Theorem 1

Assuming that time-independent parameters η i and p i are known, if all eigenvalues of matrix B lie inside the unit circle, a strict-sense stationary ergodic process { Y i t } t ∈ N will exist, where N denotes the set of natural numbers.

With this Theorem, we can first show that for a microbial community, its dynamic process { Y i t } t ∈ N has a stationary distribution by proving that all eigenvalues of matrix B lie inside the unit circle. Then, following Ives et al. (Ives et al., 2003), we consider the return rate and reactivity as two stability measures based on the variability of the stationary distribution for MAR (1) model. Specifically, return rate depends on the rate at which the perturbed microbial community approaches the stationary distribution and reactivity, and assesses how strongly population-level microbiome abundances are pulled towards the mean of the stationary distribution. Both are bounded by the largest eigenvalue of B , denoted by max ( λ B ) . In general, a smaller max ( λ B ) indicates the perturbed microbial community approaches its stationary distribution faster, or a system is less reactive, then the microbial community is more stable. The detailed proof is deferred in the Supplementary Material Section S3.2.

Based on the theory in Ives et al. (Ives et al., 2003), for a community with multiple species, the covariance matrix of the stationary distribution depends on the covariance matrix of the process error and the interactions between species captured in the matrix B . As illustrated in Figure 1, when the external perturbation(blue arrow) acts on the community, the ball(microbial community) sitting in a deep bowl in state 2 which represents a relatively stable system, will return to its stationary state faster than the ball sitting in a shallow bowl in state 1 which represents a less stable system. In a stable system, the variance of stationary distribution is only slightly greater than the variance of process error and the variance of species interaction is small. In contrast, in a less stable system, the species interaction will amplify the environmental variance and create large variance in the stationary distribution, therefore the variance of species interaction is large, assuming the process errors are similar in the compared two states. Thus, the difference between the variances of stationary distribution of different communities can be attributed to species interactions. The smaller of the variance of matrix B , the more stable of the study microbial community.

We have conducted extensive simulation studies to evaluate the performance of ARZIMM in both model fitting and variable selection by comparing it with the competing methods: penalized Poisson auto-regression (Poisson), penalized log-normal multivariate auto-regression (MAR), and extended generalized Lotka-Volterra (gLV) equations using Bayesian algorithm (MDSINE) (Bucci et al., 2016b). The brief descriptions of these methods are provided in the Supplementary Material Section S2.

We applied ARZIMM methods to the MIME study. The MIME study is an ongoing randomized trial on 80 healthy volunteers with one control group (ctrl) and two antibiotic groups (amoxicillin, amx, and azithromycin, azm); antibiotics are provided for a 1-week period at the start of the trial. The main microbiome research goal of the MIME study is to evaluate the effects of antibiotics on microbial profiles at both the community and taxonomical levels. With ARIZMM, we propose a different perspective to evaluate the effect of antibiotics through the investigation of microbial interaction and community stability across groups. Because the clinical trial is still ongoing and only partial data are available, the following data analysis is done on a subset of MIME data including only 11 subjects who were randomized to two groups: 4 ctrls and 7 azms. The main purpose of this analysis is to illustrate how to use ARIZMM, not for the scientific conclusion. For each subject, we collected two baseline microbiome samples, three samples during the course of antibiotics, and five post-antibiotic samples. The gut microbiota of these individuals were profiled using 16S rRNA gene targeted sequencing on the Illumina MiSeq platform. To obtain the microbial absolute abundances, we multiplied the relative abundances of OTUs by the sample density 1.1 g/cm³ and the number of universal 16S rRNA per gram measured using qPCR (Stein et al., 2013a). In our analysis, samples that collected before treatment in both antibiotic groups were excluded. The abundances of taxa were agglomerated at the genus level and taxa were further filtered if 1) the average relative abundances over all samples are less than 0.1%, and 2) the taxa are presented in less than 5 samples within each group.

First, Figure 4A shows a comparison of the relative abundance (top panel) and the absolute abundances determined by quantitative sequencing (bottom panel) of the dominant bacterial genera in 99 fecal samples from 11 subjects (blocks) across seven to nine time points (shown from left to right within each block) of this preliminary dataset. It is evident that the relative abundance and absolute abundance data present different information about the microbial profiles, and that the total bacterial load changes over time for each subject (i.e., within each block). Thus it is essential to study the microbial interactions using the absolute abundance data.

Then, we evaluate the model fitting of the log-normal distribution [used in MAR(1)] and zero-inflated over-dispersed Poisson distribution (used in ARZIMM) on the available subset of MIME data using chi-square goodness of fit test at 5% significance level taxon by taxon. Out of 45 taxa in the control group, 1 and 44 of their absolute abundances were fitted well (p > 0.05) by log-normal distribution and zero-inflated over-dispersed Poisson distribution respectively. The log-normal distribution fails to fit the data well when microbial taxa’s absolute abundance data are left-skewed and sparse (two examples are illustrated in Figure 4B).

Next we demonstrate how to conduct inference for microbial interactions and community stability with ARZIMM on MIME data. First, we fit ARZIMM to ctrl and azm groups separately, adjusting for age, gender, and BMI, to get their estimated interaction matrix B ^ s. Table 2 reports the characteristics of microbial interaction matrix estimates B ^ s. Defining the interaction effect as informative if its B ^ m j ’s 95% bootstrap confidence interval (based on 100 bootstrap samples) does not contain zero, we identified 125 and 105 informative interactions, respectively, in azm and ctrl groups. Their interaction effects are illustrated using networks in Figure 5. With more informative interactions, the azm groups have bigger and more complex networks than the ctrl group (first row of Figure 5), while the control group has more large estimated interaction effects than those in azm group as showed in Table 2 and the last three rows of Figure 5. This observation indicates that the antibiotic treatment reduce the strength of the interactions among the taxa and create more variations with more weak interactions among taxa, thus reduce its stability. In the last row of Table 2, based on our stability theory we report the stability properties of the studied microbial communities. The ctrl group has the lower estimates of maximum eigenvalue squared 0.11 comparing to the azm group’s maximum eigenvalue squared 0.32, which indicates that the control microbial community is more stable than the antibiotic communities.

Figure 6 provides additional information on the network feature comparison between ctrl and azm groups. Figure 6A displays the distribution of the positive and negative informative interaction estimates separately. The ratios between the numbers of positive and negative interactions are both around 1:1 in two groups. Figure 6B presents the frequency distribution of vertex degree of all the taxa in each group and they are all skew to the right. In the figure, a vertex represents a taxon in a community and its vertex degree is the number of informative interaction effect it has with the other taxa. By defining average neighbor degree as the average number of a given taxon’s neighbor vertices’ degrees, Figure 6C shows that the average neighbor degree is negatively correlated with the vertex degree in azm antibiotic treated group, but not in the control group. This indicates that there may be a group of taxa interacting with each other actively in the antibiotic group. It would be interesting to identify such sub-community with additional effort.