<figure><img src="image_1.png" /></figure>

HAL
open science

# A machine learning framework for the prediction and analysis of bacterial antagonism in biofilms using morphological descriptors

Raphaël Rubrice, Virgile Gueneau, Romain Briandet, Antoine Cornuejols, Vincent Guigue

<figure><img src="image_2.png" /></figure>

## ► To cite this version:

Raphaël Rubrice, Virgile Gueneau, Romain Briandet, Antoine Cornuejols, Vincent Guigue. A machine learning framework for the prediction and analysis of bacterial antagonism in biofilms using morphological descriptors. Artificial Intelligence in the Life Sciences, 2025, 8, pp.100137. 10.1016/j.ailsci.2025.100137. hal-05217869

HAL Id: hal-05217869

https://hal.inrae.fr/hal-05217869v1

Submitted on 25 Aug 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_4.png" /></figure>

Research article

# A machine learning framework for the prediction and analysis of bacterial antagonism in biofilms using morphological descriptors

Raphaël Rubrice a,b,c,e f,g,*, Virgile Gueneau a,b,d,f,g, Romain Briandet a,b,d,g,
Antoine Cornuejols a,b,c,e, Vincent Guigue a,b,c,e,*

## ARTICLE INFO

Dataset link: https://github.com/raphaelrubric
e/BiofilmAntagonismPrediction, https://entrep
ot.recherche.data.gouv.fr/dataset.xhtml?persis
tentId=doi:10.57745/XRXQEI

Keywords:

Biofilm

Beneficial biofilms

Pathogens

Machine learning

Artificial intelligence

Explainability

Interpretability

## ABSTRACT

Biofilms are structured microbial communities that promote cell interactions through close spatial organization, leading to cooperative or competitive behaviors. Predicting microbial interactions in biofilms could aid in developing innovative strategies to prevent the colonization of undesirable bacteria. Here, we present a machine learning approach to predict the antagonistic effects of beneficial bacterial candidates *Bacillus* and *Paenibacillus* species against undesirable bacteria (*Staphylococcus aureus*, *Enterococcus cecorum*, *Escherichia coli* and *Salmonella enterica*), based on the morphological descriptors of single-species biofilms. We trained the models using quantitative features (e.g. biofilm volume, thickness, roughness, or substratum coverage). As a proxy for antagonism, an exclusion score was used as the supervised training target. The latter was calculated based on the ratio of biofilm volume between the undesirable bacteria and the beneficial strain. We subsequently applied diverse explainability methods to analyze the resulting model and found insights highlighting the importance of biofilm formation context when predicting antagonism. Our results demonstrate that machine learning offers an efficient, data-driven tool to predict microbial interactions within biofilms and support the selection of competitive beneficial strains against pathogens. This approach enables scalable screening of microbial interactions, making it applicable to both research and biotechnological applications.

## 1. Introduction

Constructing scientific models has traditionally played an impor-
tant role in scientific discovery. Whether these models are theoretical,
experimental, or numerical, they all serve as essential tools for compre-
hending complex phenomena. Biofilms are examples of such complex
collective biological life, defined as spatially organized microbial com-
munities embedded in an extracellular matrix [1]. They colonize the
majority of ecosystems and can account for up to 80% of microbial life
forms [2].

A challenge in the field of microbial ecology is understanding the
complex dynamics of microbial interactions within biofilms, which
is crucial for better predicting their behavior. These interactions can
be cooperative, where microorganisms collaborate to exploit resources
or enhance survival. Conversely, they can be competitive, leading to

resource competition, predation, or the production of inhibitory com-
pounds [3]. The balance between cooperative and competitive inter-
actions shapes biofilm behavior, influencing critical aspects such as
growth rates, spatial distribution, and community stability. Deciphering
these intricate relationships is essential for developing effective strate-
gies to manage biofilms, which can act as reservoirs for undesirable
microorganisms.

In the One Health framework, measures to limit the use of chemical inputs such as biocides and antibiotics are encouraged. In this context, the application of microbial solutions able to guide surface microbial communities, through the formation of beneficial biofilms that prevent undesirable bacterial colonization, are considered complementary to chemical treatments, thereby supporting their reduced application [4]. Since all undesirable bacteria are pathogens in this study, they will be referred to as pathogens.

<footnote>a Paris-Saclay University, 3 rue Joliot Curie, Gif-sur-Yvette, 91190, France</footnote>

<footnote>b INRAE, France</footnote>

<footnote>c AgroParisTech, 22 place de l'Agronomie, 91123, Palaiseau, France</footnote>

<footnote>d Micalis Institute, Domaine de Vilvert, 78352, Jouy-en-Josas, France</footnote>

<footnote>e Department of Applied Mathematics and Informatics, AgroParisTech, France</footnote>

<footnote>f Lallemand SAS, 19 rue des Briquetiers, 31700, Blagnac, France</footnote>

<footnote>g LabCom, Biofilm1Health, France</footnote>

<footnote>* Corresponding authors.
E-mail addresses: raphael.rubrice@agroparistech.fr (R. Rubrice), vincent.guigue@agroparistech.fr (V. Guigue).</footnote>
The development of these microbial solutions requires the screening of beneficial strain libraries to assess their potential to establish competitive biofilms capable of reducing pathogen colonization. This requires characterizing each potential partner and then studying them in multi-strain biofilms to examine their interactions in pairs (candidate beneficial strain, pathogen). A large number of combinations (i.e. pairs) are thus to be examined. Within the experimental approach, some recent studies use confocal laser scanning microscopy (CLSM) to examine the spatial arrangement of multi-species biofilms without disrupting their structure, identifying beneficial candidate bacteria capable of excluding pathogen through biofilm formation [5]. However, these studies are time-consuming, costly, and require specialized equipment and expertise to analyze biofilms with CLSM coupled with image analysis [6].

In parallel, within the simulation approach, traditional biofilm mechanistic models like individual-based techniques or parameterized sets of equations used to model biofilm growth are quite limited when it comes to predicting interactions within biofilms [7]. Indeed, biofilms are highly complex and dynamic systems, and trying to find a set of equations or rules to simulate their inner workings and interactions is no trivial task.

Given these difficulties and the wealth of data that becomes available, artificial intelligence (AI) technologies are increasingly used to predict biological phenomena [8]. For instance, models obtained through machine learning have shown that interactions between cultivable bacteria can be predicted based on the growth of individual species using fresh-state growth data, phylogenetic distances, and metabolic profiles [9]. Similarly, microbial community profiles with plant-protective properties have been identified through machine learning [10]. By choosing the data-driven approach we want to leverage interplays between descriptors that would have been missed by mechanistic models [11]. Beyond predicting antagonism, we aim to quantify the impact of each morphological parameter on interaction. We hypothesize that single-species biofilms can be represented by specific morphological descriptors, allowing each biofilm characteristic to predict outcomes when strains are mixed.

Thus the following research question arises: Can beneficial biofilm antagonism against undesirable bacteria be accurately predicted using morphological parameters derived from single-species biofilms ?

The hypothesis that this would be possible is based on the fact that biofilm characteristics play a key role in bacterial interaction and resistance.[12-16]. To answer this question, this study aims to employ a machine learning approach to predict the antagonistic potential of beneficial biofilms against undesirable bacteria in multi-species communities. The quantitative descriptors used for prediction are derived from image stacks, including biofilm volume, thickness, and surface covered [17]. This enables the identification of beneficial candidates without requiring the study of all possible combinations of the strains. Thus, the experimental effort could be dramatically reduced. However, the high noise levels inherent in biological data, arising from both natural variability and experimental conditions, make it particularly challenging to accurately predict the outcomes of interactions between biofilms and pathogens.

This article compares different state-of-the-art machine-learning models, different protocols for filling missing data, processing replicate experiments, and attribute selection schemes to explain predictions and focus on the most promising signals. It does so by proposing a specific cross-validation protocol to better control biological homology. Given the nature and size of the data, we focus this study on machine learning techniques rather than deep-learning [18].

The main advantages and novelty of our approach are demonstrated through the following aspects: (1) A data-driven methodology is utilized, eliminating the reliance on predefined rules or equations. (2) The capability to manage multiple and heterogeneous biofilm parameters as inputs is ensured. (3) Insights into the importance and contribution of

