(C) PLOS One
This story was originally published by PLOS One and is unaltered.
 . . . . . . . . . .



Extracting host-specific developmental signatures from longitudinal microbiome data [1]

['Balázs Erdős', 'Department Of Data Science', 'Knowledge Discovery', 'Simula Metropolitan Center For Digital Engineering', 'Oslo', 'Christos Chatzis', 'Faculty Of Technology', 'Art', 'Design', 'Oslo Metropolitan University']

Date: 2026-07 

Longitudinal microbiome studies provide critical insights into microbial community dynamics and their relation to host health. Tensor decompositions offer a powerful framework for the unsupervised analysis of such data, yielding interpretable low-dimensional temporal patterns. However, existing approaches based on the CANDECOMP/PARAFAC (CP) model assume common temporal dynamics for all subjects and therefore cannot capture subject-specific trajectories. To address this limitation, we introduce a novel analytical framework based on PARAFAC2 to explicitly model subject-specific variations, such as shifts and delays in temporal patterns. Through systematic comparisons on simulated and real-world datasets—including studies of infant gut maturation and dietary interventions—we demonstrate that PARAFAC2 outperforms CP in capturing subject-specific temporal trajectories, and enables the discovery of biologically relevant patterns that are overlooked by CP. Furthermore, we introduce replicability as a robust criterion for selecting the number of model components, ensuring that the extracted patterns are replicable.

Longitudinal microbiome datasets are complex, consisting of repeated high-dimensional compositions that track changes in microbial abundance across individuals over time. While tensor decompositions are powerful tools for unraveling structure in these data, standard models like CANDECOMP/PARAFAC (CP) impose a critical limitation: they assume that temporal dynamics are identical across all individuals. This assumption often fails to capture the heterogeneity of biological processes, such as the varying pace of gut microbiome maturation or distinct individual responses to dietary changes. In this work, we introduce a novel analytical framework leveraging the PARAFAC2 model to overcome these constraints. By explicitly modeling subject-specific temporal variations, our PARAFAC2-based approach allows for the detection of individual time shifts and delays. We validated this framework using both simulation and real-world cohorts, demonstrating its ability to recover personalized trajectories that CP obscures. Additionally, we implemented a robust criterion to guide model selection, ensuring that the discovered patterns are replicable.

Funding: We acknowledge all funding received by COPSAC, listed on www.copsac.com . Specific support for this study was provided by the Lundbeck Foundation https://lundbeckfonden.com/ (grant no. R269-2017-5, COPSYCH). M.A.R. is funded by the Novo Nordisk Foundation https://novonordiskfonden.dk (grant no. NNF21OC0068517). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Data Availability: Data and code underlying the simulations are available at https://github.com/blzserdos/parafac2_microbiome and archived at https://doi.org/10.5281/zenodo.17674126 . Shotgun metagenomic sequence data from the FARMM dataset (Tanes et al. 2021) were deposited under BioProject with accession code PRJNA675301. Processed data obtained from Ma et al. (2023) are included in the GitHub repository and Zenodo archive listed above. Individual-level clinical data from the COPSAC2010 cohort are not publicly available to protect participant privacy, in accordance with the Danish Data Protection Act and European Regulation 2016/679 of the European Parliament and of the Council (GDPR) that prohibit distribution even in pseudo-anonymized form. Data can be made available under a joint research collaboration by contacting COPSAC’s data protection officer ( administration@dbac.dk ).

In this work, we introduce a novel analytical framework based on PARAFAC2 for the analysis of longitudinal microbiome data, allowing us to model subject-specific temporal trajectories. We evaluate and compare CP and PARAFAC2 on both simulated and real-world longitudinal microbiome datasets, including studies on infant gut maturation and dietary interventions. We show that PARAFAC2 more effectively captures subject-specific temporal dynamics, which in turn enables the discovery of biologically relevant patterns that are overlooked by the more restrictive CP model. We further introduce and discuss the use of replicability analyses in model selection to improve the robustness of the extracted patterns.

Unlike the CP model, PARAFAC2 [ 17 ] offers a more flexible alternative by relaxing the strict assumption of multilinearity when analyzing multi-way data by letting the factors in one mode change across tensor slices. This allows it to model shape changes, such as the shifts and delays that characterize individualized temporal trajectories. The power of this flexible approach has been well established in other domains. In chemometrics, for instance, PARAFAC2 is used to resolve and quantify chemical compounds in mixtures where chromatographic profiles exhibit retention time shifts between runs [ 18 ]. In neuroscience, it successfully identified changes in brain connectivity consistent with prior pathological knowledge of Alzheimer’s disease [ 19 ] and captured subject-specific temporal profiles from the analysis of multi-subject functional neuroimaging data [ 20 ]. It has been applied to electronic health records to construct personalized disease trajectories, capturing phenotypic evolution across clinical encounters [ 21 ]. More recently, in single-cell genomics, PARAFAC2 has been used to integrate multi-condition single-cell measurements, separating condition-specific from cell-to-cell variation [ 22 ].

