CommScores: model-informed features for microbial interaction prediction

CommScores is a Python package that computes interaction scores for pairs of genome-scale metabolic models in a defined medium, with a machine-learning layer that turns those scores into a predicted ecological class. Andrew Freiburger wrote the package and the score implementations. Chris Henry has shaped the scope, from restricting model support to ModelSEED and CarveMe reconstructions to insisting on a like-for-like comparison against SMETANA from the same genomes. The work asks whether features derived from a metabolic model add anything beyond monoculture growth and existing cross-feeding metrics. They do, and the practical result is a ranked shortlist of pairs worth measuring.

Measuring every pair in every environment is impractical

Twenty bacterial strains make 190 unique pairs. Measuring each pair on thirty carbon sources creates 5,700 pair-by-environment contexts, before varying the starting ratio, time, spatial structure or host. That experiment cannot be repeated for every new set of strains.

The alternative is to compute something from the genomes. Genome-scale metabolic models describe what an organism needs and what it can excrete, and several papers have each proposed a single-number metric for some dimension of the resulting interaction: how much two organisms compete for the same nutrients, how much one could feed the other, and how much either grows better in company. Those metrics were published separately and implemented separately. CommScores computes them in one pass over the same pair of models in the same medium, so they are comparable and a downstream model can use all of them together.

The target is the six-class ecological outcome, defined by the sign of each strain’s effect on the other: competition, amensalism, parasitism, neutralism, commensalism, mutualism.

Model-derived features beat both growth baselines and SMETANA

Four fixed feature sets isolate where the signal comes from. The learner, folds, weighting and outcome are held fixed, so only the feature information changes.

Feature set	Features	Six-class accuracy	Macro-F1	MCC
CommScores-20	20	0.51	0.36	0.34
SMETANA-only	11	0.45	0.32	0.27
SimMono, simulated monoculture growth	5	0.44	0.29	0.25
ExpMono, measured monoculture growth	10	0.41	0.27	0.21

The learner is balanced XGBoost, evaluated leave-one-carbon-out. CommScores-20 improves the F1 of every one of the six classes over SMETANA-only. Widening to all 190 order-invariant features gains little, reaching 0.530 accuracy, 0.377 macro-F1 and 0.362 MCC, so the twenty selected features carry most of the available signal.

Strain descriptors carry more of the signal than rescue or acidification

Dropping a feature family and remeasuring gives the loss attributable to that family alone. Strain descriptors, which include transport-reaction fractions, reaction-network degree, uptake capacity and dead ends, cost about 0.053 MCC when removed. That is several times more than metabolic rescue at 0.018 or acidification at 0.010. Growth response and the CommScores metabolic-potential family cost nothing measurable on their own, which reflects redundancy with the rest rather than irrelevance. Of the twenty selected features, thirteen are chosen by all five outer folds.

Most of the usable signal lives in interactions between features

The best single feature under an L2 logistic model reaches 0.135 macro-F1 and 0.192 MCC. An L1 multinomial logistic model over all twenty reaches 0.299 and 0.259. Balanced XGBoost over the same twenty reaches 0.358 and 0.337. The effect of a score depends on the partner, the carbon source and the other scores, so a linear model cannot capture it. Among tree learners XGBoost led the benchmark on macro-F1, and ExtraTrees scored slightly higher MCC.

The 0.51 figure is not comparable to published numbers from random splits

Published models on this same kChip system report about 0.88, but for a different endpoint under a random split. Matched to their protocols these models are competitive: 0.85 against a published 0.88 for one-way sign, 0.80 against 0.87 for one-way strength, and 0.81 against a published 0.71 for three-way sign. A random split lets the same pair appear in training and test under a different carbon source. The numbers above hold out a whole carbon source (LOCO, testing new media) or a whole strain (LOSO, testing new strains). The same argument, with whole clades as the held-out unit, runs through predicting what a microbe eats from its genome.

A directional model asked which of the two strains is more affected gets it right 86.6% of the time, against 82.4% for measured growth asymmetry and 65.9% for a symmetric growth lookup.

The gain depends on the carbon source

Against the measured-monoculture baseline, the improvement in Cohen’s kappa is large on ribose, maltose and xylose, moderate across most of the panel, and negative on uridine and serine. The classes themselves lie on a competition-to-cooperation continuum rather than in discrete clusters, with resource overlap rising as predicted cross-feeding falls, which is part of why exact class calls are difficult.

Coarser groupings are called more accurately than the six-way label

The same probabilities can be reported at three resolutions: exact six-class at 50.9%, four neighbourhood groups merging the harmful and the beneficial neighbours at 67.8%, and three sign groups at 76.4%. Reporting at the resolution the call supports is more useful than forcing a six-way label.

Filtering on confidence buys accuracy by keeping the majority class

Keeping only predictions with a maximum class probability of 0.70 or above retains 1,058 of the 4,407 observations, or 24% of them, and lifts exact accuracy from 51% to 74%. Competition, however, rises from 44% of all predictions to 84% of those that survive, and 83% of the retained true Amensalism cases are still called Competition.

Transfer to a second dataset failed

A kChip-scale experiment is not available for most systems, so the models were tested against AtLeaf. Exact six-class accuracy was 0.35, against 0.40 for the BIT rule, on twenty matched observations containing a single positive outcome, and the CommScores model ranked that one positive 18th of 20. Repairing the AtLeaf label pipeline moved unique pairs from 26 to 35 and positive labels from 4 to 5, which gives a sense of how thin the evaluation is. This is domain shift rather than successful transfer, and it is reported as such.

The top decile is three times richer in positives than random

Ranking pairs by predicted positive probability yields 48.1% positives in the top decile against a 15.7% random baseline, an enrichment of 3.1×, where both monoculture baselines reach 39.7%. The output is a ranked, mechanism-filtered shortlist: the candidate universe, then growth-feasible pairs, then pairs with high positive probability, then a mechanistically diverse selection. Cross-feeding breadth scales with the cooperativity of the class, so the SMETANA annotations can filter the shortlist further and attach a directional hypothesis to each pair. This supports choosing the next experiment in a context-matched pairwise setting, and says nothing about communities of more than two.

Seven steps from a genome to a prediction

A genome becomes a prediction in seven steps: reconstruct a genome-scale model, gap-fill it, simulate the kChip carbon environments, compute the scores, select the CommScores-20 feature set, and fit a balanced XGBoost against measured kChip coculture outcomes, which supply both the training labels and the evaluation.

The score set describes resource competition, cross-feeding, biosynthetic support, growth response, functional similarity and the metabolites behind each interaction. These outputs are combined with simulated growth, SMETANA dependency, acidification and strain descriptors. Because a pair has no canonical first member, every modelling feature is symmetric across the two strains.

The training table is 4,407 pair-by-carbon observations, the survivors of 190 pairs across 30 retained carbon environments after filtering and feature joins. It is lopsided, with 1,762 Competition, 967 Parasitism, 759 Amensalism, 503 Commensalism, 229 Neutralism and 187 Mutualism, and every number above carries that skew.

The modelling stack combines ModelSEEDpy, COBRApy and MSCommunity with KBase genomes and media. Its output is both a prediction table and the metabolites behind each score, so a highly ranked pair comes with a mechanistic hypothesis for the experiment.