each biofilm parameter to the prediction are provided. (4) A representative benchmark of linear and non-linear methods for this novel task is offered, providing an initial answer to the previously stated research question. (5) Offers the potential for learning from both experimental and simulated data.

By integrating experimental data with machine learning, we advance current biofilm research by introducing a powerful approach for studying and predicting microbial interactions in multi-species communities. Furthermore, this approach supports the targeted selection of bacterial strains for beneficial biofilm applications.

## 2. Material and methods

### 2.1. Data used

For each pair (Bacillus strain, pathogen), the authors of [19] plated mixed species biofilms. They also cultivated single-species biofilms as references for each organism. Two inoculation models were used to study interactions in multi-species biofilms by confocal laser scanning microscopy (CLSM). In the co-inoculation model, an equal ratios of both pathogenic and beneficial strains were determined by image analysis of CLSM observations. Following adhesion with these pre-determined ratios, biofilms were incubated 24 h before being observed with CLSM. In the recruitment model, a planktonic culture of the pathogen was added to a pre-established biofilm for 24 h prior CLSM observation. After adhesion, an incubation of 24 h was followed by CLSM observation. A staining strategy was used to distinguish candidate beneficial bacteria from pathogens. The entire biofilm was stained using SYTO61 to visualize the entire population in red (i.e. pathogens and beneficial strains), while the pathogen expresses GFP. Subtracting the green signal from the red signal allows for the estimation of the proportion of beneficial bacteria. In the present study, the co-inoculation (24 h), recruitment (0 h), recruitment(24 h) are respectively called mixing-models I, II and III.

To evaluate the effectiveness of Bacillus strains against various pathogens, we quantified the GFP biovolume (µm³/µm²) in submerged mixed biofilms, then normalized these values relative to the biovolume of biofilms of the target pathogen labeled with GFP. A normalization system was developed to represent the anti-pathogenic activity of biofilms on a standardized scale. This approach used a min-max transformation to scale the data between 0 and 1, following the formula:

$$x_{\text{normalized}} = (x - x_{\text{min}})/(x_{\text{max}} - x_{\text{min}}) \quad (1)$$

Consciously, these values were then inverted to obtain an anti-pathogenic activity score where 0 indicates the lowest activity and 1 represents the highest activity. This method allows for direct comparison between different Bacillus strains, facilitating the identification of the most promising candidates for potential applications.

This dataset provides information about:

* 18 *Bacillus* strains

* 4 pathogen species

or 18 × 4 = 72 pairs. For each of these pairs, there are 9 measurements corresponding to:

* the 3 methods underlined above for mixing 2 bacteria (a *Bacillus* and a pathogen) together in a biofilm

* and 3 biological replicates for each experiment

which all together gives 18 × 4 × 3 × 3 = 648 exclusion scores from real mixed biofilms (see [20]).

Because control samples were necessary for each pair, in addition to mixed-species biofilm, access to morphological descriptors of single-species biofilms is given. More precisely, in the dataset, there exist between 12 to 19 samples of single-species biofilms for each *Bacillus*.
## Biological data collection

<figure><img src="image_5.png" /></figure>

## AI model training

<figure><img src="image_6.png" /><figcaption>Bacillus spp.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>undesirable bacteria</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Bacillus spp. + undesirable bacteria</figcaption></figure>

**Morphological descriptors**

* biofilm volume

* substratum coverage

* thickness

* height

* roughness

* ...

**Morphological descriptors**

* biofilm volume

* substratum coverage

* thickness

* height

* roughness

* ...

biofilm volume of the undesirable bacteria

normalization to biofilm volume of undesirable bacteria alone

interaction score calculation

low undesirable bacteria exclusion by the *Bacillus* spp.
strong undesirable bacteria exclusion by the *Bacillus* spp.

**Objective:**

Interaction score prediction

using a machine learning model
on morphological descriptors of
mono-species biofilms

<figure><img src="image_9.png" /><figcaption>**Fig. 1.** Representation of the data used and the task at hand. The goal is to use morphological descriptions of single-species biofilms to predict the antagonism (represented by an interaction score close to 1) of the *Bacillus* strain regarding the pathogen.</figcaption></figure>

Respectively, there exist 36 samples for each pathogen (because for each bacteria, 12 samples were put in culture for each of the 3 mixing methods experiments).

Descriptors extracted from both single-species *Bacillus* and pathogen data will be used as input to our models. For each *Bacillus*/pathogen pair, the goal will thus be to map these single-species biofilm descriptions of both organisms to the exclusion score reflecting antagonism between them in real mixed-species biofilms.

The characteristics of each biofilm were obtained using the BiofilmQ software using a dedicated method [17].

If we can reach satisfying performance for this task, the learned model could enable a much faster selection of strains to find candidates capable of producing beneficial biofilms. Each high predicted exclusion score would indicate strains with strong potential antagonistic properties.

## 2.2. Method

The goal is to learn to predict an interaction score for each pair: (Bacillus, pathogen). However, our data set contains replicas for each organism. The question then arises as to how to use these replicas to construct a training set. There are several possible choices. In this study, we tested three methods.

In the first one, one replica was randomly selected for each of the two organisms paired. For example, the interaction between the *Bacillus* strain B1 and the pathogen *E. coli* was described by one example

made by the combination of one randomly drawn replica of *E. coli* and one randomly drawn replica of the B1 strain. Both are then linked to one randomly selected interaction score experimentally measured for this pair. This reflects the intrinsic variability that can be encountered when building the input space for a new pair.

The second method consists of using averaged and standard deviation representations. This approach is implies calculating the mean ($\mu_i$) and standard deviation ($\sigma_i$) of each feature $i$ for each organism. From this, three representations are created: one corresponding to the mean $\mu_i$ for each feature $i$, one corresponding to the $\mu_i + \sigma_i$ value for each feature $i$, and the last one to the $\mu_i - \sigma_i$ value for each feature $i$. Then, each of these three representations is associated with one exclusion score taken from the three available replicates, so that the same score is not shared between the three synthesized pairs. This method produces three training examples for each pair (*Bacillus*, *pathogen*) providing summarized information about the possible variability to the learning system.

The third method consists of describing all possible combinations of each pair (*Bacillus*, *pathogen*) in the training set. For example, given that there are $n_{B1}$ replicas for the B1 *Bacillus* and $n_{E_c}$ replicas for the pathogen *E. coli*, and given that three mixing methods were used in the experiment, this gives $n_{B1} * n_{E_c} * 3$ combinations, and thus lines, in the training set. This exhaustively explores the intrinsic variability that can be encountered when building the input space for a new pair.

This dataset is considerably larger than the two previous ones as we reach over 100,000 lines in the resulting data table.
<figure><img src="image_10.png" /><figcaption>Fig. 2. Proportion of non-missing values. “Modele” corresponds to the mixing method used to form the mixed-species Biofilms. Computed across all available samples.</figcaption></figure>

## 2.3. Detailed dataset description

Another issue concerns the management of missing values (NAs). Indeed, many learning methods cannot cope with missing values in the inputs. Each biofilm is described using 4 specific parameters:

* Height

* Volume

* SubstratumCoverage (area covered by the biofilm at substrate level)

* Roughness