However, a critical limitation of CP-based methods is the assumption of shared time loadings across all subjects. This assumption restricts these methods to modeling only population-level dynamics and prevents them from capturing subject-specific temporal variation that is commonly observed in microbiome studies [ 13 ]. This limitation is particularly relevant when individual temporal variation may reflect important biological information. For instance, in studies of infant gut microbiome maturation, subjects may follow the same underlying maturation process but at different paces, representing biological heterogeneity or environmental exposures [ 14 ]. In intervention studies involving antibiotics or dietary changes, subjects often exhibit divergent responses in the timing and magnitude of microbial adaptation (even within the same treatment group), underscoring the need for models that can capture these personalized dynamics [ 15 , 16 ].

To address these limitations, tensor decomposition-based approaches have recently emerged for the analysis of longitudinal microbiome data. Tensor decompositions, which are extensions of matrix factorizations to multi-way arrays (i.e., higher-order tensors), factorize multi-way data into a set of interpretable components, with applications across diverse fields such as chemometrics, neuroscience, and computational phenotyping [ 3 – 5 ]. In microbiome research, existing tools such as compositional tensor factorization (CTF) [ 6 ], TCAM [ 7 ], parafac4microbiome [ 8 ], microTensor [ 9 ], and temporal tensor decomposition (TEMPTED) [ 10 ] have primarily relied on the CANDECOMP/PARAFAC (CP) decomposition [ 11 , 12 ]. This formulation rests on the biological assumption that the longitudinal variation observed across taxa, time, and subjects can be approximated by the sum of a small number of latent signatures (components), each given by the outer product of a taxa, a time, and a subject loading vector—corresponding, respectively, to a microbial signature, its temporal progression, and its expression across individuals.

The human microbiome exhibits complex temporal dynamics that vary substantially both within and across individuals, making longitudinal analysis crucial for understanding its evolution through developmental stages (e.g., through childhood) or following interventions, as well as its link to host health. Advances in high-throughput sequencing now allow for dense temporal sampling, offering new opportunities to link microbial dynamics with host phenotypes. However, analyzing such data poses distinct computational challenges that call for specialized methods [ 1 ]. Traditional dimensionality reduction methods, such as principal coordinates analysis (PCoA) [ 2 ], operate at the sample level and fail to account for the inherent temporal structure and within-subject correlations present in longitudinal data.

where , , and are the r-th columns of the ground truth factor matrices, and , , and are the corresponding model estimates. In the case of PARAFAC2, and refer to vectors produced by vertically concatenating and , respectively.

where is the data tensor, is its reconstruction based on the model. A model fit of 100% indicates that the model perfectly reconstructs the data, while a fit below 100% implies that some variation in the data remains unexplained.

Robustness to mode-1 fiber missingness. To evaluate the robustness of PARAFAC2 to mode-1 fiber missingness, we performed an ablation analysis on the FARMM dataset. Starting from the original tensor (11.8%; 53/450 missing mode-1 fibers), additional mode-1 fibers were randomly set to missing. Datasets were constructed at 20%, 30%, 40%, and 50% missingness, with 10 independent realizations per rate. The 3-component PARAFAC2 model that was selected in the FARMM data analysis was applied to each ablated dataset, with 20 random initializations. For each fit, FMS was computed against the solution obtained on the original data after optimal permutation alignment of components.

The ridge penalty was added to aid convergence [ 41 ], using a penalty parameter set to . Non-negativity was imposed to resolve the sign ambiguity of PARAFAC2. No non-negativity constraints were imposed on the taxa or time modes, allowing A and to take both positive and negative values, consistent with the signed nature of CLR-transformed data. Since C is non-negative, a positive taxa loading indicates that the taxon’s abundance tracks the temporal pattern in the component, while a negative loading indicates that it varies in the opposite direction. Due to missing samples, the FARMM dataset contains completely missing mode-1 fibers (i.e., all taxon relative abundances for subject k at time point j are missing). PARAFAC2 struggles with such scenarios, as can be chosen arbitrarily, provided it conforms to the imposed constraints. To improve the robustness of the recovered factor matrices, we used a smoothness constraint in the form of a graph Laplacian penalty (with penalty parameter ) on columns of the factor matrices to penalize large differences in neighboring time points [ 32 ]. To avoid convergence to local minima, each model was fit using multiple random initializations. After confirming the uniqueness of the solution, the solution with the best fit was selected for further analysis. Finally, to address the scaling ambiguity of CP and PARAFAC2, we scaled the columns of factor matrices to unit norm.

Before applying the tensor decompositions, data tensors were divided by their Frobenius norm. For fitting the CP models, we used an alternating least squares-based algorithm [ 37 ] implemented in Tensorly (v0.9.0) [ 38 ]. The PARAFAC2 models were fit using alternating optimization with the alternating direction method of multipliers (AO-ADMM) algorithm implemented in MatCoupLy (v0.1.6) [ 32 , 39 ]. In experiments with missing data, we used an expectation maximization (EM)-based approach with both CP and PARAFAC2 [ 38 , 40 ]. We applied the following constraints and regularization terms:

The Food and Resulting Microbial Metabolites (FARMM) study collected daily fecal samples over a 15-day period from 30 adult participants, evenly distributed across three diet groups: vegan, omnivore, and exclusive enteral nutrition (EEN) [ 16 ]. The study comprised three distinct phases: the dietary phase (days 1–5), the gut microbiota purge phase (days 6–8), and the recovery phase (days 9–15). All participants received antibiotic and polyethylene glycol (Abx/PEG) treatment during days 6–8. For a detailed description of the study design, sample collection, and sequencing procedures, we refer to the original publication [ 16 ]. Preprocessed shotgun metagenomics-based relative abundance data were obtained from [ 9 ]. Taxa detected in fewer than five samples with a relative abundance of at least were filtered out. Additionally, time point zero was removed because no subjects in the vegan group had samples at that time point. Subsequently, the CLR transformation was applied after adding a pseudo-count of 0.5. The data were arranged as a third-order tensor of size with modes taxa, time, and subjects, respectively. The data tensor contained 11.8% missing entries.

The Copenhagen Prospective Studies on Asthma in Childhood 2010 (COPSAC 2010 ) is a population-based birth cohort of 700 children recruited in pregnancy and followed prospectively with deep clinical phenotyping [ 33 ]. The data used here comprises fecal samples collected at 1 week, 1 month, 1 year, 4 years, and 6 years after birth and analyzed via 16S rRNA gene amplicon sequencing. Details of the laboratory workflow, sequencing and data processing have been described previously [ 34 , 35 ]. Additional covariates include mode of delivery (vaginal delivery/Cesarean section) and maternal intrapartum antibiotics (yes/no). Children with missing samples were excluded from the analysis. Species not detected in ≥ 10% of subjects in at least one time point were filtered out. To account for the compositional nature of the data, the centered log-ratio (CLR) transformation was applied, using a pseudo-count of 0.5 [ 36 ]. After pre-processing and data transformation, the data were arranged as a third-order tensor of size with modes taxa, time, and subjects, respectively. The data tensor contained no missing entries.

Top panel: ground truth factors including the taxa ( ), time ( ), and subjects ( ) loadings visualized after scaling the time loadings by the corresponding subject loading ( ). The microbial signature is present in subject k over time according to their scaled time loadings . Group differences characterized by a presence/absence of the microbial signature in the subject is indicated by dashed colored/continuous gray lines. Thick borders highlight the component 2 subpanels containing the patterns of interest. These comprise subject-specific time loadings for as denoted by the respective line colors. Middle and bottom panels: factors recovered by CP and PARAFAC2, respectively.

As a proof-of-concept, and in order to illustrate the use of PARAFAC2 to capture subject-specific temporal trends in longitudinal microbiome data, we constructed a simulated data tensor including variation in the temporal patterns across subjects. The simulation is intended as a controlled, illustrative example rather than a benchmark of realistic microbial dynamics, with the known ground truth allowing CP and PARAFAC2 to be compared directly with respect to their ability to recover subject-specific temporal patterns. The data were generated in a way that does not fully conform to the structural assumptions of either CP or PARAFAC2, to better represent real-world scenarios where the true data-generating process is unknown. The simulated data tensor with I = 36 microbes, J = 21 time points, and K = 12 subjects was constructed according to the ground truth factors shown in Fig 2 . Each of the three patterns, corresponding to the R = 3 components, describes the evolution of a microbial signature over time according to the temporal profile with strength in subject k. All three components are constructed to encode group-level differences between subsets of subjects through the subject loadings ( ). In addition, the second component includes subject-specific delays in the temporal profiles ( ), through variation introduced to the time loadings across a subset of subjects. In particular, the first pattern is the bell-shaped trend in time of taxa 1–12 and is unique to the group of subjects {1,2,6,7,11,12}. The third pattern is the dual-peaking temporal trend of microbes 19–30 and is unique to the group of subjects {7,8,9,10,11,12}. The second pattern is the subject-specific saturation of microbes 10–21 that is unique to each subject within the group {4,5,6,7,8,9}. Microbes 31–36 have low feature loadings in all subjects at all time points. Of interest are the between-group differences, with only one group of subjects exhibiting the taxa signatures in each pattern, as well as the between-subject differences that are specific to a group of subjects in component 2. Simulated tensors were constructed both in the absence and presence of noise. To introduce noise, we constructed a noisy tensor as

The number of subsets F depends on the characteristics of the dataset in use, such as the number of subjects K and the presence of relevant grouping labels. We consider an R-component model to be replicable if at least 90% of the pooled FMS AB (or FMS A and FMS C*B in the case of PARAFAC2) values exceed 0.9. To maximize explanatory power, we select the largest number of components that still yield a replicable model. Beyond replicability, we adopt an exploratory approach to assess biological relevance [ 28 ]. While PARAFAC2 comes with estimating many more parameters compared to CP, and thus an increased risk of overfitting, the replicability-based approach serves as a safeguard against this, as components reflecting idiosyncratic variation in a particular sample will fail to stably reappear across random subsamples, and are therefore flagged rather than interpreted.

The factor match score is calculated over the scaled time loadings to avoid overemphasizing poorly estimated time loadings, which can arise when the corresponding subject loading is close to zero [ 32 ]. Because of the permutation ambiguity in CP and PARAFAC2 decompositions, we first align the order of components to the optimal permutation before computing FMS. The FMS ranges from 0 to 1, with higher values indicating greater similarity. A visual illustration of the model selection using replicability is shown in S2 Fig . Specifically, we assess the replicability of an R-component CP and PARAFAC2 model as follows.

