MARSH: Drift-Robust Microbial Source Tracking via Tree-Wasserstein Optimal Transport

1. Introduction

Microbial community assembly reflects interacting processes that include dispersal, selection, diversification, and ecological drift [1]. An observed community may also contain microbial contributions from several distinct sources. Understanding these contributions is important in settings where the origin of a community carries biological, environmental, or public-health significance. Microbial source tracking (MST) provides a framework for this task by estimating the relative contributions of candidate source communities to a target, or “sink,” community [2,3]. MST has been used to identify sources of fecal contamination in recreational and environmental waters, where distinguishing human, animal, wastewater, and other inputs can inform pollution control and water-quality management [4,5]. It has also been applied to trace microbial contamination and dispersal in hospitals, laboratories, offices, and other built environments [2,6], and to study the formation of host-associated microbiomes, including the contributions of maternal and early-life communities to the developing infant gut [3]. More broadly, source attribution provides a way to move beyond describing which microbes are present and ask where an observed community came from.

In a typical MST analysis, each source and sink is represented by its relative-abundance profile, and the goal is to estimate the fraction of the sink attributable to each candidate source [2,3]. These profiles are compositional because sequencing data describe the abundance of each taxon relative to the total observed community rather than providing independent absolute abundances [7]. Source attribution is further complicated by the fact that the reference library may not contain every community that contributed to the sink. An environmental or host-associated sample, for example, may contain input from an unsampled reservoir or from a source that was not included when the reference panel was constructed. MST methods must therefore estimate the contributions of the available sources while also allowing part of the sink to originate from an unknown source [2,3].

Beyond incomplete source libraries, MST can also be affected when the microbial composition of the sink no longer aligns exactly with that of the reference sources. We refer to this source–sink mismatch as drift. Drift can arise biologically as microbial communities change over time [8], or technically when source and sink samples are characterized using different 16S regions, primers, or sequencing protocols [9]. In either case, taxa observed in the sink may differ from those in the references even when they are phylogenetically related. Existing MST approaches based on Bayesian mixture models, expectation–maximization, and matrix factorization operate on a shared feature space and match abundance across corresponding taxa [2,3,10]. They therefore cannot directly use the relationship between a taxon present in the sink and a closely related taxon present in a source. As drift increases, this mismatch can lead source contributions to be misestimated or attributed to the unknown component.

Phylogeny provides a natural way to incorporate this relationship. Phylogenetic community distances were developed to distinguish changes among close relatives from changes across deeper lineages [11,12]. Optimal transport makes this distinction explicit by treating microbial communities as probability measures on a phylogenetic tree. In this setting, the tree-based Wasserstein-1, or Kantorovich–Rubinstein, distance coincides with raw weighted UniFrac and has a closed-form edge-sum representation on the tree [12,13]. Transport between nearby taxa therefore costs less than transport between taxa separated by long phylogenetic distances. MARSH uses this geometry directly in the source-mixture estimation problem, allowing phylogenetically local differences between the sink and its references to be treated as drift rather than automatically as evidence for a different or unknown source. Our contributions are threefold:

1.

We formulate MST as a constrained optimal-transport problem on the phylogenetic tree, estimating source mixing weights by minimizing the tree-Wasserstein distance between the source mixture and the sink. The resulting problem is solved as an exact least-absolute-deviations linear program (Section 3.4).

2.

We introduce a fixed locality rule that limits the loss to shallow phylogenetic structure. This preserves the ability to match locally drifted taxa while preventing broad clade-level imbalances from dominating the fit (Section 3.3).

3.

We introduce an estimator of the unknown-source fraction based on the residual mismatch after fitting the known sources. It is designed to distinguish unexplained source contribution from mismatch that can be explained by local phylogenetic drift (Section 3.5).

We evaluate these components against FEAST, SourceID-NMF, and FastST on simulated and real 16S data, including unknown sources, induced phylogenetic drift, and naturally occurring temporal drift.

2. Related Work

Most quantitative microbial source tracking methods begin from the standard mixture assumption. Let $\mathbf{x}\in {\mathbb{R}}_{+}^n$ denote the observed sink, let ${\mathbf{s}}_1,\dots ,{\mathbf{s}}_K$ denote candidate source profiles over the same n taxa, and let $\mathbf{w}=(w_1,\dots ,w_K)$ denote their mixing proportions. In its simplest form, MST assumes that the expected sink composition can be written as

$$\mathbf{p}\left(\mathbf{w}\right)={\displaystyle \sum _{k=1}^K}{w}_k{\mathbf{s}}_k,w_k\ge 0,{\displaystyle \sum _{k=1}^K}{w}_k=1.$$

(1)

Methods differ primarily in how the source profiles are modeled, how an unobserved source is represented, and what statistical discrepancy is used to compare the mixture with the sink. The methods considered below all use a common taxonomic feature space: taxon j in the sink is modeled through abundance associated with taxon j in the sources. Some methods couple residuals across taxa statistically, but none of these formulations assigns a smaller mismatch cost merely because two different taxa are phylogenetically close. We summarize the main mathematical approaches to MST below.