(For more details on each parameter calculation, please refer to the following website https://drescherlab.org/data/biofilmQ/docs/ and paper [17]).

The statistics summed up in Fig. 2 show that both Substratum Coverage and Roughness had missing values. Moreover, **Roughness was never acquired on pathogen biofilms** which is why it is not presented in the plot. The only categorical feature, **(named “Modele” on all figures)**, specifies the mixing model. It was One Hot Encoded [21] thus using three (boolean) features corresponding to the three mixing methods.

## 2.4. Evaluation protocol

Performance evaluation must be carried out on a test set that differs from the training set, otherwise, the predictive performance of the learned model may be overestimated, corresponding to what is called over-fitting of the training data. The way training and testing samples are separated is key when developing models.

### 2.4.1. Homology-controlled cross-validation

To that end, it is important to note that homology exists in every biological data because of common ancestors, living conditions, or functional homology of processes. The homology signal is thus often trivial and significant between data points. Using traditional random train/test split poses the issue of homology learning rather than causal approximation by the model. Indeed, because of this usually strong signal, machine learning algorithms can easily learn homology-based relationships instead of more causal ones. This is an issue because such models will only perform well on data points that share homologies with training data which limits the generalization power of both the model and the insights extracted from it.

In our approach, the goal is to achieve a model that could be used on new strains of positive bacteria and undesirable ones. These new organisms will likely be of various phyla and may not share as much homology as there is in our current dataset. Thus, controlling homology to train our models becomes crucial during the evaluation phase to identify which models are most promising for generalization.

Two types of homology need to be taken into account: organism-wise and interaction-wise. **Organism-wise homology** can be controlled by training on every interaction except the ones involving a specific organism and by then testing the model on interactions involving the specifically ignored organism. This needs to be done for all organisms **Interaction-wise homology** between Org1 and Org2 is thus controlled by removing Org1 and Org2 from the training set and evaluating the model exclusively on interactions between Org1 and Org2. This needs to be done for all possible pairs.

To measure the prediction performance, we use an evaluation protocol based on cross-validation where folds are designed to control homology instead of random train/test splits. It will be referred to as **Homology-Controlled Cross-Validation (HCCV)**. Since there are 18 *Bacillus* strains and 4 pathogens in the dataset, there are 18 + 4 hold-out sets for organism-wise homology and 18 + 4 hold-out sets for interaction-wise homology which leads to a total of 94 folds (see Fig. 3). This protocol allows: (i) to control for homology, (ii) to compute an overall generalization error approximation (metric average overall folds just as in classical cross-validation [22]) and (iii) to carry out a thorough analysis of performances and feature importances.

A weighting strategy will be applied when computing the average performance to take into account test set sizes as all different hold-out sets are not equal in size. This strategy consists of giving more weight to larger test sets because the metric evaluation will be more robust as the test set size increases. The weight applied to each fold i is $\frac{n_i}{\sum_{j=1}^{n_i} n_j}$ where F is the number of folds (94) and n_i the test set size for fold i.

### 2.4.2. Evaluation metrics

Regarding evaluation metrics, as is generally done in regression tasks, we used the Root Mean Squared Error (RMSE) score and the Mean Absolute Error (MAE) score.

The RMSE score formula is:

$$RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2} \qquad (2)$$

Where $y_i$ is the real score value and $\hat{y}_i$ is the predicted value. As for the MAE score, the formula is as you would expect:

$$MAE = \frac{1}{n} \sum_{i=1}^{n} |y_i - \hat{y}_i| \qquad (3)$$

With $y_i$ and $\hat{y}_i$ being the same as for the RMSE formula. MAE is the most interpretable metric out of the two but RMSE has the advantage of penalizing more large errors than MAE. Because both metrics are valuable, some results will be obtained by referring to the average of both metrics. This will be called the Cross Mean (for mean across metrics) and is simply the average of both weighted RMSE and MAE scores. Recall that the target variable for this task (the exclusion score), takes its values between 0 and 1 thus this is the expected range for acceptable errors for our models.
<figure><img src="image_11.png" /><figcaption>Fig. 3. Comparison of our evaluation protocol, homology-controlled cross-validation (or HCCV) with a classical 5-fold cross-validation (CV) scheme.</figcaption></figure>

### 2.4.3. Confidence intervals

Regarding this evaluation protocol, the last important piece of information to convey is the method of computation of confidence intervals. We essentially are evaluating 3 similar but still different tasks: predicting interactions involving new *Bacillus*, new pathogen, or both (respectively *Bacillus*, *Pathogens*, *Interactions*). These may or may not present the same distributions and because we do not wish to make strong assumptions on the distribution of evaluation metrics across folds, classical methods based on normal or Student's t distributions were considered not appropriate. We therefore chose the widely used Bias-Corrected acceleration (BCa) method which leverages bootstrap and bias correction to compute confidence intervals [23]. This was done using the *arch* library [24].

### 2.5. Model selection

Once the evaluation protocol and metrics are defined, the next step is to choose which model to use. We tested several models on the same datasets. During this phase, we aim to identify the most promising algorithm to achieve satisfying predictive power. The one with the best performance will be kept for further optimization.

### 2.6. Performance optimization

After Model selection, we performed several strategies to try performance optimization. We tried to identify which combination of hyperparameters led to the best performances. We also conducted feature selection and feature engineering with the goal of reaching better performances.

### 2.7. Model analysis

As our objective is to have a somewhat understandable model, we thoroughly analyzed our models. Specifically, we analyzed **feature importances** of the model obtained after applying our protocol. Details on how feature importance was obtained will be given later on. Another method to identify key features is to make an **ablation study**. This allows one to very clearly demonstrate how each feature or category of features contributes to making a better model. We also performed **error analysis** to try to understand where, why, and how much the model is making mistakes. This is truly important as it can allow a better judgment of the model's capabilities and sometimes reveal patterns in the errors made. If such patterns are found it might then be possible to address them and correct the model. Finally, the use of **SHAP values** to explain machine learning predictions and model behavior is widely accepted as a valuable approach to explainable machine learning, we thus used the *shap* library to do so [25–27].

## 3. Model formulations and specificities

We consider the following statistical learning problem: given a dataset $X \in \mathbb{R}^{n \times p}$, where each row $X_i = (X_{i, bacillus}, X_{i, pathogen})$ represents a vector of values for *p* single species biofilm parameters (height, volume, etc.) for a pair of *Bacillus* and pathogen strains, and a vector $Y \in \mathbb{R}^n$, where each element $Y_i$ represents the exclusion score computed from image analysis of biofilms for the **same pair of strains when mixed in reality**. We seek to estimate a function $f: \mathbb{R}^p \to \mathbb{R}$ such that $f(X_i) \approx Y_i$ for all $i \in \{1, \dots, n\}$. For simplicity, we make the assumption that the same function $f$ is good for predicting the exclusion scores in all three mixing methods (see Section 2.2) which may not be the case in reality.

To estimate the function $f$, which predicts an interaction score from the descriptive features of biofilms, we have tested several linear hypotheses learned using Ridge Regression, Lasso Regression and Elastic Net in addition to non-linear ones using models such as the Support Vector Regressor (SVR), Gradient Boosting Trees (GBT) and Random Forest (RF) [28–30]. We then compared their performances and characteristics.

Linear models (classic Linear regression, Ridge Regression, Lasso Regression) are models that fit a linear function of the form $f(X) = X\beta + \epsilon$, where $\beta \in \mathbb{R}^p$ is the vector of parameters to estimate, and $\epsilon \in \mathbb{R}^n$ is the vector of errors. The estimation of the parameters is done by minimizing a cost function $J(\beta)$, which is a combination of the sum of squared errors and a penalty on the norm of $\beta$ (for Ridge and Lasso). The penalty allows one to reduce overfitting and improve the fitting stability of the model. The form of the penalty depends on the model: Ridge uses a $L_2$ penalty, and Lasso uses a $L_1$ penalty. The cost functions are the following:

$$J_{\text{Ridge}}(\beta) = \frac{1}{n} \sum_{i=1}^{n} (Y_i - X_i \beta)^2 + C \|\beta\|_2 \quad (4)$$

$$J_{\text{Lasso}}(\beta) = \frac{1}{n} \sum_{i=1}^{n} (Y_i - X_i \beta)^2 + C \|\beta\|_1 \quad (5)$$

Where $C$ is a hyperparameter that controls the degree of regularization. In ElasticNet, both regularization methods are used to learn the regression estimator [31–33].

Non-linear models can capture non-linear relationships between the variables and are often performing better than linear models on complex problems. The SVR uses a kernel function to transform the data into a feature space of higher dimension, where it is easier to find a linear regression function. GBT combines several decision trees into a model and uses a gradient descent optimization algorithm to learn the ensemble estimator. Unlike RF, the trees in GBT models are not built independently. Each new tree in the ensemble is built to correct the errors of the previous ones. This iterative process results in a meta-model which, often, minimizes prediction errors effectively. By contrast, the RF builds and combines several strong decision trees where all trees are built independently and in parallel.
<figure><img src="image_12.png" /><figcaption>Fig. 4. Model selection campaign results. Lower is better. All error bars are overlapping. However, on average, ensembles of trees stand out as the best candidates. LGBMRe regressor on the combinatoric dataset and RandomForest Regressor yield the best performances when evaluated respectively using weighted MAE and RMSE scores. Error bars are BCa 95% Confidence Intervals (CI) computed across all HCCV folds. All models were used after applying Mean Imputation and Standard Scaling to input variables.</figcaption></figure>