In the case of PARAFAC2, we compute FMS A using the relevant terms (i.e., taxa mode loadings ) of Eq 3 . Since the time mode factors are specific to the set of subjects in the data tensor, they can only be compared using the loading vectors corresponding to the set of subjects common to any two submodels. To do this, we define maps such as and from a position in the complete tensor to its positional index in the submodel. Let and be the vectors of vertically concatenated scaled time loadings of the subjects shared between submodels one and two:

where denotes the vector 2-norm, and are the r-th columns of the taxa and time factor matrices from the CP model of (first submodel), and are the corresponding columns from the CP model of (second submodel).

Selecting the appropriate number of model components R is crucial for ensuring that the resulting model is both biologically meaningful and valid. Several approaches have been proposed to determine the number of components [ 26 – 28 ]; nevertheless, this task remains an active topic of research. In this study, we determine R based on the replicability and interpretability of the extracted components [ 28 , 29 ]. Replicability refers to the consistency of identifying similar patterns across random subsamples of the dataset and extends the concept of split-half analysis [ 30 ]. We adopt a replicability-based model selection approach—by subsampling along the subject mode—to ensure that the extracted patterns are robust and representative of the underlying study population, as previously done to select the number of components for CP models [ 31 ]. Let be the set of all slice indices in the subject mode of the complete data tensor. As an example, consider the subsampled tensors, and , formed by stacking the slices and , respectively, where and are randomly selected, partially overlapping subsets of indices and is the cardinality of the set . To quantify the similarity between factors extracted from different subsampled tensors (e.g., and ), we use the factor match score (FMS). For two CP models, e.g., CP models of and , FMS over the taxa and time modes (FMS AB ) is computed as:

Taxa, time, and subject loadings of the CP and PARAFAC2 models are visualized in Fig 1 , either using all three modes (panel i), or in two modes (panel ii), by scaling the time loadings with their corresponding subject loading, referred to as scaled time loadings, i.e., .

where is the set defining the constant cross-product constraint of PARAFAC2, which ensures unique decomposition up to scaling and permutation ambiguities under certain conditions [ 25 ]. See S1 Fig for a geometric illustration of this constraint. The added flexibility leads to the additional challenge of sign ambiguity in PARAFAC2. Reformulating Eq (2) as , where , , and , we observe that may arbitrarily flip signs together with [ 17 ]. We resolve this ambiguity by imposing non-negativity on C [ 25 ]. Given a taxa by time by subjects tensor , PARAFAC2 relaxes the strict assumption of the CP model by allowing temporal patterns to vary across subjects k. In contrast, CP requires identical temporal patterns for all subjects, i.e., ( Fig 1 ). Therefore, the PARAFAC2 model may be better suited for datasets in which subjects exhibit heterogeneous temporal patterns.

PARAFAC2 decomposes into taxa factors , subject factors , and host-specific time factors , where K is the number of subjects. PARAFAC2 allows the factor matrices to vary across tensor slices under the constraint that their cross product is constant. This is more flexible than CP, which finds common time factors describing all subjects. Component r comprises taxa loadings characterizing a microbial signature (a set of co-varying taxa), time loadings in case of CP (or with PARAFAC2) that capture the temporal profile along which this signature is expressed, and subject loadings that act as per-subject weights indicating how strongly each individual exhibits the component. Example visualization of a component across its three modes (i) , and in two modes, after deriving scaled time loadings by scaling the time loadings with their corresponding subject loading, i.e., (ii) . Scaled time loadings in case of CP are calculated as .

where is an diagonal matrix with the kth row of C on the diagonal, and ⊤ denotes the matrix transpose. The CP model is unique up to permutation and scaling ambiguities under mild conditions [ 24 ], ensuring reliable interpretation of the resulting factors. More specifically, columns of the factor matrices can be permuted in the same way across modes (permutation ambiguity), and can be scaled as long as the component norm stays the same, e.g., (scaling ambiguity). Under these ambiguities, the interpretation of the CP factors remains unchanged. For a taxa by time by subjects tensor , each column of A (taxa loadings), i.e., , may reveal a microbial signature—a set of co-varying taxa. The corresponding temporal profile, describing how this signature is expressed over time, is captured in the columns of B (time loadings), i.e., . Columns of C (subject loadings), i.e., , act as per-subject weights indicating how strongly each individual exhibits the component ( Fig 1 ). In this way, the loadings are weights that quantify the contribution of each taxon, time point, or subject to a component.

The CP model [ 11 , 12 , 23 ] approximates a higher-order tensor as the sum of rank-one tensors. An R-component CP model of a third-order tensor is defined as where ∘ denotes the vector outer product, , , and correspond to the r-th column of factor matrices , , and , respectively. Using the CP model, the kth frontal slice of , denoted , is approximated as:

This study analyses data from two previously published longitudinal microbiome studies; no new data were collected. The COPSAC 2010 cohort study was conducted in accordance with the guiding principles of the Declaration of Helsinki and approved by the Local Ethics Committee of the Danish Capital Region (H-B-2008–093) and the Danish Data Protection Agency (2015-41-3696). Both parents provided written informed consent before enrollment. The FARMM study data were obtained in processed form from a publicly available source [ 9 ].