One influential approach is the Bayesian mixture model used by SourceTracker [2], which formulates MST as a hierarchical Bayesian mixture model at the sequence level. Let ${\theta }_k\in {\Delta }^{n-1}$ denote the taxon distribution of source environment k, and let $\mathbf{w}$ denote the proportions of the source environments in the sink. For sink read r, the model introduces a latent source assignment $z_r$:

$$\begin{array}{cc}\hfill \mathbf{w}& \sim \mathrm{Dirichlet}\left(\beta \right),\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill {\theta }_k& \sim \mathrm{Dirichlet}\left({\alpha }_k\right),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill z_r\mid \mathbf{w}& \sim \mathrm{Categorical}\left(\mathbf{w}\right),\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill x_r\mid z_r=k,{\theta }_k& \sim \mathrm{Categorical}\left({\theta }_k\right).\hfill \end{array}$$

(5)

Reference-source observations update the corresponding ${\theta }_k$, while Gibbs sampling integrates uncertainty in both source-specific taxon frequencies and sink mixing proportions. An additional source category represents taxa attributable to an unobserved environment.

The categorical likelihood has taxon-specific support. If sink read r is observed as taxon j, then under source k its likelihood is

$$p(x_r=j\mid z_r=k)={\theta }_{kj}.$$

(6)

Thus, the model can use the probability assigned to taxon j by a source, but it contains no term involving the relationship between j and another taxon $j^{\prime }$. A source with abundance at a close phylogenetic relative receives no special credit unless that relationship is already reflected in the shared taxon frequencies themselves.

FEAST [3] instead uses a maximum-likelihood mixture model. It retains the probabilistic mixture formulation but replaces Markov-chain Monte Carlo inference with expectation–maximization. Let $\mathbf{x}=(x_1,\dots ,x_n)$ be sink counts with total count $C={\sum }_j{x}_j$, and let ${\mathbf{y}}_k=(y_{k1},\dots ,y_{kn})$ be the observed counts for known source k with total $C_k$. FEAST introduces source-specific relative-abundance vectors ${\gamma }_k$ and mixing proportions $\alpha$, including an additional unobserved source. The sink probabilities are

$${\beta }_j={\displaystyle \sum _{k=1}^{K+1}}{\alpha }_k{\gamma }_{kj},{\displaystyle \sum _{k=1}^{K+1}}{\alpha }_k=1,{\displaystyle \sum _{j=1}^n}{\gamma }_{kj}=1.$$

(7)

The generative model is

$${\mathbf{y}}_k\sim \mathrm{Multinomial}\left(C_k,{\gamma }_k\right),\mathbf{x}\sim \mathrm{Multinomial}\left(C,\beta \right).$$

(8)

Up to terms independent of the unknown parameters, the joint log-likelihood is

$$\ell (\alpha ,\gamma )={\displaystyle \sum _{j=1}^n}{x}_{j}log\left({\displaystyle \sum _{k=1}^{K+1}}{\alpha }_k{\gamma }_{kj}\right)+{\displaystyle \sum _{k=1}^K}{\displaystyle \sum _{j=1}^n}{y}_{kj}log{\gamma }_{kj}.$$

(9)

The first term fits the source mixture to the sink, while the second anchors the known latent source profiles to their observed reference counts. For fixed source profiles, maximizing the sink term is equivalent, up to an additive constant and a factor of C, to minimizing

$$D_{\mathrm{KL}}\left(\widehat{\mathbf{x}}\parallel \beta \right),{\widehat{x}}_j={\displaystyle \frac{x_j}{C}},$$

(10)

where $D_{\mathrm{KL}}$ is the Kullback–Leibler divergence.

The taxon alignment becomes particularly clear in the EM responsibility for assigning abundance of taxon j to source k:

$$r_{kj}=Pr(z=k\mid j)={\displaystyle \frac{{\alpha }_k{\gamma }_{kj}}{{\sum }_{\ell =1}^{K+1}{\alpha }_{\ell }{\gamma }_{\ell j}}}.$$

(11)

Only the abundance of the same taxon j enters this assignment. Abundance at a phylogenetically nearby taxon $j^{\prime }$ cannot directly explain the observation at j.

SourceID-NMF [10] uses non-negative matrix factorization to relax the assumption that the taxonomic profile of a source in the sink must be identical to its observed reference profile. Let $X\in {\mathbb{R}}_{+}^{n\times 1}$ denote the sink, let $Y\in {\mathbb{R}}_{+}^{n\times (K+1)}$ contain the reference source profiles, and let $W\in {\mathbb{R}}_{+}^{n\times (K+1)}$ and $H\in {\mathbb{R}}_{+}^{(K+1)\times 1}$ denote the latent source profiles and their proportions, respectively. The final column corresponds to an unknown source and is masked from the reference penalty by A. SourceID-NMF solves

$$\begin{array}{cc}\hfill \underset{W,H\ge 0}{min}& \frac{1}{2}{\parallel X-WH\parallel }_F^2+\frac{1}{2}{\parallel A\circ (W-Y)\parallel }_F^2\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \sum _k{H}_k=1,\sum _i{W}_{ik}=1\forall k.\hfill \end{array}$$

(12)

where ∘ denotes the Hadamard product. The problem is optimized using the alternating direction method of multipliers.