Mathematically, the prediction of a single decision tree T for an input vector x can be written as:

$$T(x) = \sum_{l=1}^{L} c_l \mathbb{I}(x \in R_l) \qquad (6)$$

Where L is the number of terminal nodes (leaves), R_l is the region of the feature space corresponding to the path to the lth terminal node, c_l is the constant value associated with the lth terminal node, and I is the indicator function.

The prediction of the XGBRe regressor is then the weighted sum of the predictions of the N trees:

$$GBT(x) = \sum_{i=1}^{N} w_i T_i(x) \qquad (7)$$

where T_i(x) is the prediction of the ith decision tree for the input vector x, and w_i is the weight of the ith tree, which is learned during the training process. In GBT, weights w_i are determined by the gradient boosting algorithm, which aims to minimize the loss function by adding new trees that correct the residual errors left by the existing ensemble [29].

Compared to linear methods, some of the non-linear models can present limitations, such as increased computational time, especially when dealing with large datasets or complex tree structures (an issue mainly for SVR models defined with non-linear kernels).

Both linear and non-linear model implementations are taken from existing powerful libraries [34–37].

## 4. Results

In this section, we will present results to motivate dataset, model, and preprocessing choices followed by a performance optimization campaign (feature selection, feature engineering, and hyperparameter grid search).

### 4.1. Dataset and model choice

First, we need to determine which dataset construction method (see Section 2.2) and which model is going to be best suited for the task at hand. To do that we validated each model on each dataset using the protocol described in Fig. 3. Random sampling, 'averaged +/- standard deviation' and combinatoric datasets are respectively noted as 'random', 'avg' and 'combinatoric'.

Results shown in Fig. 4 lead to conclude that there is no significant difference between model performances as all 95% Confidence Intervals (CI) are quite large and overlapping. However, by considering only the average, we can note that ensemble of tree-based methods tend to outperform linear ones, especially on the average and combinatoric datasets. Within ensemble methods, results are even closer between both types of methods depending on the dataset construction used, and no clear implementation of ensemble models is identifiable.

Using MAE as an evaluation metric, the best configuration is obtained using the LGBMRe regressor (GBT) on the combinatoric dataset with a weighted score of 0.155. However, using RMSE, the best performances are achieved using the RandomForest Regressor on the average dataset with a weighted score of 0.205. Recall that the target variable
for our learning task ranges from 0 to 1. Thus, results show acceptable errors: our machine learning models did learn patterns between input features and the output variable.

The average dataset is a simplified version of the problem whereas the combinatoric dataset is an exhaustive display of potential biological variability. Based on that, we decided to value biological variability over simplicity. For this reason, we chose the best ensemble model on the combinatoric dataset: LGBMRegressor. Further details on why this choice was made can be found in Appendix A.

## 4.2. Preprocessing method selection

Handling of missing values and normalization was thoroughly investigated by evaluating the LGBMRe regressor using a variety of Imputer/Normalizer strategy pairs (cf. Appendix B). We found that different strategies lead to similar performances with only slight changes in both model variance and performance (the same random state was used for all runs which were conducted on the same device). We chose K-nearest neighbors imputing and the Robust scaling strategy for a compromise between simplicity and robustness to outliers.

## 4.3. Input feature space and hyperparameters optimization

In this section, we apply feature selection and engineering as well as hyperparameter optimization. Note that we distinguish native features (original features) and feature-engineered ones.

### 4.3.1. Native feature selection

After setting the dataset, model, and preprocessing choices, we were able to start model optimization. First, we tried to identify the best combination of features. Because the HCCV protocol involves a lot of folds and because the combinatoric dataset is large (115,020 rows), each experiment is computationally costly. Thus, instead of using ablation studies to determine which features are useful and which can be removed at each selection step, we used permutation feature importances (PFI) to select such features.

Indeed, the latter has the advantage of not requiring to re-train a model to compute the feature importance of a single feature. The PFI of a feature *j* is defined as $ε_j - ϵ$ where $ε_j$ is the error of the model when shuffling values of feature *j* and $ϵ$ is the error the model makes when using normal features. By randomizing feature values we break the link between the feature space and the target distribution for a specific feature. This means that if the feature is useful for predictions, breaking this link will lead to worse performances. The opposite, for a feature that is of no use to the model, means that no matter the value this feature takes, no impact is observed on the predictions which will lead to little to no changes in model performance. A truly useless feature is defined as having a PFI of 0 or below [30,38].

During this process, we compute the average of both metrics (weighted RMSE and MAE, named cross mean) to have a single score to base our decisions on.

Fig. 5 shows that after the first step, substratum coverage for the pathogens was removed from the input space since its PFI was the lowest at 0.001. During step 2, a model without this feature was trained and used as baseline (0.1849 Cross Mean) and the new least important feature was biofilm roughness with a PFI of 0.002. Finally, after the last step, the reference was a model without biofilm roughness. In this final step, the cross mean increased by 0.0002 going to 0.1851 compared to the last step. Though this change was small it still constituted an overall increase in error. For that reason, the process was stopped and the only feature that will be removed is the substratum coverage of pathogens.

It can be noted that the 'Modele' and 'pathogen Height' variables are statistically significantly more important than all the others, with the former being the most important feature by quite a margin. This interesting result will be further discussed in Section 5 on Model Analysis.

<figure><img src="image_13.png" /><figcaption>Fig. 5. Native Feature selection. Permutation feature Importances (Permuted - Reference) for each feature. At each step, we remove all features with PFI strictly inferior to 0.001. If no feature satisfies this condition, we remove the feature having the lowest PFI. At step *i* > 1 the baseline is a model trained on all features except those removed following the *i* - 1 previous steps. Cross-Mean refers to the average of weighted RMSE and MAE scores.</figcaption></figure>

### 4.3.2. Feature engineering and selection

We then tried to obtain better performances by adding engineered features based on the selected native variables. The following were added: Products, Differences, Powers of 2 and 3, Density, and FakeScore.

Density corresponds to the ratio between the Height and the Volume of the Biofilm. It was calculated for each *Bacillus* and each pathogen. FakeScore is an attempt at constructing a proxy of the real exclusion score by using input space features and reflecting on the nature of the target score. It corresponds to the following formula :

$$ \text{FakeScore} = 1 - \frac{P_{\text{BiofilmHeight}}^3}{P_{\text{BiofilmVolume}}} \qquad (8) $$

Results are presented in Fig. 6.

Following the same logic as for native feature selection, out of all the engineered and native features, the final set of input variables is
Top 10 Feature Importances with 95% CI - Step 1 - Baseline Cross Mean: 0.1816

<figure><img src="image_14.png" /></figure>

Top 10 Feature Importances with 95% CI - Step 2 - Baseline Cross Mean: 0.1804

<figure><img src="image_15.png" /></figure>

Top 10 Feature Importances with 95% CI - Step 3 - Baseline Cross Mean: 0.1816

<figure><img src="image_16.png" /><figcaption>Fig. 6. Top 10 features (ranked by PFI) at each step. More detailed results are given in Appendix C.</figcaption></figure>

made of: *Bacillus* and Pathogen Height, *Bacillus* SubstratumCoverage, *Bacillus* Volume, Modele, product of Height and SubstratumCoverage for *Bacillus*, product of Height and Volume for *Bacillus*, product of SubstratumCoverage and Volume for *Bacillus*, difference between Height and Roughness for *Bacillus*, difference between *Bacillus* Roughness and Pathogen Height, *Bacillus* Volume cubed, *Bacillus* Density and *Bacillus* FakeScore.

Using this combination of features yielded a Cross Mean of 0.180 which is an improvement compared to using only the best native features. These results show the importance of feature engineering and selection when building a machine learning model and highlight the usefulness of adding interaction and domain-specific features.

4.3.3. Hyperparameter search

Our last step for performance optimization was to identify the best hyperparameters for our model on the combinatoric dataset. Using the optuna library we conducted a campaign of 50 trials to identify the best hyperparameters by exploring ranges of regularization, bagging, and tree structure hyperparameters. This search led to a marginal increase in performance upon using the found hyperparameters going from 0.184 to 0.1787. Details about the campaign results are available in Appendix D.