To further assess robustness to mode-1 fiber missingness, we performed an ablation analysis on the FARMM dataset. Starting from the original data (11.8% missingness; 53/450 fibers), additional mode-1 fibers were incrementally introduced to obtain datasets at 20%, 30%, 40%, and 50% missingness (10 random realizations per rate; see Materials and Methods). Factor recovery against the original-data solution, quantified by FMS, remained high across all rates (median FMS ≥ 0.94 at 50% missingness; S15 Fig ), and convergence on a representative 50% missing realization was comparable to the original-data fit ( S16 Fig ).

To illustrate the spectrum of uncertainty, we visualized the scaled time loadings from all relevant submodels for a selection of representative subjects from both cohorts ( Fig 5 ). Overall, the profiles showed good agreement across submodels, confirming their replicability at the subject level, even within the same exposure group. For example, panel (ii) shows the mean (± 2 standard deviations) for three subjects (A, B, C) in the VD - no IAP exposure group from the COPSAC 2010 dataset, all of which show high consistency (subject-specific average FMS C*B ≥ 0.99). Similarly, panel (iv) shows robust profiles for subjects 9003 and 9024 of the EEN study group from the FARMM dataset, with subject-specific average FMS C*B ≥ 0.99. The comparatively larger uncertainty observed for subject 9032 (subject-specific average FMS C*B = 0.97) reflects the underlying missingness pattern of the data, as this subject was missing samples on study days 1–2 and 7–9 (see Discussion). Results displaying the complete set of model components are shown in S13 Fig and S14 Fig , for the COPSAC 2010 and FARMM studies, respectively.

Panels (i, iii) show histograms of the average standard deviation, calculated for each subject by first finding the standard deviation of the scaled time loadings at each time point across submodels, then taking the average. Panels (ii, iv) show the mean (lines) ±2 standard deviations (shaded area) of scaled time loadings at each time point across submodels for three representative subjects from the VD - no IAP group in the COPSAC 2010 cohort (panel ii) and the EEN group in the FARMM (panel iv) study.

Beyond the global replicability (FMS C*B ) used in model selection, we also assessed the uncertainty of subject-specific temporal patterns. To do this, we first identified, for each subject, the specific submodels from the replicability analysis in which they were included (89/100 and 39/50 submodels per subject for COPSAC 2010 and FARMM, respectively). We then calculated the standard deviation of their scaled time loadings at each time point across the submodels and averaged the standard deviations to use as an indicator of the variance in their inferred profiles. The average standard deviations were generally an order of magnitude smaller than the scaled time loadings, indicating good agreement across submodels and robust recovery of individual temporal patterns ( Fig 5 , panels i, iii).

Both models identified broadly similar temporal trajectories of the microbial signatures, as reflected in the scaled time loadings; however, PARAFAC2 revealed subject-specific dynamics that were not apparent in the CP results. PARAFAC2 uncovered distinct microbial response to the exclusive enteral nutrition (EEN) diet compared to the vegan and omnivore diets in component 1, which was missed by CP ( Fig 4 , panel v). In particular, PARAFAC2 resolved a decrease in the abundance of the microbial signature during the dietary phase for the EEN group. This dynamic was not observed in the other groups. The improved resolution is further reflected in the corresponding subject loadings, where PARAFAC2 successfully distinguished between the EEN and vegan groups ( Fig 4 , panel vi). These findings align with those reported by Tanes et al. (2021) [ 16 ], highlighting the capability of PARAFAC2 to capture complex, individualized temporal patterns. Furthermore, while the CP model oversimplified the recovery phase as a constant difference across time, PARAFAC2 accurately modeled the evolving divergence between the EEN group and others ( S11 Fig , panels ii, v).

Shown are component 1 of the selected CP (panels i–iii) and PARAFAC2 (panels iv–vi) models, associated with the dietary (days 1–5) and antibiotics phases (days 6–8) of the trial. Panels (i, iv) display the 15 species with the largest absolute weights in each component, representing the extracted microbial signatures. Scaled time loadings (ii, v) at time points were averaged per study group. Shaded area indicate standard errors of group means. Panels (iii, vi) present subject loadings compared across exposure groups using the Mann–Whitney U test with Bonferroni correction.

The overall patterns of microbial features suggested a general concordance between the models across phases, as reflected in the substantial overlap of top-ranked species across components. In the first dietary phase–related component (comp. 1, both models; Fig 4 ), 87% of the top 15 species with the largest absolute weights overlapped between models. Overlap was complete (100%) for the recovery phase–related components ( S11 Fig panels i, iv), with species typically associated with exposure to antibiotics such as Veillonella and Klebsiella being identified. The second dietary phase–related component showed greater divergence, with 60% overlap between the two models ( S12 Fig ). Here, the microbial signature in both models indicated an increased relative abundance of Clostridium species, while CP attributed larger absolute weights to the simultaneously decreased relative abundance of species, such as Eubacterium eligens and Bifidobacterium adolescentis ( S12 Fig , panels i, iv). These patterns reflected the difference in microbial composition between the EEN group compared to the other diet groups which were consistent between PARAFAC2 and CP ( S12 Fig , panels iii, vi).