This formulation is more flexible than fixing the source profiles to Y: the latent profile $W_k$ can deviate from the observed reference $Y_k$. However, both terms remain Euclidean in the taxon coordinates. In particular,

$${\parallel X-WH\parallel }_F^2={\displaystyle \sum _{i=1}^n}{\left(X_i-{\displaystyle \sum _{k=1}^{K+1}}{W}_{ik}{H}_k\right)}^2$$

(13)

and

$${\parallel A\circ (W-Y)\parallel }_F^2={\displaystyle \sum _{i=1}^n}{\displaystyle \sum _{k=1}^{K+1}}{A}_{ik}^2{(W_{ik}-Y_{ik})}^2.$$

(14)

Changing a source profile from taxon i to taxon j is therefore not assigned a cost based on the biological distance between i and j. A transfer between two closely related taxa and a transfer between distant clades are distinguished only through the resulting numerical changes in the abundance vectors.

FastST [14] takes a generalized least-squares approach, starting from the same multinomial mixture relation

$${\beta }_j={\displaystyle \sum _{k=0}^K}{\alpha }_k{\gamma }_{kj},1\le j\le n,$$

(15)

where $k=0$ denotes an unobserved source. Instead of jointly estimating the source taxon distributions and mixing proportions iteratively, it treats the known ${\gamma }_{kj}$ as given, or estimates them from the observed source frequencies. Multiplying (15) by the sink sequencing depth C gives

$$\mathbb{E}\left[X_j\right]={\displaystyle \sum _{k=0}^K}C{\alpha }_k{\gamma }_{kj}.$$

(16)

Because both the abundance profile and proportion of the unobserved source are unidentified, FastST assumes that ${\alpha }_0{\gamma }_{0j}$ is constant across taxa. Defining the resulting unknown-source term as an intercept and setting ${\tilde{\alpha }}_k=C{\alpha }_k$ for the observed sources gives the linear model

$$X_j={\tilde{\alpha }}_0+{\displaystyle \sum _{k=1}^K}{\tilde{\alpha }}_k{\gamma }_{kj}+{\epsilon }_j.$$

(17)

Under the multinomial model,

$$Cov\left(\epsilon \right)=C\left[diag\left(\beta \right)-\beta {\beta }^{\top }\right].$$

(18)

FastST estimates this covariance from the sink and uses generalized least squares. Writing $\Gamma =[\mathbf{1},{\gamma }_1,\dots ,{\gamma }_K]$ and the estimated covariance as $\Sigma$, the estimator takes the form

$$\widehat{\tilde{\alpha }}={\left({\Gamma }^{\top }{\Sigma }^{+}\Gamma \right)}^{-1}{\Gamma }^{\top }{\Sigma }^{+}\mathbf{x},$$

(19)

where ${\Sigma }^{+}$ denotes an inverse after removing the multinomial singularity or, equivalently, a Moore–Penrose pseudoinverse.

Unlike a simple coordinate-wise squared-error loss, GLS accounts for the negative covariance among multinomial taxon counts through $\Sigma$. Nevertheless, the design matrix remains taxon-aligned: ${\gamma }_{kj}$ predicts sink abundance at taxon j. The covariance structure models statistical dependence among residuals but does not encode phylogenetic distance between taxon j and another taxon $j^{\prime }$.

Taken together, existing MST methods provide distinct probabilistic and optimization-based formulations for estimating source contributions, including Bayesian inference, maximum-likelihood estimation, matrix factorization, and generalized least squares. Although these approaches differ in their treatment of uncertainty, unknown sources, and computation, they represent sources and sinks on a shared taxonomic feature space and do not explicitly account for phylogenetic distance between nonmatching taxa. Conversely, phylogenetic distances such as weighted UniFrac incorporate evolutionary relatedness when comparing microbial communities [11,12,13], but do not by themselves provide a source-mixture estimator. This separation between source attribution and phylogenetic community comparison motivates a formulation that brings the two together.

3. Methods

MARSH estimates the contributions of candidate microbial sources to an observed sink by fitting a mixture model under a phylogeny-aware transport loss. We first define the source-tracking problem and its unknown-source component, then introduce the phylogenetic discrepancy used by MARSH, its locality modification, and the resulting estimators.

3.1. Problem Formulation

Let the observed sink community contain counts $\mathbf{y}\in {\mathbb{Z}}_{\ge 0}^n$ over n taxa, with total sequencing depth $N={\sum }_i{y}_i$. We work with the corresponding relative-abundance profile

$$\mathbf{p}={\displaystyle \frac{\mathbf{y}}{N}}\in {\Delta }^{n-1},{\Delta }^{n-1}=\mathbf{x}\in {\mathbb{R}}_{\ge 0}^n:\sum _i{x}_i=1.$$

(20)

The reference panel contains K candidate source communities with relative-abundance profiles ${\mathbf{s}}_1,\dots ,{\mathbf{s}}_K\in {\Delta }^{n-1}$. The basic MST task is to estimate the proportion of the sink contributed by each of these sources.

Because the reference panel may not contain every community that contributed to the sink, we also allow an unknown source with profile ${\mathbf{b}}_0\in {\Delta }^{n-1}$ and contribution $w_0$. We model the sink as