5. Model analysis

After the optimization campaign, we now have determined that the best performances are reached using the LGBMRegressor with searched
<figure><img src="image_17.png" /><figcaption>Fig. 7. Ablation study on individual and grouped features using MAE score. A gradient of red to green is used to respectively denote the increase or decrease in error.</figcaption></figure>

parameters on the combinatoric dataset with a refined set of descriptors for the *Bacillus* and the pathogen. A growing concern in the field of Machine Learning and AI is related to model explainability which is key when proposing data-driven approaches to any industry. In this section, we will discuss our model and analyze its properties using diverse explainability methods.

## 5.1. Feature importances

First, we need to understand feature importances, i.e. which features are the most important to achieve satisfying predictive power. This is investigated through the previously presented PFIs and an ablation study to further confirm the analysis.

### 5.1.1. Permutation feature importances

The method of Permutation Feature Importances described in Section 4.3.1 to select features already gave a substantial insight into feature importances in our dataset. The mixing model and the Height of pathogen biofilms stand out as the most crucial ones (see Fig. 6 and Appendix C).

The identification of the mixing method (Modele) as most critical for enhancing prediction accuracy, reflects the importance of experimental conditions. Pathogen biofilm height being the second most important also is striking as it shows that structural information can indeed be linked to antagonism and thus, underlying mechanisms. This insight is consistent with prior experimental observations on biofilm interactions [12–16]. Although other features have modest PFIs, this does not mean that they are useless, but rather suggests they contribute to incremental improvements rather than driving major prediction shifts or that they are important for specific cases.

Overall, this initial analysis brings interesting results as it shows that the model uses features in a way that is both insightful and reflects experimental observations.

### 5.1.2. Ablation study

To complement the feature importances computed via PFIs, an ablation study (see Fig. 7) was conducted to further confirm the model's reliance on key features.

Results are consistent with the previous analysis as removing the 'Modele' feature severely impaired performance by increasing the absolute error on average from 0.154 to 0.198. Moreover, the second most impactful ablation consists of removing all input features related to pathogen descriptors (Height, the difference between Height and *Bacillus* Roughness), and the third exaequo is the sole removal of 'pathogen Height'. These results strongly confirm our previous conclusions.

Interestingly, complete removal of either all *Bacillus* or pathogen parameters also plummets performances. This shows that it is truly the interaction between both pathogen and *Bacillus* biofilm parameters and the mixing-condition features that allow good predictions. Finally, the range of 95% Confidence Intervals demonstrates a high variability in feature importances across validation folds, confirming that different organisms and interactions showcase different antagonism mechanisms.

Those results demonstrate that the model captured biologically interesting insights as it takes advantage of all different features in a way that reflects experimental observations and theoretical expectations.

### 5.1.3. Does the approach yield good estimators ?

Nonetheless, our goal is to obtain a useful tool for screening strains in the process of making new positive biofilms.

Since scores range from 0 to 1, on average, the trained model makes a 15.4% on the exclusion score. This corresponds to good performances as this type of model would be used to help choose between various strains, the strains predicgted to have good antagonistic properties destined to be more thoroughly tested in real experimental conditions afterwards.

Thus, this application does not require extreme precision, only a model that is able to accurately indicate when an antagonistic activity might be strong and when it might be weak. As said previously, a model that makes predictions of the actual antagonistic properties that are close enough to the truth would then be incredibly helpful for screening, which seems to be the case here.

## 5.2. Error analysis

To further study the errors made by the model and check the quality of predictions, we investigated error distributions and profiles.
<figure><img src="image_18.png" /><figcaption>**Fig. 8.** Error profile (upper figure) and Error by prediction category (lower figure). Error bars correspond to BCa 95% confidence interval around the mean. The **upper figure** shows the average proportion of several error ranges across all validation sets. The **lower figure** shows the average error the model makes when predicting a score in a specific range. **Note that we compare the predicted score to the error made, not the true score.** This enables us to evaluate which types of predictions the model is most trustworthy on. In this same figure, bars are colored based on the average proportion of the score category in validation sets. For example the color red for the left-most bar in the bottom plot shows that almost 50% of true exclusion score values are above 0.8 on average, indicating that in most test folds, around half of the samples show strong antagonism. This thus shows that the model is of quality (0.151 MAE) on a lot of points.</figcaption></figure>

### 5.2.1. Error distributions

In Fig. 8, we can see in the upper figure that 68% of predicted scores present an error of less than 0.2, and slightly less than half of predictions have an error of less than 0.1. Since almost 32% of predictions are wrong by a margin that exceeds 0.2, which can be considered as big errors, the apriori on model quality may be moderate.

However, this result is counterbalanced by the precision of the model when it predicts a strong antagonistic activity (the most interesting type for our biofilm experts). Indeed, the second graph of Fig. 8 reveals a very important result: when the model predicts a strong exclusion score (predicted score superior to 0.8) it makes an error of 0.151 on average, which is very satisfying for its application.

Since the main goal of the process of making positive biofilms is to identify which strains can have strong antagonistic effects on pathogens, having a model that can be this accurate when it predicts a strong interaction score is very valuable. Another important note is the fact that across all validation folds, almost half of the instances are scores with true scores above 0.8: Most samples show strong antagonism and the model is most accurate when predicting them. Moreover, this second plot shows that errors stay around 0.2 on average for less confident ranges which is reassuring in assessing model quality.

### 5.2.2. Stratified regression

Although results show that the model achieves satisfying performances, low and medium exclusion score ranges display less accurate

predictions. In an attempt to improve this, we applied stratified regression by stratifying according to exclusion score ranges. The goal was to determine whether or not a boost in performances for such ranges could be obtained by using range-specific submodels.

Results are discussed in Appendix E. They did not show an improvement in performances but rather a slight bias accumulation preventing the prediction accuracy from improving. This could be due to the fact that the stratified regressor used to make the final prediction depends on the output of a previously biased estimator (see protocol details in Appendix E).

Because stratification did not yield better results, the following sections are based on the non-stratified model.

### 5.2.3. Conformal prediction

Another aspect of model analysis corresponds to the study of aleatoric uncertainty at the sample level, i.e quantifying the uncertainty due to the noise in data [39]. In order to evaluate such uncertainty one approach consists of using conformal predictions (CP) [40,41]. In practice conformal prediction is often achieved using the Split/Inductive CP method [42] where in addition to the training and test sets, a conformalization set is created in order to compute conformality scores. These conformality scores allow the definition of intervals around a model predictions giving an estimation of aleatoric uncertainty.

To that end, we used the MAPIE Python package [43] and its implementation of Residual Normalized conformity scores which provide
<figure><img src="image_19.png" /><figcaption>**Fig. 9.** Average absolute error by Pathogen (first graph) and by *Bacillus* (second graph). Error bars BCa 95% confidence intervals. This reveals the heterogeneity in model performances across HCCV folds.</figcaption></figure>

adaptative conformal prediction intervals (see [44]). Results are shown in Appendix F. They suggest that the data is inherently quite noisy which could also explain why trained estimators cannot reach lower average errors easily.

### 5.2.4. Error analysis by organism and interactions

Given the heterogeneous structure of errors by prediction ranges and the extent of error bars in Fig. 7, results strongly suggested heterogeneous results across HCCV folds. Thus this led us to specifically look at errors by organisms and by interactions.

Fig. 9 shows a valuable insight: **the model is, overall, very good at predicting antagonism** for interactions involving completely new Pathogen or *Bacillus*. When the model ‘discovers’ a new organism in the test set of the HCCV protocol, the hardest organisms to predict are E. coli for pathogens and the 1234 *Bacillus* strain with respectively 0.186 and 0.247 MAE. Although the latter is quite high, these results reveal the **usefulness of our model for a variety of organisms** in a setting that is the closest to a real use case where potentially new strains or pathogens need to be assessed.

Moreover, when looking at model performances on totally new interactions Appendix G, the heatmap reveals even more heterogeneity across interactions. Notably, some interactions are extremely well described by general input feature patterns such as between the *Bacillus* strain 11457 and E. coli with an average absolute error of 0.079. On the other hand, we can also see an extremely badly predicted interaction: *Bacillus* strain 1234 and *S. enterica* with an MAE of 0.455. Results in Fig. 9 already pointed out the difficulty of predicting both organisms when removed individually, so it is expected that when removing both simultaneously, errors increase. Sur le fold ou la souche 1234 est