We applied PARAFAC2 and CP to shotgun metagenomics-based longitudinal gut microbiome data from adults in the FARMM study to explore the effects of dietary and antibiotic/polyethylene glycol intervention (Abx/PEG) [ 16 ]. The replicability-based model selection approach resulted in a 3-component PARAFAC2 and a 3-component CP model ( S8 Fig ). Due to the smaller sample size (N = 10 in each of 3 study groups), the replicability analysis was conducted with F = 5 subsets, with stratified sampling for study group. The PARAFAC2 model explained a larger proportion of the variance in the data (51.0%) compared to CP (48.3%). Both models identified two components related to the dietary phase of the study (days 1–5) and one component related to the recovery phase (days 9–15) following the Abx/PEG perturbation ( S9 Fig and S10 Fig ).

The second early-life component (PARAFAC2 comp. 5 and CP comp. 7) captured a microbial signature with higher abundances of Bacteroides species in both models, accompanied in CP by a simultaneously decreased abundance of taxa such as Clostridium perfringens ( Fig 3 , panel ii; S7 Fig , panel ii). The subject scores for these signatures differed significantly by mode of delivery in both models ( Fig 3 , panel vi; S7 Fig , panel vi), consistent with known effects of Cesarean section [ 45 , 46 ]. The CP model suggests that this signature diminishes over time in all groups ( S7 Fig , panel iv), implying a uniform replacement of early colonizers [ 47 , 48 ]. In contrast, the PARAFAC2 scaled time loadings show that the signature remains relatively stable in CS–born infants (characterized by consistently lower initial abundance), whereas in vaginally delivered infants, the initially high abundance decreases over time ( Fig 3 , panel iv). This suggests that the uniform trajectory in CP may be an artifact of the model’s rigidity. Crucially, by examining subject-specific rather than group-mean profiles, PARAFAC2 unmasked substantial heterogeneity in gut microbiome maturation ( Fig 3 , panel vii; vs S7 Fig , panel vii). Unlike CP, which forces a shared temporal profile, PARAFAC2 successfully captured individual temporal shifts and delays, providing a more accurate representation of the subject-specific pace of microbiome development.

Both PARAFAC2 and CP recovered a component representative of a general early-life microbial signature (PARAFAC2 comp. 3 and CP comp. 2), characterized by typical early colonizers such as Escherichia species and members of the Enterobacteriaceae family ( Fig 3 , panel i; S7 Fig , panel i) [ 44 ]. In both models, the scaled time loadings show that the prevalence of this signature decreases gradually over time ( Fig 3 , panel ii; S7 Fig , panel ii). However, the added flexibility of PARAFAC2 revealed more differential dynamics across exposure groups that CP failed to resolve. While CP subject loadings indicated a difference between VD-no IAP and VD-IAP ( S7 Fig , panel v), PARAFAC2 distinguished VD (both without and with IAP) from CS ( Fig 3 , panel v). Furthermore, the PARAFAC2 scaled time loadings captured a temporal interaction: the microbial signature was more dominant in CS compared to VD specifically at 1 week ( Fig 3 , panel ii), a detail obscured in the CP model.

The early-life representations of the PARAFAC2 and CP models were broadly consistent in terms of the uncovered microbial signatures ( Fig 3 , panels i, iii; S7 Fig , panels i, iii). The top 15 taxa with the largest absolute weights in each of the early-life-related components (PARAFAC2 comp. 3 vs. CP comp. 2, and PARAFAC2 comp. 5 vs. CP comp. 7) showed an 87% and 80% overlap between the models at the ASV level. When considering the union of top 15 ASVs across the components, the overlap between models is 89%, highlighting the consistency in the extracted microbial signatures. However, while the taxa profiles were similar, the models differed in their ability to resolve group-level and subject-specific temporal dynamics ( Fig 3 , panels ii, iv; S7 Fig , panels ii, iv).

Components 3 (i–ii) and 5 (iii–iv) capture patterns in gut microbiome composition during the first year of life. Panels (i, iii) show taxa loadings with the 15 ASVs of largest absolute weight in the component, colored by genus. Scaled time loadings (ii, iv) at time points were averaged per exposure group. Shaded area indicate standard errors of group means. Panels (v, vi) present subject loadings in components 3 and 5 compared across exposure groups using the Mann–Whitney U test with Bonferroni correction. Panel (vii) illustrates subject-specific scaled time loadings of component 5, highlighting heterogeneity in the trajectories.

We selected a 5-component PARAFAC2, and a 7-component CP model based on our model selection approach, focusing on replicability and interpretation of the components ( S4 Fig ). The replicability analysis was conducted using F = 10 subsets, generated through sampling stratified by delivery mode to account for its known impact on gut microbiome development [ 42 , 43 ]. For PARAFAC2, the dimensions of the dataset limited the maximum number of extractable components to five, in line with the model’s uniqueness conditions [ 25 ]. Notably, PARAFAC2 explained slightly more variance (50.0%) with fewer components compared to CP (49.8%), reflecting its ability to capture temporal misalignment more efficiently. Both models identified components that capture early-life patterns—specifically, variation in samples collected at 1 week, 1 month, and 1 year of age (PARAFAC2: comp. 3 & 5, Fig 3 ; CP: comp. 2 & 7, S7 Fig ). The complete models, including all components, are visualized in S5 Fig and S6 Fig .