$$\mathbf{p}\approx (1-w_0){\displaystyle \sum _{k=1}^K}{w}_k{\mathbf{s}}_k+w_0{\mathbf{b}}_0,\mathbf{w}\in {\Delta }^{K-1},0\le w_0\le 1,$$

(21)

where $\mathbf{w}=(w_1,\dots ,w_K)$ gives the relative proportions of the named sources within the explained component of the sink. The final contribution of named source k is therefore $(1-w_0)w_k$.

The inputs to MARSH are the sink profile $\mathbf{p}$, the candidate source profiles ${\left\{{\mathbf{s}}_k\right\}}_{k=1}^K$, and, when available, a phylogenetic tree T relating the observed taxa. The outputs are the estimated named-source contributions and, when unknown-source estimation is enabled, the unknown fraction ${\widehat{w}}_0$ and estimated profile ${\widehat{\mathbf{b}}}_0$.

The sink and its reference sources need not agree exactly at the taxon level. Under phylogenetically local drift, abundance associated with a source taxon may instead appear at a nearby relative in the sink. A discrepancy based only on matching corresponding taxa treats such a change as an ordinary mismatch and ignores the relationship between the two taxa. MARSH therefore fits the mixture in Eq. (21) using a phylogenetic transport loss, so that mismatches between closely related taxa incur a smaller cost than mismatches between distantly related taxa.

3.2. Phylogenetic Transport Loss

Let T be a rooted phylogenetic tree whose leaves correspond to the n taxa, with edge set E and edge lengths ${\ell }_e>0$. For edge e, let $L\left(e\right)$ denote the set of descendant taxa and define the cumulative mass below that edge as

$$m_e\left(\mathbf{q}\right)=\sum _{i\in L\left(e\right)}{q}_i.$$

(22)

Optimal transport compares two distributions by asking how much mass must be moved, and how far, to transform one distribution into the other [15,16]. When the transport cost between taxa i and j is their patristic distance $d_T(i,j)$ on the phylogeny, the Wasserstein-1 distance is

$$W_{1,T}(\mathbf{p},\mathbf{q})=\underset{\Pi \in \mathcal{U}(\mathbf{p},\mathbf{q})}{min}{\displaystyle \sum _{i=1}^n}{\displaystyle \sum _{j=1}^n}{d}_T(i,j){\Pi }_{ij},$$

(23)

where $\mathcal{U}(\mathbf{p},\mathbf{q})$ is the set of non-negative transport plans with marginals $\mathbf{p}$ and $\mathbf{q}$.

On a tree, this optimization has the closed form

$$W_{1,T}(\mathbf{p},\mathbf{q})=\sum _{e\in E}{\ell }_e\left|m_e\left(\mathbf{p}\right)-m_e\left(\mathbf{q}\right)\right|.$$

(24)

This quantity is the phylogenetic Kantorovich–Rubinstein distance and is equivalent to raw weighted UniFrac for normalized communities [12,13]. Equation (24) provides the key property needed for source tracking: a mismatch between nearby taxa contributes only the length of the short path separating them, whereas a mismatch between distant taxa incurs a larger cost.

The tree representation is also computationally useful. Since each $m_e\left(\mathbf{q}\right)$ is linear in $\mathbf{q}$, the tree-Wasserstein distance is a weighted $L_1$ norm of a linear transformation of $\mathbf{p}-\mathbf{q}$. This allows the source-mixture problem to be written directly as a linear program.

3.3. Local Phylogenetic Transport

Full tree-Wasserstein distance incorporates structure at every phylogenetic scale. This is desirable when source–sink differences occur among nearby taxa, but deep edges can also give substantial weight to differences between large clades. For source attribution, we therefore restrict the phylogenetic contribution of the loss to local tree structure.

Let

$$h=round\left(\sqrt{n}\right)$$

(25)

and retain only edges whose descendant subtree contains at most h taxa. Define

$${\tilde{\ell }}_e={\ell }_e\mathbf{1}\left\{\left|L\right(e\left)\right|\le h\right\}.$$

(26)

The resulting local tree loss is

$$D_T(\mathbf{p},\mathbf{q})=\sum _{e\in E}{\tilde{\ell }}_e\left|m_e\left(\mathbf{p}\right)-m_e\left(\mathbf{q}\right)\right|.$$

(27)

This modification retains all leaf edges and shallow internal edges. Exact taxonomic mismatches therefore remain penalized, while shallow internal edges provide partial tolerance when abundance shifts among nearby phylogenetic relatives. Deep edges are excluded so that broad clade-level imbalances do not dominate the source-mixture fit. The cutoff depends only on the number of taxa represented by the tree and introduces no drift-dependent tuning parameter.

Unless stated otherwise, a phylogenetic tree is supplied to MARSH, which uses the local loss in Eq. (27). The unrestricted tree-Wasserstein loss in Eq. (24) is used as a comparison in our experiments.

3.4. Estimating the Named-Source Proportions

For each edge e and source k, define

$$M_{ek}=m_e\left({\mathbf{s}}_k\right),c_e=m_e\left(\mathbf{p}\right).$$