nouvelle et sur le fold où *S. enterica* est nouvelle, les erreurs étaient déjà parmi les plus élevées. La tâche de prédire deux organismes complètement nouveaux en intéraction étant plus dure, il n'est pas surprenant de voir que le fold ou l'interaction entre ces deux la est nouvelle est le plus mauvais.

### 5.2.5. Focus on the worst case

To push the analysis further, we investigated what distinguished the worst predicted interaction (1234 x *S. enterica*) from the best one (11457 x *E. coli*). To do that, the widely used SHAP values are very useful [26,27]. When computing such values, we must pay attention to the ordering of features (from most important to least important) and the relationship between feature values and contributions to prediction (SHAP value of a feature for a given observation). The beeswarm plot allows one to see both of the aforementioned simultaneously (see Fig. 10).

Before looking at the relationship between feature values and SHAP values, it is critical to verify the validity of the SHAP method. Indeed, they can be of paramount value but can also be misleading if used or analyzed poorly because of the theoretical assumptions they rely on, such as feature independence [26,45,46]. This extends to all explainability methods: each one is based upon a set of assumptions that may or may not be verified by data. Thus, a good approach is to compare different methods and extract from that comparison robust insights. Following our analysis of Permutation Feature Importances and Ablation study insights, we will compare the results obtained using SHAP values with the previous methods.

We can see that the order of features in both beeswarm plots is highly similar to that of PFI and Ablation study experiments (Modele and Pathogen Biofilm Height features are the first ones) which furthermore confirms previous insights.
<figure><img src="image_20.png" /><figcaption>**Fig. 10.** Beeswarm plots of model SHAP values for the worst (upper) and best (lower) predicted interactions, respectively (1234 x *S. enterica*) and (11457 x *E. coli*).</figcaption></figure>

By looking at the relationship between pathogen biofilm height values and their contribution to predictions, we can see a rather **heterogeneous pattern for the worst predicted interaction** whereas a **quite clear one is seen for the best predicted one**. In the latter, higher biofilm heights lead to weaker contributions revealing an interesting interaction mechanism. The fact that a highly heterogeneous pattern is observed for the worst interaction suggests that the mechanism between the two organisms is different and more intricate. Based on this and its performances, overall input features do not describe this interaction very well.

## 6. Discussion

Our study explored the use of machine learning, specifically the LGBMRegressor model, to predict microbial exclusion in biofilms based

on various biofilm properties, namely Mixing model for biofilm formation, Pathogen Biofilm Height, all *Bacillus* biofilm descriptors and engineered features (refer to Section 4.3.2 for details). Thanks to the HCCV protocol, model analysis revealed that our approach led to overall very satisfying results, although some interactions and organisms are poorly predicted.

It is important to note that our dataset was limited in genetic diversity, particularly among the *Bacillus* strains. This narrow scope limits the robustness of insights. Future studies would benefit from a larger, more genetically and phenotypically diverse dataset to capture a broader range of biofilm interactions. Additionally, while biofilm parameters such as volume, height, and roughness proved useful, they may not fully encapsulate the complexity of biofilm interactions. Adding features like EPS composition, metabolic activity, and molecular signaling could potentially enhance the model's ability to predict exclusion scores more accurately and homogeneously across interactions.

From a methodological perspective, exploring alternative models, such as neural networks or Bayesian frameworks, could be interesting on a substantially more genetically diverse dataset. supprimer le VPC par défautAssociation with wet lab experiments to further validate the usefulness of the approach is also a desired future direction.

Finally, while our machine learning approach is primarily data-driven, integrating biological knowledge (such as mechanistic models or ecological theories) could improve robustness. In several scientific fields, it is becoming more and more clear that combining data-driven and theory-driven approaches can often yield more accurate and insightful, suggesting that a hybrid approach could be beneficial here as well.

## 7. Conclusion

We developed a novel machine learning approach to study microbial interactions within biofilms using the LGBMRegressor algorithm and a specific protocol called HCCV. We achieved a weighted mean absolute error of 0.154 on a curated dataset of biofilm parameters. This approach provides valuable insights regarding key morphological features for prediction of antagonism, such as mixing condition and biofilm Height. It also presents high prediction accuracy when predicting strong exclusion scores, which could significantly help in selecting beneficial biofilm strains. Despite some limitations, such as dataset size and feature scope, our approach can contribute to the current knowledge and applications of biofilm research, by providing a new and powerful way of studying and predicting microbial interactions in biofilms, and by facilitating the selection and design of beneficial microbial communities for biofilm control.

### CRediT authorship contribution statement

**Raphaël Rubrice:** Writing – original draft, Visualization, Software, Resources, Methodology, Investigation, Formal analysis. **Virgile Guéneau:** Writing – review & editing, Resources, Data curation, Conceptualization. **Romain Briandet:** Writing – review & editing, Resources, Data curation, Conceptualization. **Antoine Cornuejols:** Writing – review & editing, Validation, Supervision, Methodology. **Vincent Guigue:** Writing – review & editing, Validation, Supervision, Methodology.
## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