To characterize infant gut microbiome maturation while accounting for subject-specific differences in the temporal evolution of dominant microbial signatures, we applied the PARAFAC2 decomposition to 16S rRNA longitudinal gut microbiome data from 267 children of the COPSAC 2010 cohort. To highlight the distinctions in the extracted patterns, we also compared the results with those obtained using the CP decomposition. We focus on the model components capturing microbial variation in early life and related to mode of delivery and maternal intrapartum antibiotic exposure to illustrate the approaches. For downstream analysis and visualization, subjects were grouped into exposure groups according to mode of delivery and maternal intrapartum antibiotics: (i) vaginal delivery without intrapartum antibiotic prophylaxis (VD-no IAP), (ii) vaginal delivery with intrapartum antibiotic prophylaxis (VD-IAP), and (iii) Cesarean section (with intrapartum antibiotic prophylaxis; CS).

PARAFAC2, however, successfully captured the between-subject variability in comp. 2. In particular, it faithfully recovered the temporal patterns of subjects 4–9 ( Fig 2 , bottom panel, highlighted). While CP forced a shared temporal profile, PARAFAC2 accurately resolved the saturation timing differences, achieving a mean absolute error in of 0.0013 versus 0.0068 for CP. This demonstrates the superior performance of PARAFAC2 in capturing subject-specific dynamics. Results were similar, with PARAFAC2 outperforming CP at 25% noise ( S3 Fig ).

While both models achieved high overall fit (CP: 98.8%, PARAFAC2: 99.9%) as well as factor match scores (CP: 0.994, PARAFAC2: 0.997) and successfully recovered the common taxa factors ( ) and temporal dynamics in components 1 and 3, they differed substantially in handling subject-specific variation. CP did not recover the subject-specific time loadings of comp. 2 due to its assumption of shared time loadings between subjects ( Fig 2 , middle panel, highlighted with solid borders); the recovered profiles misrepresented the trajectories as identical across subjects, with only scale differences due to subject loadings .

To illustrate the utility of PARAFAC2 in microbiome data analysis on a controlled example with known ground truth, we compared its performance with that of CP in recovering the underlying true factors from a simulated dataset specifically designed to reflect subject-specific variation in the temporal patterns. The ground truth as well as the recovered factors are shown in Fig 2 .

4. Discussion

In this study, we introduced the PARAFAC2 model as a powerful tool for the unsupervised analysis of longitudinal microbiome data, specifically to capture subject-specific temporal trajectories. Furthermore, to ensure the robustness of the uncovered patterns, we proposed a replicability-based framework to determine the number of components, ensuring the biological relevance and stability of the extracted signatures.

Results from the simulated data expose a critical limitation of CP-based tensor decompositions for longitudinal microbiome data analysis: the inability to capture subject-specific microbial dynamics (such as accelerations or delays), phenomena frequently observed in real populations [50,51]. This limitation was explicitly demonstrated by the peak timing of the microbial signature in component 2. While the ground truth exhibited a broad temporal shift, peaking early at for subject 4 and late at for subject 9, PARAFAC2 precisely resolved these distinct timings. In contrast, the CP model masked this heterogeneity, rigidly forcing a shared temporal profile that peaked at an intermediate for all subjects. Consequently, CP compensated for this timing misalignment by introducing variability in peak magnitude not supported by the data. In real-world contexts, where shifts in saturation timing reflect distinct environmental exposures or host states, accurate recovery of these individual trajectories is essential.

While PARAFAC2 outperformed CP in terms of fit and FMS on this illustrative example, it is important to note that perfect scores are not expected, as the simulated data did not fully conform to the structural assumptions of either model. Beyond this, the simulation was not designed to reproduce the full complexity of real microbial community dynamics, where the true data-generating process is generally unknown. In such settings, the added flexibility of PARAFAC2 may make it the more suitable modeling choice when individual variation in temporal profiles is expected or of biological interest.

In real data, similarly to the simulation, PARAFAC2 outperformed CP in terms of fit, and the two models revealed broadly consistent microbial signatures. However, the scaled time loadings of PARAFAC2 provided a more fine-scale temporal evolution of these signatures at the individual level, as well as additional group differences not found by CP. The subject-specific PARAFAC2 patterns revealed heterogeneity among Cesarean-section-born infants, highlighting that some show an earlier change toward a composition resembling that of vaginally born infants. Identifying these subject-specific maturation trajectories is important, as they may reflect differences in postnatal exposures, including environmental conditions or feeding practices that help to compensate for the early-life perturbations [49,52]. In the FARMM study, PARAFAC2 resolved a compositional change in the EEN group, characterized by decreasing relative abundances of species such as Faecalibacterium prausnitzii, Parabacteroides merdae, and Bacteroides ovatus during the dietary phase—dynamics consistent with [16], which were obscured by the CP model. The assumptions of CP require that all subjects follow the same temporal trajectory, allowing only for differences in the magnitude of how strongly each individual expresses the microbial signature according to the trajectory. Because of this assumption, inter-individual differences in the timing or progression of microbial changes could not be captured using CP, hindering the model’s ability to represent subject-specific dynamics and overlooking biologically meaningful patterns. Nevertheless, the increased complexity of PARAFAC2 may not always be necessary. For datasets or research questions where the primary goal is to characterize dominant, shared temporal trends rather than individual heterogeneity, the simpler CP model remains a robust and effective choice.