For a candidate mixture ${\sum }_k{w}_k{\mathbf{s}}_k$, linearity gives

$$m_e\left(\sum _k{w}_k{\mathbf{s}}_k\right)=\sum _k{w}_k{M}_{ek}.$$

The named-source weights are therefore estimated by

$$\begin{array}{cc}\hfill \widehat{\mathbf{w}}& =\underset{\mathbf{w}\in {\Delta }^{K-1}}{arg\; min}\sum _{e\in E}{\tilde{\ell }}_e\left|{\displaystyle \sum _{k=1}^K}{w}_k{M}_{ek}-c_e\right|.\hfill \end{array}$$

(28)

Introducing one auxiliary variable $t_e$ for each retained edge gives the least-absolute-deviations linear program

$$\begin{array}{cc}\hfill \underset{\mathbf{w},\mathbf{t}}{min}& \sum _{e\in E}{\tilde{\ell }}_e{t}_e,\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& -t_e\le \sum _k{w}_k{M}_{ek}-c_e\le t_e,e\in E,\hfill \\ \hfill & \sum _k{w}_k=1,w_k\ge 0.\hfill \end{array}$$

(29)

The problem contains $K+\left|E\right|$ variables and a number of constraints linear in the number of tree edges, and can therefore be solved directly with a standard linear-programming solver.

When no phylogenetic tree is supplied, the taxa are attached to a flat star scaffold. For the fully observed estimator, this makes the tree-Wasserstein objective equal to coordinate-wise $L_1$ distance. When the unknown-source estimator is enabled without a phylogenetic tree, its weight update instead uses simplex-constrained least squares. Define ${\tilde{\mathbf{s}}}_0={\mathbf{b}}_0$ and ${\tilde{\mathbf{s}}}_k={\mathbf{s}}_k$ for $k=1,\dots ,K$. Conditional on the current unknown profile, the update is

$$\widehat{\tilde{\mathbf{w}}}=\underset{\tilde{\mathbf{w}}\in {\Delta }^K}{arg\; min}{∥{\displaystyle \sum _{k=0}^K}{\tilde{w}}_k{\tilde{\mathbf{s}}}_k-\mathbf{p}∥}_2^2.$$

(30)

This weight update alternates with a soft responsibility update of ${\mathbf{b}}_0$ until the weights converge. The residual-based procedure below is used when a phylogenetic tree is supplied. An optional reweighted sparsity penalty may be added when the candidate source library is large; all default experiments use $\gamma =0$ unless otherwise specified.

3.5. Residual-Based Unknown-Source Estimation with a Phylogenetic Tree

Estimating an unknown source introduces an identifiability problem. A discrepancy between the sink and the named-source mixture can arise either because a genuine source is absent from the reference panel or because a known source has changed relative to its reference profile. Treating every unexplained taxon as unknown therefore risks confusing drift with a missing source.

We first fit only the named sources and define

$$={\displaystyle \sum _{k=1}^K}{\widehat{w}}_k{\mathbf{s}}_k,R=D_T(,\mathbf{p}),$$

(31)

where R is the residual discrepancy that cannot be eliminated by reweighting the named sources.

We construct a candidate unknown-source profile from the positive taxonomic residual,

$${\mathbf{b}}_0=\Pi \left[{(\mathbf{p}-)}_{+}\right],$$

(32)

where $\Pi$ renormalizes a non-zero vector to the simplex. If the residual is zero, we set $w_0=0$ and no unknown profile is required.

The transport loss has a useful scaling property. For any $0\le a\le 1$,

$$D_T\left(,(1-a)+a{\mathbf{b}}_0\right)=a{D}_T(,{\mathbf{b}}_0).$$

(33)

Thus, if the residual were generated entirely by mixing an unknown profile ${\mathbf{b}}_0$ into the named mixture, its fraction would naturally scale as $R/D_T(,{\mathbf{b}}_0)$. Under heavy drift, however, the residual profile may lie very close to the named mixture, making this denominator small and causing the ratio to overestimate the unknown fraction. We therefore introduce the intrinsic source-separation scale

$$\overline{D}={\displaystyle \frac{1}{K}}{\displaystyle \sum _{k=1}^K}{D}_T({\mathbf{s}}_k,)$$

(34)

and estimate

$${\widehat{w}}_0=min{\displaystyle \frac{R}{max\left(D_T({\mathbf{b}}_0,),\overline{D}\right)}},1.$$

(35)

The floor prevents a small residual-to-mixture distance from producing a large unknown estimate and is determined entirely by the observed source geometry rather than by a tuned threshold.

After estimating ${\widehat{w}}_0$ and ${\mathbf{b}}_0$, we remove the estimated unknown contribution,

$${\mathbf{p}}^{\prime }=\Pi \left[{(\mathbf{p}-{\widehat{w}}_0{\mathbf{b}}_0)}_{+}\right],$$

(36)

and refit the named-source weights. In our implementation this update is repeated for three rounds, after which the reported source proportions are

$$\left((1-{\widehat{w}}_0){\widehat{w}}_1,\dots ,(1-{\widehat{w}}_0){\widehat{w}}_K,{\widehat{w}}_0\right).$$