This study was made possible thanks to the 'La Recherche Et Moi' program at AgroParisTech Engineering School. We thank all collaborators for their role in study conceptualization, manuscript writing, and editing. We also thank the MIMA2 platform (Microscopie et Imagerie des Microorganismes, Animaux et Aliments, https://doi.org/10.15454/1.5572348210007727E12) for microscopic observations. The French Agency for Research is acknowledged for the funding of the ANR LabCom project Biofilm1Health. Fig. 1 was designed with BioRender and Fig. 3 with Canva.

## Appendix A. Details on model selection

Because no statistical difference was observed between methods, we decided to base our choice on the best-performing model on average on the combinatoric dataset. By looking at top methods using weighted MAE and RMSE scores Fig. A.1, we can see that Microsoft's LGBMRegressor is selected.

## Appendix B. Details on preprocessing selection

To identify the best preprocessing techniques for imputation and normalization, we tested a variety of strategies. Results showed no effect of the choice of preprocessing, (see Fig. B.1). This can be explained by the nature of the LGBMRegressor. It is an implementation of Gradient Boosting Trees which are known to be robust by definition to varying ranges of input features and missing values. Indeed trees are threshold-based so no matter the scale of the descriptive variable chosen at a node, the result of thresholding operation is simply the same but scaled to another range of values once normalization is applied.

Regarding missing values, it can possibly mean two things: the features with missing values (i.e SubstratumCoverage and Roughness, see Fig. 2) are not very important for predicting exclusion scores so, no matter the imputed value, prediction is not affected, or, that diverse imputation strategies lead to similar imputation values which leads to little to no differences between them. The first hypothesis is confirmed by the results shown in Figs. 5, 6 and C.1 since those features have very low PFIs.

Using the mean would not make sense even though performances are similar because each strain or interaction is subject to different inner relationships. However, using the RandomForestImputer was considered overcomplicated with respect to the fact that no actual performance increase was observed. Thus, we chose the KNNImputer which we deemed to be a good compromise between conceptual validity and complexity. The RobustScaler was chosen due to its ability to mitigate the effects of outliers when using normalization which can be an interesting property when dealing with quite noisy data as is the case here.

<figure><img src="image_21.png" /><figcaption>**Fig. A.1.** Detailed results of best methods on the combinatoric dataset. Error bars are 95% Confidence intervals (CI) computed across all HCCV folds. All models were used after applying Mean Imputation and Standard Scaling to input variables.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Fig. B.1.** LGBMReessor RMSE scores using different preprocessing strategies. Error bars are the 95% Confidence Intervals. Colors show the average computed for a specific task (either *Bacillus* hold-outs, pathogen hold-outs, or interaction hold-outs).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. B.2.** LGBMReessor RMSE scores using the KNNImputer (Impute) or without (NoImpute). Error bars are the 95% Confidence Intervals.</figcaption></figure>
<table><caption>**Table D.1**

Hyperparameter ranges explored in the optuna study.</caption><thead><tr><th>Hyperparameter</th><th>Type</th><th>Range/Values</th><th>Additional Information</th></tr></thead><tbody><tr><td>n_estimators</td><td>Integer</td><td>100–1000 (step 10)</td><td>–</td></tr><tr><td>λ<sub>11</sub></td><td>Float</td><td>10<sup>-8</sup>–10.0</td><td>Logarithmic scale</td></tr><tr><td>λ<sub>12</sub></td><td>Float</td><td>10<sup>-8</sup>–10.0</td><td>Logarithmic scale</td></tr><tr><td>num_leaves</td><td>Integer</td><td>2–64</td><td>–</td></tr><tr><td>feature_fraction</td><td>Float</td><td>0.4–1.0</td><td>–</td></tr><tr><td>bagging_fraction</td><td>Float</td><td>0.4–1.0</td><td>–</td></tr><tr><td>bagging_freq</td><td>Integer</td><td>0–10</td><td>–</td></tr><tr><td>min_child Samples</td><td>Integer</td><td>5–100</td><td>–</td></tr><tr><td>metric</td><td>Categorical</td><td>{11, 12, huber, quantile}</td><td>–</td></tr></tbody></table>

<table><caption>**Table D.2**

Best trial results from the optuna study.</caption><thead><tr><th colspan="2">Best Trial Value</th></tr><tr><th>Metric</th><th>Value</th></tr></thead><tbody><tr><td>Cross Mean</td><td>0.1787</td></tr><tr><th colspan="2">Best Hyperparameters</th></tr><tr><td>n_estimators</td><td>160</td></tr><tr><td>λ<sub>11</sub></td><td>3.73546 × 10<sup>-8</sup></td></tr><tr><td>λ<sub>12</sub></td><td>8.64680 × 10<sup>-5</sup></td></tr><tr><td>num_leaves</td><td>23</td></tr><tr><td>feature_fraction</td><td>0.8953</td></tr><tr><td>bagging_fraction</td><td>0.8629</td></tr><tr><td>bagging_freq</td><td>3</td></tr><tr><td>min_child Samples</td><td>33</td></tr><tr><td>metric</td><td>huber</td></tr></tbody></table>

Because imputation can introduce significant bias, we also experimented without any kind of imputation. Identifying the imputation bias was done by comparing the performances with and without the KNNComputer. Results are shown in Fig. B.2 and show that we can indeed observe an imputation bias of 3.776e-03 with MAE scores of 0.150 and 0.153 obtained respectively for the approach without and with missing values replaced. This is however marginal and the 95% bootstrap confidence intervals show that both approaches lead to similar performances. Since the imputation bias is marginal, we decided to continue experimenting with imputation because in a real use case this allows much more flexibility.

## Appendix C. Details on feature engineering

To further improve performances, a model was trained on a dataset presenting an enriched feature space with additional Products between variables, Differences, Exponentiation, and two designed features (Density and FakeScore). Figs. C.1 and C.2 presents the heatmaps of PFIs for all features at each step.

## Appendix D. Details on hyperparameter search

Our hyperparameter search allowed us to slightly improve performances on average going from 0.1804 to 0.1787 (these values correspond to cross metric evaluation, i.e. the average of weighted MAE and RMSE). Details regarding the exploration of the hyperparameter space are provided in Table D.1 and the optimal hyperparameters are specified in Table D.2.

## Appendix E. Details on stratified regression

Startified regression was achieved by separating exclusion scores in 5 classes (low, medium low, medium high, high and very high) using the following protocol: (1) The training set is divided in 6 subsets (1 to train a submodel on all ranges of exclusion scores, called the oracle, and 5 range-specific training sets). (2) The oracle is trained to predict

<figure><img src="image_24.png" /><figcaption>**Fig. C.1.** Detailed results of step 1 of feature selection via Permutation Feature Importances on the feature engineered dataset.</figcaption></figure>

exclusion scores (3) Range-specific submodels are trained on samples falling into the specific range according to the oracle's predictions. (This is done in order to be able to use the approach at inference time when the true score is not known).

Two types of range-specific training protocols were explored. The first one consists of simply using the true scores while, in the second, the submodel learns to predict both the true scores and the scores predicted by the oracle. The latter was attempted in order to try to make the range-specific submodels closer to the oracle in hope of mitigating bias propagation. Both protocols are respectively referred to as 'Default' and 'Mixed' range-specific training. An illustration of the stratified regression protocols is shown in Fig. E.1.
<figure><img src="image_25.png" /><figcaption>Fig. C.2. Detailed results of steps 2 and 3 of feature selection via Permutation Feature Importances on the feature engineered dataset.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Fig. E.1. Illustration of the stratified regression protocol. Oracle refers to the model that is learned on a non-altered distribution of exclusion scores and which is used at inference time to redirect samples to range-specific models.</figcaption></figure>

Exclusion score ranges were defined either using the 20th, 40th, 60th and 80th distribution quantiles or by using the domain-knowledged informed threshold values of 0.2, 0.4, 0.6 and 0.8. These approaches are respectively referred to as ‘Quantile’ or ‘Custom’ range specifications, in both of them, 5 classes are defined.

Comparison of the different stratification protocols with one another as well as with respect to the non-stratified approach are presented in Fig. E.2. Results show that stratification does not improve performances compared to its non stratified counterpart. When considering imputation, the best approach remains the one presented in the main body of this research article with a Weighted MAE of 0.153 (the slight difference from 0.154 stated earlier was observed after a changed of hardware on which the stratified experiments were conducted). The lack of improvement could be attributed to the accumulation of bias at inference time: the oracle plays a key role but is already biased, the submodels are also biased thus if the benefit of using the specified model is not substantial, bias accumulation stands out.

Nonetheless, it is worth noting that the imputation bias remains similar even when using stratification. Another note is that the ‘Mixed’ training approach does indeed help to reduce bias propagation since all experiments using it led to better results on average than the ‘Default’ counterparts. Moreover, using carefully designed thresholds to define ranges based on biofilm expert knowledge (‘Custom’ thresholds) yields slightly better performances than simply relying on quantiles. This underscores the importance of domain experts in the development of machine learning models.

## Appendix F. Details on uncertainty quantification

In order to also investigate aleatoric uncertainty we analyzed average coverage and conformal prediction interval widths for the main model presented in this research article but also for the model without imputation and the best performing stratified model. Coverage

indicates the proportion of instances where the conformal interval contained the true value. Interval width can be useful to interpret coverage results as larger intervals are correlated with larger coverage. Results are presented in Fig. F.1 and Fig. F.2. We can observe that coverage is above 0.7 with the best value obtained by the model with imputation at 0.762. The latter displays a better coverage than the stratified model while having smaller interval width (respectively 0.805 and 0.823). Nonetheless, these interval widths are massive when put in perspective of the range of exclusion scores. This indicates that the aleatoric uncertainty remains quite substantial in our data which could further explain why trained models do not display better accuracy.

## Appendix G. Details on model performances across prediction tasks

Thanks to the HCCV protocol, we can confidently assess the ability of our approach to be useful in real use cases since each fold is designed to confront the model on partially new interactions (One organism involved in the interaction was removed from training data) or totally new ones (both organisms involved were removed from training data), see Fig. 3. Fig. 9 presents the insight brought by looking at the former and below are presented the results for the latter (see Fig. G.1).

## Data availability

The data and code used for this paper are available at https://github.com/raphaelrubrice/BiofilmAntagonismPrediction. The broader dataset used is provided by (Guéneau et al., 2024) and available at https://entrepot.recherche.data.gouv.fr/dataset.xhtml?persistentId=doi:10.57745/XRXQEI.
<figure><img src="image_27.png" /><figcaption>**Fig. E.2.** Performances using stratified regression using diverse protocols. Imputation status is indicated by “Impute” or “NoImpute”, the method for defining exclusion score ranges is either “Quantile” (20th, 40th, 60th, 80th) or “Custom” (0.2, 0.4, 0.6, 0.8) and the training protocol of submodels is specified in each row by “Default” or “Mixed”. “Normal” indicates that the estimator is not stratified.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Fig. F.1. Model conformal prediction interval coverage (using a target coverage of 90%).</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. F.2. Model conformal prediction interval widths.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Fig. G.1. Model performances on completely new interactions.</figcaption></figure>

## References

[1] Sauer K, Stoodley P, Goeres DM, Hall-Stoodley L, Burmølle M, Stewart PS, et al. The biofilm life cycle: expanding the conceptual model of biofilm formation. Nat Rev Microbiol 2022;20(10):608-20.

[2] Flemming HC, Wuertz S. Bacteria and archaea on earth and their abundance in biofilms. Nat Rev Microbiol 2019;17(4):247-60.

[3] Nadell CD, Drescher K, Foster KR. Spatial structure, cooperation and competition in biofilms. Nat Rev Microbiol 2016;14(9):589-600.

[4] Guéneau V, Plateau-Gonthier J, Arnaud L, Piard JC, Castex M, Briandet R. Positive biofilms to guide surface microbial ecology in livestock buildings. Biofilm 2022;4:100075.

[5] Gueneau V, Guillier L, Berdous C, Noirot-Gros M-F, Jimenez G, Plateau J, et al. Guided assembly of multispecies positive biofilms targeting undesirable bacteria. 2024, BioRxiv, 2024-10.
[6] Guéneau V, Charron R, Costache V, Bridier A, Briandet R. Spatial analysis of multispecies bacterial biofilms. In: Methods in microbiology, vol. 53, Elsevier; 2023, p. 275-307.

[7] Dzianach PA, Dykes GA, Strachan NJ, Forbes KJ, Pérez-Reche FJ. Challenges of biofilm control and utilization: lessons from mathematical modelling. J R Soc Interface 2019;16(155):20190042.

[8] Greener JG, Kandathil SM, Moffat L, Jones DT. A guide to machine learning for biologists. Nature Rev Mol Cell Biol 2022;23(1):40-55.

[9] Nestor E, Toledano G, Friedman J. Interactions between culturable bacteria are predicted by individual species' growth. Msystems 2023;8(2):e00836-22.

[10] Emmenegger B, Massoni J, Pestalozzi CM, Bortfeld-Miller M, Maier BA, Vorholt JA. Identifying microbiota community patterns important for plant protection using synthetic communities and machine learning. Nat Commun 2023;14(1):7983.

[11] Barton DL, Chang YR, Ducker W, Dobnikar J. Data-driven modelling makes quantitative predictions regarding bacteria surface motility. PLoS Comput Biol 2024;20(5):e1012063.

[12] Dergham Y, Le Coq D, Bridier A, Sanchez-Vizuete P, Jbara H, Deschamps J, et al. Bacillus subtilis ndmed, a model strain for biofilm genetic studies. Biofilm 2023;100152.

[13] Li Z, Wang X, Wang J, Yuan X, Jiang X, Wang Y, et al. Bacterial biofilms as platforms engineered for diverse applications. Biotech Adv 2022;57:107932.

[14] Sharma S, Sharma S, Tiwari V. Role of biofilm in host-pathogen interaction. In: A complete guidebook on biofilm study. Elsevier; 2022, p. 227-72.

[15] Zhou J, Cai Zh. Microbial social interactions in biofilm. Implicit Quor Sens Syst Biofilm Form Virulence 2018;29-46.

[16] Molloy S. Biofilm microanatomy. Nat Rev Microbiol 2013;11(5). 301-301.

[17] Hartmann R, Jeckel H, Jelli E, Singh PK, Vaidya S, Bayer M, et al. Quantitative image analysis of microbial communities with BiofilmQ. Nat Microbiol 2021;6(2):151-6.

[18] Grinsztajn L, Oyallon E, Varoquaux G. Why do tree-based models still outperform deep learning on typical tabular data? Adv Neural Inf Process Syst 2022;35:507-20.

[19] Guéneau V, Berdous C, Noirot-Gros M-F, Serror P, Castex M, Briandet R. Confocal laser scanning microscopy (CLSM) images to study microbial competition in multi-strains biofilms. Recherche Data Gouv; 2024, http://dx.doi.org/10.57745/XRXQEI.

[20] Guéneau V. Vers l'utilisation rationnelle de biofilms positifs pour aider à maîtriser les microorganismes indésirables dans les bâtiments d'élevage [Ph.D. thesis], Université Paris-Saclay; 2023.

[21] Alkharusi H. Categorical variables in regression analysis: A comparison of dummy and effect coding. Int J Educ 2012;4(2):202.

[22] Asnicar F, Thomas AM, Passerini A, Waldron L, Segata N. Machine learning for microbiologists. Nat Rev Microbiol 2023;1-15.

[23] Efron B. Better bootstrap confidence intervals. J Amer Statist Assoc 1987;82(397):171-85.

[24] Sheppard K, Khrapov S, et al. basitage/arch: Release 7.2. Zenodo; 2024, http://dx.doi.org/10.5281/zenodo.14035889.

[25] Shapley LS, et al. A value for n-person games. Princeton University Press Princeton; 1953.

[26] Molnar C. Interpretable machine learning. 3rd ed.. 2025, URL: https://christophm.github.io/interpretable-ml-book.

[27] Lundberg SM, Lee SI. A unified approach to interpreting model predictions. In: Guyon I, Luxburg UV, Bengio S, Wallach H, Fergus R, Vishwanathan S, et al., editors. Advances in neural information processing systems 30. Curran Associates, Inc.; 2017, p. 4765-74, URL: http://papers.nips.cc/paper/7062-a-unified-approach-to-interpreting-model-predictions.pdf.

[28] Boser BE, Guyon IM, Vapnik VN. A training algorithm for optimal margin classifiers. In: Proceedings of the fifth annual workshop on computational learning theory. 1992, p. 144-52.

[29] Friedman JH. Greedy function approximation: a gradient boosting machine. Ann Stat 2001;1189-232.

[30] Breiman L. Random forests. Mach Learn 2001;45:5-32.

[31] Hoerl AE, Kennard RW. Ridge regression: applications to nonorthogonal problems. Technometrics 1970;12(1):69-82.

[32] Tibshirani R. Regression shrinkage and selection via the lasso. J R Stat Soc Ser B Stat Methodol 1996;58(1):267-88.

[33] Zou H, Hastie T. Regularization and variable selection via the elastic net. J R Stat Soc Ser B Stat Methodol 2005;67(2):301-20.
[34] Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, et al. Scikit-learn: Machine learning in python. J Mach Learn Res 2011;12:2825-30.

[35] Chen T, Guestrin C. Xgboost: A scalable tree boosting system. In: Proceedings of the 22nd ACM Sigkdd international conference on knowledge discovery and data mining. 2016, p. 785-94.

[36] Ke G, Meng Q, Finley T, Wang T, Chen W, Ma W, et al. Lightgbm: A highly efficient gradient boosting decision tree. Adv Neural Inf Process Syst 2017;30.

[37] Raschka S, Patterson J, Nolet C. Machine learning in python: Main developments and technology trends in data science, machine learning, and artificial intelligence. 2020, arXiv preprint arXiv:2004.04803.

[38] Fisher A, Rudin C, Dominici F. All models are wrong, but many are useful: Learning a variable's importance by studying an entire class of prediction models simultaneously. J Mach Learn Res 2019;20(177):1-81.

[39] Hüllermeier E, Waegeman W. Aleatoric and epistemic uncertainty in machine learning: An introduction to concepts and methods. Mach Learn 2021;110(3):457-506.

[40] Saunders C, Gammerman A, Vovk V. Transduction with confidence and credibility. 1999.

[41] Vovk V, Gammerman A, Shafer G. In: Algorithmic learning in a random world, vol. 29, Springer; 2005.

[42] Papadopoulos H, Proedrou K, Vovk V, Gammerman A. Inductive confidence machines for regression. In: Machine learning: ECML 2002: 13th European conference on machine learning helsinki, Finland, August 19-23, 2002 proceedings 13. Springer; 2002, p. 345-56.

[43] Cordier T, Blot V, Lacombe L, Morzadec T, Capitaine A, Brunel N. Flexible and Systematic Uncertainty Estimation with Conformal Prediction via the MAPIE library. In: Conformal and probabilistic prediction with applications. 2023.

[44] Lei J, G'Sell M, Rinaldo A, Tibshirani RJ, Wasserman L. Distribution-free predictive inference for regression. J Amer Statist Assoc 2018;113(523):1094-111.

[45] Fryer D, Strümke I, Nguyen H. Shapley values for feature selection: The good, the bad, and the axioms. Ieee Access 2021;9:144352-60.

[46] Aas K, Jullum M, Løland A. Explaining individual predictions when features are dependent: More accurate approximations to Shapley values. Artificial Intelligence 2021;298:103502.