A limitation of PARAFAC2 lies in the constant cross-product constraint it imposes on the factor matrices. Although the constraint is essential for obtaining a unique solution, which is required for interpretability, the assumption that loading vectors of should have constant angles for all k may impose biologically unwarranted structure [53], potentially limiting the ability of the model to fully capture the diversity of subject-specific temporal dynamics.

In practice, the constraint implies that the different microbial signatures’ temporal profiles should be similar (more specifically, have constant angles) across subjects (see S1 Fig for a visual illustration). This formulation leads to more expressive time patterns with PARAFAC2, admitting shifts and delays in subjects’ microbial signatures as opposed to the restrictive assumption of CP, which forces the temporal patterns to be identical across subjects. If subjects have structurally different temporal profiles, i.e., when the temporal profiles of the microbial signatures are similar (based on the above criteria) in some subjects but dissimilar in others, the constraint may be unreasonable. In such cases, however, PARAFAC2 may still be more realistic than CP. Additionally, the number of components PARAFAC2 can extract is limited by its uniqueness conditions. In particular, for datasets with only a few time points, a CP-based approach may be a better choice, as it may recover more microbial signatures. More generally, the choice between CP and PARAFAC2 should be guided by the research question, as well as whether the subject-specific patterns extracted by PARAFAC2 add biologically meaningful structure beyond CP; if they do not, the simpler CP model may be preferable. Notably, PARAFAC2 has emerged as a robust compromise between models of differing flexibility in the analysis of functional magnetic resonance imaging (fMRI), effectively capturing connectivity and spatial components while accounting for subject variability [54]. In addition, the model has been shown to perform well even when the PARAFAC2 constraint is only approximately satisfied, as demonstrated in our results on simulated data and in previous studies [32,25].

An open challenge when fitting PARAFAC2 is the presence of completely missing measurements at a time point for a subject (i.e., a missing mode-1 fiber, such as in the case of the FARMM dataset), which makes the corresponding entries of the time loadings arbitrary as long as the imposed constraints are satisfied. To overcome this indeterminacy and improve the robustness of our results, we employed a smoothness constraint on the temporal mode. Our ablation analysis introducing additional mode-1 fiber missingness in the FARMM cohort (S15 Fig, S16 Fig) supports the robustness of this approach at missingness rates above those encountered in our real data examples. However, at the individual level, in cases where such missingness is extensive or affects consecutive missing measurements (e.g., 9032, panel iv, Fig 5), this approach may be inadequate, highlighting the need for alternative strategies for the handling of missing mode-1 fibers.

Validating the subject-specific patterns of PARAFAC2 is challenging, owing to its added flexibility. In the absence of biological replicates, it is not possible to directly assess whether the subject-specific trajectories reflect true biological variation or are partly driven by model assumptions and data idiosyncrasies. To account for this, we focused on group-level trends, aggregating subject-specific loadings within predefined exposure groups (e.g., by delivery mode and diet group). While this approach captures consistent between-group differences, it may mask individual-level heterogeneity. Therefore, to more directly assess replicability at the subject level, we inspected the heterogeneity of subject-specific scaled time loadings between the submodels fitted during replicability analysis. By leveraging the partial overlap of subjects between submodels, we could compare subject-specific trajectories across submodels, which confirmed their replicability and supported our model selection approach. Nevertheless, parameters of the resampling, such as the number of subsamples, depend on dataset characteristics. Finally, replicate samples, when available, can aid in both determining the number of model components and assessing the biological validity of the discovered patterns. This approach has been demonstrated in metatranscriptomic analyses of microbial communities using CP [55], and the study of functional neuroimaging data via PARAFAC2 [20].

Beyond the unsupervised formulation considered here, an interesting direction for future work is to incorporate group structure directly into the model. A group-aware variant of PARAFAC2, in which the subject-specific matrices are shared within predefined groups or shrunk toward group-level centroids, would provide a middle ground between the shared dynamics of CP and the fully subject-specific formulation used here. Related ideas have been explored in the matrix factorization setting, most notably by MEFISTO [56]; an analogous tensor-decomposition formulation that preserves subject-level temporal flexibility remains an open methodological direction.

In summary, PARAFAC2 presents a flexible and interpretable approach to analyze longitudinal microbiome data while explicitly modeling subject-specific temporal trajectories. Accommodating timing differences, shifts and delays, PARAFAC2 enables a more accurate characterization of temporal changes. The integration of replicability as a criterion for model selection enhances the robustness of the extracted patterns. Together, these features advance the methodological toolkit for studying microbiome dynamics and open avenues for future applications across diverse longitudinal omics datasets.

[END]
---
[1] Url: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1014486

Published and (C) by PLOS One
Content appears here under this condition or license: Creative Commons - Attribution BY 4.0.

via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/