The distinction between drift and an unknown source follows directly from the geometry of the loss. If an amount of mass $a_r$ is displaced from taxon $i_r$ to taxon $j_r$, then

$$D_T(,\mathbf{p})\le \sum _r{a}_r{\tilde{d}}_T(i_r,j_r),$$

(37)

where ${\tilde{d}}_T$ is the path cost induced by the retained tree edges. Drift between nearby taxa therefore produces a small residual, whereas source mass that cannot be reconciled with the named mixture produces a larger residual. This distinction is not absolute: an unknown source that is itself phylogenetically similar to the named sources may also produce a small transport residual. We return to this identifiability limit in Section 6.

3.6. Uncertainty

For fits without an unknown source, we quantify sampling uncertainty with a multinomial bootstrap. Given sink counts $\mathbf{y}$ at sequencing depth N, bootstrap replicate b is drawn as

$${\mathbf{y}}^{\left(b\right)}\sim \mathrm{Multinomial}(N,\mathbf{p}),b=1,\dots ,B.$$

(38)

Each reference source is independently resampled in the same way at its observed sequencing depth. The resampled profiles are renormalized and the estimator is rerun to obtain ${\mathbf{w}}^{\left(b\right)}$. A $100(1-\alpha )\%$ confidence interval for each estimated proportion is given by the empirical percentile interval

$$\left[Q_{\alpha /2},Q_{1-\alpha /2}\right],$$

(39)

where $Q_q$ denotes the empirical q-quantile of the corresponding B bootstrap estimates. The current implementation reports point estimates, but not bootstrap intervals, when the unknown-source estimator is enabled.

4. Experimental Design

We evaluate MARSH on simulated and real 16S data under three settings: fully observed sources, unknown sources, and source–sink drift. The simulated experiments provide controlled ground truth, while the real-data experiments use observed microbial communities and taxonomy.

4.1. Simulated Data

Each simulated replicate is generated using our export-sim tool. We first draw a random birth–death phylogeny with n leaves and positive edge lengths. For $K=5$ sources, the source profiles are sampled as

$${\mathbf{s}}_k\sim \mathrm{Dirichlet}\left(\alpha \mathbf{1}\right),k=1,\dots ,K.$$

(40)

Ground-truth mixing weights ${\mathbf{w}}^{★}$ are then selected. When an unknown source is included, we also generate an unknown profile ${\mathbf{b}}_0$ with fraction $w_0^{★}$. The sink distribution is formed using Eq. (21), and observed counts are sampled as

$$\mathbf{y}\sim \mathrm{Multinomial}(N,\mathbf{p}).$$

(41)

Thus, all methods operate on finite-depth count data rather than the underlying probability vectors.

Three controlled evaluations use this generator. The fully observed scaling evaluation (Section 5.1) uses $n\in \{50,500,5,000,50,000\}$ and $w_0^{★}=0$. The unknown-source evaluation without drift (Section 5.2) fixes $n=500$ and varies the unknown fraction $w_0^{★}\in \{0.2,0.4,0.6,0.8\}$. The simulated drift evaluation (also Section 5.2) fixes $n=200$ and $w_0^{★}=0.30$. Each condition is repeated over five independent replicates.

4.2. GlobalPatterns

We use the GlobalPatterns 16S dataset [17], processed using QIIME 2 and phyloseq [18]. The dataset contains eight environment types: Feces, Skin, Tongue, Freshwater, Freshwater (creek), Ocean, Sediment (estuary), and Soil. Samples from each environment are pooled to form a source profile ${\mathbf{s}}_e$.

After filtering low-prevalence taxa, $6,231$ taxa remain. The accompanying phylogeny is midpoint-rooted. For the drift experiment, we use a $1,000$-taxon subtree to reduce the cost of repeated tree-based optimization. We use a leave-one-environment-out design. For each environment U, we remove U from the candidate source set and construct the sink from two retained environments, A and B, together with the withheld environment:

$$\mathbf{p}=0.40{\mathbf{s}}_A+0.30{\mathbf{s}}_B+0.30{\mathbf{s}}_U.$$

(42)

The true sink contributions are ${\rho }_A^{★}=0.40$, ${\rho }_B^{★}=0.30$, and ${\rho }_0^{★}=0.30$. The environments A and B are selected using a fixed random seed. The remaining candidate set contains seven environments, including five zero-weight distractors. Since U is withheld from the candidate set, its contribution can only be recovered through the unknown-source estimate. Repeating this procedure for all eight environments gives eight paired trials. The no-drift and genus-drift evaluations are reported in Section 5.3.

4.3. Drift Models

We evaluate two forms of drift. The simulated drift evaluation (Section 5.2) uses the TADA generative model [19] to introduce structured phylogenetic perturbations in the sink. For the GlobalPatterns evaluation (Section 5.3), we define taxonomy-based genus turnover. Each eligible reference taxon i is paired with a source-absent taxon $\pi \left(i\right)$ from the same genus. At drift level f, a fraction of the abundance at taxon i is transferred to its congener:

$$\begin{array}{cc}\hfill p_{\pi \left(i\right)}& \leftarrow p_{\pi \left(i\right)}+f{p}_i,\hfill \\ \hfill p_i& \leftarrow (1-f)p_i.\hfill \end{array}$$

(43)

The profile is renormalized after all transfers. At $f=0$, no drift is applied, while at $f=0.5$ half of each affected reference taxon’s mass is transferred to its congener. The mapping $\pi$ is defined using genus annotations rather than distances on the phylogenetic tree used by MARSH. The drift construction is therefore independent of the tree used by the estimator. For both simulated and real data, drift is evaluated at $f\in \{0.1,0.2,0.3,0.4,0.5\}$; we do not evaluate or claim behavior beyond $f=0.5$.

4.4. Baselines and Evaluation

We compare MARSH with the published implementations of FEAST [3], SourceID-NMF [10], and FastST [14]. FEAST is run using its R implementation, SourceID-NMF using its Python implementation with 300 iterations, and FastST using its R implementation.

Accuracy is measured using the $L_1$ error between the estimated and true source proportions:

$$L_1={\displaystyle \sum _{k=0}^K}\left|{\widehat{w}}_k-w_k^{★}\right|,$$

(44)

where $k=0$ denotes the unknown-source fraction. Lower values indicate more accurate source recovery.

Simulated results are reported descriptively as means over five independent replicates. Real-data results use paired observations from the eight leave-one-environment-out trials or the 27 Moving Pictures sinks. We compare real-data errors using two-sided paired Wilcoxon signed-rank tests [20]. The two external-method comparisons in each single-condition real-data evaluation are adjusted together using the Benjamini–Hochberg procedure [21]. For the five-level GlobalPatterns drift sweep, the five p-values for each of the two external comparators are combined into one family of ten tests for Benjamini–Hochberg adjustment. Adjusted values are reported as q-values. Exact permutation p-values are computed from the signed-rank statistic using midranks for tied absolute differences; differences smaller than ${10}^{-12}$ are treated as numerical ties. In the real-data tables, superscripts indicate p- or q-value thresholds: ^*$<0.05$, ^**$<0.01$, and ^***$<0.001$. No superscript indicates a non-significant difference. The leave-one-environment-out tests quantify consistency across the eight GlobalPatterns environments rather than population-level uncertainty. The Moving Pictures tests treat sinks as the paired analysis units but do not model temporal dependence among repeated samples from the same body site, so those p-values are interpreted as exploratory.

5. Results

5.1. Simulated Data with Fully Observed Sources Across Taxon Counts

No phylogenetic tree is supplied to MARSH in this experiment. With all sources known and no unknown, MARSH is competitive with the best baseline on small tables and attains the lowest mean error as dimensionality grows, while FEAST and SourceID-NMF degrade (Table 1). MARSH and FastST trade the lowest mean error across the sweep: FastST leads at 50 and 500 taxa, while MARSH leads at $5,000$ and $50,000$ taxa.

5.2. Simulated Data with an Unknown Source

We first consider the simulated unknown-source setting without drift. No phylogenetic tree is supplied to MARSH for this evaluation. As the true unknown fraction increases from $0.2$ to $0.8$, MARSH attains the lowest error at every level (Table 2). FEAST’s error grows sharply as it misattributes the missing mass, while SourceID-NMF trails without leading.

We then introduce TADA drift while holding the true unknown fraction at $0.30$. A phylogenetic tree is supplied to MARSH. Its error remains near $0.15$ and is lowest at every drift level (Table 3, Figure 1(a)). The transport loss carries drifted mass back to the correct source, while the floored denominator limits inflation of the estimated unknown fraction when the residual profile lies close to the named-source mixture.

5.3. GlobalPatterns Data with an Unknown Source

We first evaluate GlobalPatterns without drift. No phylogenetic tree is supplied to MARSH in this setting. With one environment withheld as an unknown source with a true fraction of $0.30$, MARSH recovers the unknown fraction as $0.292$ and attains an order of magnitude lower $L_1$ error than the best competitor (Table 4). Two-sided paired Wilcoxon signed-rank tests give $q=0.008$ against each competitor, with all eight trials favoring MARSH.

We next introduce genus drift. A phylogenetic tree is supplied to MARSH for this evaluation. Across $f\in \{0.1,0.2,0.3,0.4,0.5\}$, MARSH’s error rises from $0.06$ to $0.28$ but remains the lowest throughout (Table 5, Figure 1(b)). SourceID-NMF degrades from $0.12$ to $0.66$ as exact matching fails, while FEAST’s error changes from $0.33$ to $0.38$. After Benjamini–Hochberg correction across all ten external-method comparisons in the drift sweep, the two-sided paired Wilcoxon tests give $q=0.0087$ for both external methods at $f=0.1$–$0.4$. At $f=0.5$, the tests give $q=0.039$ for FEAST and $q=0.0087$ for SourceID-NMF.

5.4. Moving Pictures Data with Natural Temporal Drift

We evaluate naturally occurring drift using the QIIME 2 Moving Pictures dataset [8], in which the same body sites (gut, tongue, palm) are sampled repeatedly over $\sim 15$ months. We take each body site’s day-0 community as its source and attribute each later sample (days 84–168) to its body site; the sink has drifted from its own day-0 reference through real turnover, not a synthetic perturbation. Scoring each attribution against its true body site by $L_1$ (Table 6, $n=27$ sinks), MARSH attains the lowest mean error. Sink-level two-sided paired Wilcoxon tests give $q=0.015$ versus FEAST and $q<0.001$ versus SourceID-NMF; these comparisons are exploratory because temporal dependence within body sites is not modeled. Unlike the induced genus drift, natural gut turnover is dominated by phylogenetically distant novel taxa rather than local congener movement. A phylogenetic tree is supplied to MARSH for the primary evaluation, which produces an error of $0.92$; omitting the tree produces an error of $0.99$, and the difference is not significant (two-sided paired Wilcoxon, $p=0.812$). The locality-bounded metric (Section 3.3) avoids the higher mean error observed with an unrestricted tree-Wasserstein loss, which raises MARSH’s $L_1$ error to $1.18$, while the locality-bounded loss preserves its large advantage under phylogenetically local drift (Table 3 and Table 5).

5.5. Runtime and Peak Memory on Simulated and GlobalPatterns Data

Figure 2 summarizes computational scaling on the fully observed simulated data. Without a phylogenetic tree, MARSH is the fastest method through $5,000$ taxa. At $50,000$ taxa its runtime rises to $333.31$ seconds because the coordinate-wise $L_1$ linear program grows with the number of taxa; it remains faster than FastST but is slower than FEAST and SourceID-NMF. MARSH has the lowest measured peak memory at every table size, using $149.9$ MB at $50,000$ taxa.

On the full GlobalPatterns input with $6,231$ taxa and seven candidate sources, MARSH was run with a phylogenetic tree and unknown-source estimation. It took $86.88$ seconds, compared with $29.92$–$59.49$ seconds for the external methods, while using less peak memory: $49.2$ MB compared with $131.6$–$1,852.2$ MB.

6. Discussion

The experiments show that phylogenetic information is most useful when source and sink profiles differ through local taxonomic replacement. In the simulated TADA and GlobalPatterns genus-drift evaluations (Table 3 and Table 5), tree-Wasserstein distance maps abundance discrepancies along local evolutionary paths [12,13]. The fixed locality cutoff ($h=round\left(\sqrt{n}\right)$, Section 3.3) limits the distance calculation to shallow subtree edges, allowing drifted congener counts to remain associated with their reference sources. In the Moving Pictures temporal-drift evaluation (Section 5.4), the paired comparison does not detect an error difference between fits with and without a phylogenetic tree, while MARSH retains its advantage over the external baselines. Under heavier drift, the denominator floor $\overline{D}$ prevents collapsed transport distances from inflating the estimated unknown fraction. Together, these results support tree-based optimal transport as a practical, geometry-aware basis for community deconvolution [15,16].

The results also expose several operational trade-offs. When an unobserved source is phylogenetically adjacent to a candidate reference, transport partially absorbs missing mass as local drift, leading to a modest underestimation of the unknown fraction. Integrating multi-sample temporal tracking or lineage-specific priors could resolve this boundary by separating static ancestral divergence from rapid phenotypic drift. Because tree transport assigns low cost to short paths, it can also blur fine differences between near-identical strain variants; hierarchical loss schedules that enforce exact coordinate matching at the leaves while using tree transport across internal nodes could preserve both strain-level precision and clade-level drift tolerance. The present drift evaluations use synthetic perturbations and a congener-enriched panel (92% coverage) to probe the local-transport regime. Benchmarking on field cohorts with diverse, uncurated environmental drift will further map performance across broad taxonomic scales.

The largest simulated evaluation indicates that runtime, rather than memory, is the immediate scaling constraint. A matrix-free primal–dual solver could exploit the convex weighted-$L_1$ form without constructing the full linear-program constraint matrix; such methods are established for large structured nonsmooth convex problems [22].

Future development will focus on adapting the optimization framework to incorporate relevant biological priors and information, including phylogenetic tree structure, lineage-specific evolutionary rates, sequencing depth, and temporal or environmental context. These additions could make the model more responsive to differences among microbial systems while preserving a common source-tracking framework. Evaluating MARSH across longitudinal environmental systems, including wastewater surveillance, infant gut colonization, and industrial bioreactor tracking, will test its utility across diverse ecological succession dynamics.

7. Conclusion

MARSH reframes microbial source tracking as constrained optimal transport under the tree-Wasserstein metric, coupling phylogenetic geometry with a residual-anchored unknown-source estimator. Across synthetic drift benchmarks and real 16S rRNA genus turnover, this approach attains lower mean $L_1$ attribution errors than the evaluated coordinate-wise mixture and matrix-factorization baselines at every tested drift level. Limiting the loss to local tree structure improves attribution under the evaluated feature-level drift without the higher mean error observed for unrestricted tree-Wasserstein loss under broad turnover. The denominator floor limits inflation of the unknown-source estimate, and the no-drift evaluations support recovery of withheld-source contributions. Future development will focus on adapting the optimization framework to incorporate relevant biological priors and information, including phylogenetic tree structure, lineage-specific evolutionary rates, sequencing depth, and temporal or environmental context.

Source code and scripts reproducing all tables and figures are available at https://github.com/advaitb/MARSH.