# Revealing the dynamics and mechanisms of bacterial interactions in cheese production with metabolic modelling

Maxime Lecomte, Wenfan Cao, Julie J. Aubert, David James Sherman, Hélène Falentin, Clémence Frioux, Simon Labarthe

► To cite this version:

Maxime Lecomte, Wenfan Cao, Julie J. Aubert, David James Sherman, Hélène Falentin, et al. Revealing the dynamics and mechanisms of bacterial interactions in cheese production with metabolic modelling. 2023. hal-04088301v2

HAL Id: hal-04088301

https://hal.inrae.fr/hal-04088301v2

Preprint submitted on 22 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Revealing the dynamics and mechanisms of
bacterial interactions in cheese production
with metabolic modelling

Maxime Lecomte¹,², Wenfan Cao¹, Julie Aubert³, David James
Sherman², Hélène Falentin*¹, Clémence Frioux*², and Simon
Labarthe*⁴,²

¹ *Univ. Rennes, INRAE, STLO, F-35042 Rennes, France*

² *Inria, Univ. Bordeaux, INRAE, F-33400 Talence, France*

⁴ *Univ. Bordeaux, INRAE, BIOGECO, Cestas, France*

³ *Univ. Paris-Saclay, AgroParisTech, INRAE, UMR MIA Paris-Saclay, 91120, Palaiseau, France*

* Co-last authors

Corresponding Author: Simon Labarthe - simon.labarthe@inrae.fr

CC-BY 4.0
Abstract

Cheese taste and flavour properties result from complex metabolic processes occurring in microbial communities. A deeper understanding of such mechanisms makes it possible to improve both industrial production processes and end-product quality through the design of microbial consortia. In this work, we caracterise the metabolism of a three-species community consisting of *Lactococcus lactis*, *Lactobacillus plantarum* and *Propionibacterium freudenreichii* during a seven-week cheese production process. Using genome-scale metabolic models and omics data integration, we modeled and calibrated individual dynamics using monoculture experiments, and coupled these models to capture the metabolism of the community. This model accurately predicts the dynamics of the community, enlightening the contribution of each microbial species to organoleptic compound production. Further metabolic exploration revealed additional possible interactions between the bacterial species. This work provides a methodological framework for the prediction of community-wide metabolism and highlights the added value of dynamic metabolic modeling for the comprehension of fermented food processes.

**Keywords** Systems biology - Fermentation - Metabolic modeling - Microbial community - Cheese - Dynamics - Flux Balance Analysis
# 1 Introduction

Fermentation of food and beverages by microbial consortia is one of the oldest food processing technologies, associated to health benefits, and used today in both traditional and industrial contexts (Tamang et al., 2016b,a). The transformation of raw food material by micro-organisms is the result of metabolic processes and interactions between the members of the microbial community. The content of the latter can be controlled, such as in industrial starters (Somerville et al., 2021), for food safety and for the predictability of organoleptic properties (Galimberti et al., 2021), i.e. pertaining to the taste and flavour. While functional properties of microbial strains associated to fermented food are better and better understood (Tamang et al., 2016a), the precise molecular mechanisms at stake are only beginning to be elucidated (Blasche et al., 2021). In particular, the understanding of the individual metabolism of consortia members is not sufficient to capture or predict the global behaviour of the community system, which is impacted by microbe-microbe interactions.

Systems biology approaches combined with omics data acquisition are a promising avenue for characterising the functional landscape and interactions occurring in microbial consortia of fermented foods (Mannaa et al., 2021). Working with empirical communities of unknown composition motivates culture-independent inventories of metabolite and species as well as cultures of isolates (Blasche et al., 2021). For controlled community starters,
a comparison between mono- and co-cultures can inform on the impact of
interactions on the dynamics and metabolic profiles of the community Öz-
can et al. (2020). The use of genome-scale metabolic models (GEMs) (Fang
et al., 2020) seems particularly opportune as a means to associate the strains
of the consortia to the functions harboured by the systems and ultimately to
the organoleptic properties of the food (Özcan et al., 2020). Such models of
microbial consortia can be useful for numerous applications in food biotech-
nology, from the assembly of strains in starter design, to the optimization of
processes for desired organoleptic properties (Rau and Zeidan, 2018).

The nature of the food matrix, the complexity of the bacterial commu-
nity, and the highly strain-specific secondary metabolic pathways involved in
the organoleptic features make it difficult to built mechanistic models of the
biological system. Therefore, working with small-scale, controlled systems
of reduced complexity is particularly relevant (Rau and Zeidan, 2018). Yet
even simple systems must be considered and modeled over time, as fermen-
tation is an intrinsically dynamic process. In this work, we are interested
in modeling the metabolism of cheese making, from milk to ripening using
a bacterial community starter mimicking industrial cheese production, com-
posed of two Lactic Acid Bacteria (LAB), *Lactococcus lactis* and *Lactiplan-
tibacillus plantarum*, and one propionic acid bacterium, *Propionibacterium
freudenreichii*. Previous work (Cao et al., 2021) described how a fine-tuning
of cheese standardised process parameters could modulate the production of
aroma compounds. Here, using omics data, we build and exploit metabolic
networks and community models in order (i) to evaluate the metabolic con-
tribution of each microbial species during cheese making, (ii) to simulate the
dynamics of the microbial metabolisms and the effect of the community over
the individual behaviours, and (iii) propose possible interactions among the
bacterial species.

Our work identifies the metabolic pathways responsible for nutrient trans-
formation in milk. We optimized our individual models of metabolism and
were able to reproduce the dynamics of pure cultures. From these models,
we constructed a dynamic model of the bacterial community that accurately
predicts the production of aroma compounds in the cheese.

2 Results

## 2.1 Individual metabolic models identify specific metabolic pathways in a milk environment

The microbial community starter was composed of *Propionibacterium freudenreichii* CIRM-BIA122, *Lactiplantibacillus plantarum* CIRM-BIA465 and *Lactococcus lactis* biovar *diacetylactis*. These species have been selected for their metabolic complementarity, their ability to produce organoleptic compounds and the extended literature available regarding their metabolism (Cao et al., 2021). *L. lactis* and *P. freudenreichii* are common species in cheese. *L. plantarum* is an ubiquitous species found in cheese and in a large variety of fer-
mented food of vegetal origin. In addition, there were existing metabolic
networks available for strains of *L. lactis* (Flahaut et al., 2013) and *L. plantarum* (Teusink et al., 2006), providing support for model refinement.

Full-length total DNA was obtained for each strain. Assembling bacte-
rial genomes from a combination of short and long reads produced accurate
and complete assemblies revealing one circular chromosome and the presence
of 7 plasmids for the *L. lactis* strain and 4 plasmids for the *L. plantarum*
strain. The *P. freudenreichii* strain did not harbour any plasmids. Anno-
tated genome can be retrieved at ENA under accession ENA:PRJEB42478.
The number of protein encoding genes, ribosomal sequences, tRNA and GC
percent and accession numbers are gathered in Table 1. Annotation identified
a list of 8234 coding sequences including those from chromosomal (2460 for
*P. freudenreichii*, 2417 for *L. lactis* and 3016 for *L. plantarum*) and plasmid
origins (139 for *L. lactis* and 82 for *L. plantarum*).

Genome-scale metabolic networks were reconstructed for each individ-
ual strain starting from their respective annotated genomic sequences. Af-
ter an automatic first step of reconstruction, we performed manual refine-
ments to the central carbon metabolism based on literature and biological
knowledge, in order to account for each strain's specificity. The resulting
metabolic networks satisfy the quality expectations in the field: a low frac-
tion of universally-blocked reactions (< 3%), and a low fraction of reactions
without gene-proteins-reactions rules (< 30%) (Lieven et al., 2020). Table 2
describes the characteristics of the three models. We observe that *L. lactis*
has a smaller number of reactions than the other two strains, a qualitative
difference that is also observed for other strains of the same species in other
databases (Noronha et al., 2018).

We built genome-scale metabolic models (GEMs) from the three metabolic
reconstructions using Flux Balance Analysis (FBA) (Palsson and Varma,
1994; Orth et al., 2010), ensuring the production of biomass nutrients on a
nutritional environment reflecting milk composition (see Materials and Meth-
ods). We obtained growth rates ranging from 0.187 mmol.gDW⁻¹.hr⁻¹ to
2.17 mmol.gDW⁻¹.hr⁻¹, the former obtained for *P. freudenreichii* and the
highest value obtained for *L. lactis*. Particular attention was given to the
role of oxygen in the three metabolisms. Since bacterial growth on cheese
corresponds to microaerobic conditions (Miyoshi et al., 2003), we altered
the upper bounds of the oxygen transport reactions to limit their cytosolic
import.

The central carbon metabolism produces NADH and NADPH reducing
equivalents, ATP and precursors needed for the biosynthesis of molecules es-
sential for the growth and the metabolic activities of bacterial cells e.g. amino
acids, purine, pyrimidine, glycerol-3-phosphate, fatty acids, N-acetyl-glucosamine,
vitamins. The three studied species shared almost the same glycolysis and
pentose phosphate pathways. Carbon sources, in this case lactose, citric and
lactic acids from milk, are converted into pyruvic acid, which is then metab-
<table><caption>Table 1: Genomic overview of *P. freudenreichii*, *L. lactis* and *L. plantarum* after assembly.</caption><tr><td>Bioproject Number</td><td colspan="3">PRJEB42478</td></tr><tr><td>Genome name</td><td>Propionibacterium freudenreichii</td><td>Lactiplantibacillus plantarum</td><td>Lactococcus lactis</td></tr><tr><td>Strain name</td><td>CIRM-BIA122</td><td>CIRM-BIA465</td><td>lactis biovar diacetyllactis CIRM-BIA1206</td></tr><tr><td>Taxon id</td><td>1744</td><td>1590</td><td>1358</td></tr><tr><td>Accession number</td><td>ERS5564379</td><td>ERS5564378</td><td>ERS5564377</td></tr><tr><td>number of plasmids</td><td>Ø</td><td>4</td><td>7</td></tr><tr><td>size of plasmid (bases)</td><td>Ø</td><td>p1: 40748, p2: 30463, p3: 9152, p4: 2012</td><td>p1: 40748, p2: 30463, p3: 9152, p4: 2012, p5: 30463, p6: 9152, p7: 2012</td></tr><tr><td>Size of the circular genome (bases)</td><td>2622405</td><td>3121980</td><td>2365039</td></tr><tr><td>number of rRNA</td><td>6</td><td>16</td><td>19</td></tr><tr><td>number of tRNA</td><td>45</td><td>67</td><td>62</td></tr><tr><td>percentage of GC (genome, min & max in plasmids)</td><td>67</td><td>44 (36-41)</td><td>35 (30-36)</td></tr><tr><td>number of CDS</td><td>2472</td><td>3112</td><td>2563</td></tr><tr><td>number of genes in each replicon</td><td>chromosome: 2557</td><td>p1: 41, p2: 27, p3: 13, p4: 2<br>chromosome: 3157</td><td>p1: 69, p2: 43, p3: 14, p4: 9, p5: 5, p6: 4, p7: 3, chromosome: 2542</td></tr></table>
<table><caption>Table 2: Characteristics of the three genome-scale metabolic network reconstructions and their associated FBA models in milk medium at the inoculation point.</caption><thead><tr><th></th><th><i>P. freudenreichii</i></th><th><i>L. plantarum</i></th><th><i>L. lactis</i></th></tr></thead><tbody><tr><td>Number of genes</td><td>1473</td><td>1433</td><td>1272</td></tr><tr><td>Number of metabolites</td><td>1284</td><td>1045</td><td>939</td></tr><tr><td>Number of reactions</td><td>1790</td><td>1523</td><td>1337</td></tr><tr><td>Percentage of reactions associated to genes</td><td>76.6</td><td>79.0</td><td>82.1</td></tr><tr><td>Growth rate in mmol.gDW<sup>-1</sup>.hr<sup>-1</sup></td><td>0.187</td><td>0.645</td><td>2.172</td></tr></tbody></table>

olized *via* specific pathways. The specificity of each strain's central carbon metabolism is illustrated in Figure 1. We describe below the characteristics of each strain together with the main refinements performed to their GEMs, paying specific attention to secondary pathways essential for organoleptic compound production.
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Metabolic maps of specific and shared pathways for *L. lactis* (a), *L. plantarum* (b) and *P. freudenreichii* (c). The glycolysis pathway (grey) is shared by the three species, while Leloir (salmon), acetolactate (pink) and citrate degradation pathways are shared by the two lactic acid bacteria. The TCA cycle (light blue) is shared by *P. freudenreichii* and *L. lactis*. Tagatose (green), transketolase (blue) and Wood-Werkman (purple) pathways are specific to *L. lactis*, *L. plantarum* and *P. freudenreichii* respectively. Orange and green metabolites are respectively inputs and outputs of the metabolic models. *Compound Abbreviations:* 13dpg: 3-Phospho-D-glyceroyl phosphate; 6pgc: 6-Phospho-D-gluconate; 6pgl: 6-phospho-D-glucono-1,5-lactone; ac: acetate; acald: acetaldehyde; accoa: acetyl-coa; actn_R: acetoine; akg: 2-Oxoglutarate; alac_S: (S)-2-Acetolactate; btd_RR: butanediol; cit: citrate; dgal6p: d-galactose-6-phosphate; dhap: Dihydroxyacetone phosphate; diact: diacetyl; f6p: fructose-6-phosphate; fdp: D-Fructose 1,6-bisphosphate; fum: fumarate; g1p: glucose-1-phosphate; g3p: Glyceraldehyde 3-phosphate; g6p: glucose-6-phosphate; gal: galactose; gal1p: galactose-1-phosphate; glc_D: glucose; lac_D, lac_L: lactate; lac6p: lactose-6-phosphate; lcts: lactose; mmcoa_R: (R)-Methylmalonyl-CoA; mmcoa_S: (S)-Methylmalonyl-CoA; oaa: Oxaloacetate; ppa: propionate.; ppap: Propanoyl phosphate; ppcoa: propionyl-coa; pyr: pyruvate; ru5p_D: D-Ribulose 5-phosphate; succ: succinate; succoa: succinyl-coa; tag6p_D: D-tagatose-6-phosphate; tagdp_D: D-Tagatose 1,6-biphosphate; udpg: UDP-glucose; xu5p_D: D-Xylulose 5-phosphate. *Reaction Abbreviations:* 2131pyrpp: Methylmalonyl-CoA carboxyltransferase 5S subunit; 6PGALSZ: 6-phospho-beta-galactosidase; ACLD: Acetolactate decarboxylase; ACTD2: Acetoin dehydrogenase; ACLDC: Acetolactate decarboxylase; ACLS: Acetolactate synthase; BTDD_RR: R R butanediol dehydrogenase; CITL: Citrate lyase; CITt4_1: Citrate transport via sodium symport; CS: Citrate synthase; D_LACt2: D lactate transport via proton symport; ENO: Enolase; FBA: Fructose-bisphosphate aldolase; FEDCabc: FEDCabc; FRD2rpp: Fumarate reductase / succinate dehydrogenase (irreversible) (periplasmic, membrane potential dissipating); G6PDH2r: Glucose 6-phosphate dehydrogenase; GAL6PI: Galactose-6-phosphate isomerase; GALkr: Galactokinase; GALUi: UTP-glucose-1-phosphate uridylyltransferase (irreversible); GAPD: Glyceraldehyde-3-phosphate dehydrogenase; GND: Phosphogluconate dehydrogenase; HEX1: Hexokinase (D-glucose:ATP); L_LACt2r: L lactate reversible transport via proton symport; LACpts: Lactose transport via PEP:Pyr PTS; LACZ: B-galactosidase; LCTSt3ipp: Lactose transport via proton aniport (periplasm); LDH_D: D-lactate dehydrogenase; LDH_L: L-lactate dehydrogenase; MME: Methylmalonyl-CoA epimerase; MMM2: Methylmalonyl-CoA mutase; OOR3r: 2-oxoglutarate synthase (rev); PC: Pyruvate carboxylase; PFK: Phosphofructokinase; PFK_2: Phosphofructokinase; PFL: Pyruvate formate lyase; PGI: Glucose-6-phosphate isomerase; PGK: Phosphoglycerate kinase; PGL: 6-phosphogluconolactonase; PGM: Phosphoglycerate mutase; PGMT: Phosphoglucomutase; phosphoketolase: phosphoketolase; PPAKr: Propionate kinase; PPCSCT: Propanoyl-CoA: succinate CoA-transferase; PTA: Phosphotransacetylase; PTA2: Phosphate acetyltransferase; PYK: Pyruvate kinase; RPE: Ribulose 5-phosphate 3-epimerase; TGBPA: Tagatose-bisphosphate aldolase; UGLT: UDPglcD: hexose-1-phosphate uridylyltransferase</figcaption></figure>
2.1.1 LAB metabolism

Both *L. lactis* and *L. plantarum* utilize lactose and citric acid, the two main
carbon sources in milk (Widyastuti et al., 2014). Lactose is degraded by
a beta-galactosidase into glucose and galactose. Glucose feeds glycolysis,
while galactose and citric acid metabolisms (Palles et al., 1998) are strain-
dependent. Galactose is converted into glucose-6-phosphate *via* Leloir path-
ways before entering in glycolysis. Citric acid is converted into pyruvate by
the citrate lyase multimeric enzyme encoded by the *citDEF* operon. Both
species lead to the acidification of milk through the production of lactic and
acetic acid (de Roissart, 1994). The production of butanediol and diacetyl,
responsible for a buttery flavour, occurs in both models through the aceto-
lactate pathway from the transformation of pyruvate. All of these pathways
are present and functional in the model (Fig. 1 a, b).

**L. lactis** The production of butanediol required curation of the GEM (Fig. 1 a). The acetoin-dehydrogenase flux was blocked (*ACTD2*) and bounds of the acetolactate decarboxylase (*ACLDC*), acetolactate synthase (*ACLS*) were refined permitting the activation of the acetolactate pathway and the production of butanediol, in accordance with the literature (Carroll et al., 1999; Swindell et al., 1996; Makhlouf, 2006). Finally, the consumption of lactose was regulated by refining the flux bounds of *LACZ* reaction which enabled a lactate production flux for *L. lactis* (Fig. 1 a).

To validate this model, we considered an example of pathway activation
during cheese production. *L. lactis* consumes lactose either through the tagatose pathway, encoded in the *lacABCD* plasmidic operon, or by the Leloir pathway, encoded in the *galKTE* chromosomal operon. Metatranscriptomic data shows that *L. lactis* expresses the genes of the tagatose pathway at the beginning of cheese production, when lactose concentration is high, and those of the Leloir pathway during ripening, when lactose concentration is low. This strongly suggests that lactose concentration determines which pathway is employed, and predicts that a higher availability of lactose would produce a higher lactose flux through the tagatose pathway. The FBA model of *L. lactis* indeed leads to a greater flow of lactose through the tagatose pathway at the inoculation, and therefore meets the expectations (Fig. 2 a).

*L. plantarum* Glucose and galactose degradation lead to pyruvic acid, which contributes to the release of lactic acid only in a homolactic metabolism or ethanol, acetic and lactic acids in a heterolactic metabolism. As *L. plantarum* is a heterolactic bacteria, it involves the specific transketolase enzyme (EC 4.1.2.9), previously described in (Abedi and Hashemi, 2020); thus, acetic acid can be produced via the acetate kinase pathway. The production of butanediol required similar curation as for *L. lactis*. The conversion of pyruvate into acetic acid was ensured by blocking the pyruvate formate lyase (PFL) (Quatravaux et al., 2006). To reproduce the degradation of galactose and glucose through Leloir and glycolysis pathways by the β-galactosidase LACZ, the *LCTSt3ipp* reaction was knocked out. To complete the galactose
degradation, the UTP-glucose-1-phosphate uridylyltransferase (*GALUi*) was
made reversible.

2.1.2 *P. freudenreichii* metabolism

The CIRM-BIA122 strain of *P. freudenreichii* consumes either lactic acid, its main source of carbon (Thierry et al., 2011; Dank et al., 2021) or lactose (Loux et al., 2015). Lactic acid and lactose are catabolized through the glycolysis and the pentose-phosphate pathways, leading to pyruvic acid. One part of pyruvic acid fuels the tricarboxylic acid cycle (TCA) and the Wood-Werkman pathway (Fig. 1 c). Production of propionic acid, specific to the *Propionibacterium* genus, is ensured by routing succinic acid from TCA. The production of acetic acid and CO2 (Turgay et al., 2020) from pyruvate is permitted through the activity of two enzymes: the pyruvate dehydrogenase and the decarboxylative oxidation.

Ensuring the correct propionate-to-acetate ratio required further curation
of the GEM. The bounds of *PPAKr*, *PTAr*, *2131pyrpp* reactions were refined
in order to drive propionate through the predominant production pathway,
i.e. Wood-Werkman, and *POX2* reaction was blocked to regulate flux of
acetate.

To validate the model, we computed the ratio of propionate/acetate
and obtained a value of 2.19, consistent with biological observations of 1:2
(Vaughan L. Crow, 1986; Dank et al., 2021; Cao et al., 2021; Turgay et al.,
2020).
After refinement, FBA models of the three species predicted the pro-
duction of several aroma compounds such as lactic, acetic, propionic acids,
diacetyl, and butanediol (Smid and Kleerebezem, 2014; Cao et al., 2021)
(Fig. 1). They furthermore accurately reproduce the metabolic specificities
of each strain with respect to their primary metabolism, i.e. Wood-Werkman
pathway for *P. freudenreichii*, tagatose pathway for *L. lactis* and transketo-
lase pathway for *L. plantarum*.

## 2.2 Individual dynamic model optimization

Each microbial strain was cultivated in pure culture: monitoring of pH for
LAB, and growth, for all strains, was performed in milk under microaerophilic
conditions (Supp. Table S1). To mimic lactate availability in co-culture,
the medium was supplemented with lactate and a hydrolyzate of milk pro-
teins simulating proteolysis for *P. freudenreichii*. Under these conditions,
all strains reached a culturability above 9 log10 CFU g⁻¹ and pH initially
at 6.7 at inoculation time reached 5.7 and 5.1 in *L. plantarum* and *L. lactis*
cultures respectively. The plateau of culturability corresponding to the entry
into stationary phase was reached in 48 hours for *P. freudenreichii*, 14 hours
for *L. plantarum* and 6 hours for *L. lactis*. We note the FBA models at the
inoculation predicted the highest growth rate for *L. lactis*, and the lowest for
*P. freudenreichii*, consistent with these experimental observations (Table 2).
In a second experiment, metabolite dosage was performed for *P. freudenre-
ichii* at the final time point, showing a large production of propionate (8g/L
of milk), a strong consumption of lactate (9g/L of milk) and a lower production of acetate and succinate (respectively 3,07 g/L and 0.371g/L of milk) (Supp. Table S2).

We used the experimental cultures of the strains to calibrate dynamic models of their respective metabolism. Experimentally, we observed that the metabolism of LAB is still active at the steady state. As lactic acid is produced from lactose, and lactose was not a limiting substrate, we therefore regulated the flux of lactic acid and acetate production by the flux of lactose consumption. In addition, following the approach of (Özcan et al., 2020), lactose consumption was further regulated by undissociated lactic acid (see Supplementary material §B.5). For *P. freudenreichii*, additional bounds on metabolite production were estimated on the acid dosage data described above in order to constraint production fluxes to physiological values (see Supplementary material §B.7). Furthermore, as carbon sources were not limited in the growth media, a logistic growth was assumed to model the plateau phase.

To calibrate the metabolic models, we estimated dynamic model parameters on these mono-culture data by solving the minimization problems (9)-(10). To avoid over-fitting, a reduced number of free parameters was kept for each model. Namely, only two parameters were fitted for *L. lactis*, tuning pH-dependent regulation, three parameters for *L. plantarum* involved in pH-dependent regulation and logistic growth, and one parameter for *P. freudenreichii* tuning a logistic regulation (see Material and Methods §4.3 for de-
tails).

After optimization, the dynamic FBA (dFBA) simulations fitted with the experimental growth data (Fig. 2, 3 and 4, (b)). The pH curves (Fig. 2, 3, (c)) are also correctly rendered, showing a good agreement of the pH model, together with *P. freudenreichii* metabolite production (Fig. 4, (c)). An overall coefficient of determination of 0.99 was recovered (Supplementary data, Fig. S1, left panel).

In order to validate that optimized FBA models correctly render known
metabolism, fluxes of pathways of interest described in Figure 1 are repre-
sented in Figure 2, 3 and 4 (a) at inoculation, when no metabolite is limiting.

The *L. lactis* model showed preferential fluxes of consumption of lactose *via* the tagatose pathway rather than Leloir pathway and production of L- and D-lactic acid *via* the glycolysis pathway, which explains the acidification of the milk and curd during cheese making (Fig. 2a). Production fluxes of diacetyl and butanediol occurred through an activation of acetolactate metabolism. An uptake flux of citrate by the *CITL* reaction encoded by the *citL* gene was observed.

The *L. plantarum* model presented a consumption flux of lactose and production fluxes of L- and D-lactic acid through the glycolysis pathway and acetic acid *via* the acetatekinase pathway, indicating that *L. plantarum* contributes to cheese acidification (Fig. 1b). The model shows different values of fluxes at different time points through the heterolactic (transketolase pathway) and homolactic (glycolysis pathway) metabolism. The activation
of the transketolase pathway in a dairy context was, to our knowledge, ob-
served for the first time. The consumption of citrate is strain dependant
for *L. plantarum* (Palles et al., 1998) and was observed in our case through
the *FEDCabc* transport reaction, leading citrate degradation to increase the
pyruvate pool.

The *P. freudenreichii* model showed a higher consumption flux of lactose
compared to D-lactic acid (Fig. 4a). TCA and Wood-Werkman cycles are
activated simultaneously (Deborde and Boyaval, 2000): the former allows
the regeneration of carbon sources necessary for the activation of the Wood-
Werkman cycle and thus, the production of propionate.
<figure><img src="image_3.png" /><figcaption>**Figure 2:** *L. lactis* GEM fitted on pure culture data. (a) FBA optimization with optimal parameters applied to the central carbon metabolism at the inoculation step. The color scale represents the reaction flux values predicted by the FBA and normalized by the highest flux value of the illustrated pathways. (b) Dynamics of *L. lactis* population in pure culture after parameter inference in the model (green line) and in the experimental data (orange lines, experiment replicates). (c) pH in pure-culture of *L. lactis* in the model (green line) and in the experimental data (orange lines, experimental replicates). *Com pound abbreviations:* ac: acetate; actn__R: acetoine; alac__S: (S)-2-Acetolactate; btd_RR: butanediol; cit: citrate; dgal6p: d-galactose-6-phosphate; dhap: Dihydroxyacetone phosphate; diact: diacetyl; f6p: fructose-6-phosphate; fdp: D-Fructose 1,6-bisphosphate; g1p: glucose-1-phosphate; g3p: Glyceraldehyde 3-phosphate; g6p: glucose-6-phosphate; gal: galactose; gal1p: galactose-1-phosphate; glc__D: glucose; lac__D, lac__L: lactate; lac6p: lactose-6-phosphate; lcts: lactose; oaa: Oxaloacetate; pyr: pyruvate; tag6p__D: D-tagatose-6-phosphate; tagdp__D: D-Tagatose 1,6-biphosphate. *Reaction abbreviations:* 6PGALSZ: 6-phospho-beta-galactosidase; ACLD: Acetolactate decarboxylase; ACTD2: Acetoin dehydrogenase; ACLDC: Acetolactate decarboxylase; ACLS: Acetolactate synthase; BTDD_RR: R R butanediol dehydrogenase; CITL: Citrate lyase; CITt4_1: Citrate transport via sodium symport; CS: Citrate synthase; D_LACt2: D lactate transport via proton symport; ENO: Enolase; FBA: Fructose-bisphosphate aldolase; GAL6PI: Galactose-6-phosphate isomerase; GALkr: Galactokinase; GAPD: Glyceraldehyde-3-phosphate dehydrogenase; HEX1: Hexokinase (D-glucose:ATP); L_LACt2r: L lactate reversible transport via proton symport; LACpts: Lactose transport via PEP:Pyr PTS; LACZ: B-galactosidase; LDH_D: D-lactate dehydrogenase; LDH_L: L-lactate dehydrogenase; PC: Pyruvate carboxylase; PFK: Phosphofructokinase; PFK_2: Phosphofructokinase; PGI: Glucose-6-phosphate isomerase; PGK: Phosphoglycerate kinase; PGM: Phosphoglycerate mutase; PYK: Pyruvate kinase; TGBPA: Tagatose-bisphosphate aldolase; UGLT: UDPglucose-hexose-1-phosphate uridylyltransferase</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3:** *L. plantarum* GEM fitted on pure culture data. (a) FBA optimization with optimal parameters applied to the central carbon metabolism at the inoculation step. The color scale represents the reaction flux values predicted by the FBA and normalized by the highest flux value of the illustrated pathways. (b) Dynamics of *L. plantarum* population in pure culture after parameter inference in the model (green line) and in the experimental data (orange lines, experiment replicates). (c) pH in pure-culture of *L. plantarum* in the model (green line) and in the experimental data (orange lines, experiment replicates). 6pgc: 6-Phospho-D-gluconate; 6pgl: 6-phospho-D-glucono-1,5-lactone; ac: acetate; acald: acetaldehyde; accoa: acetyl-coa; actn__R: acetoine; alac__S: (S)-2-Acetolactate; btd_RR: butanediol; cit: citrate; dhap: Dihydroxyacetone phosphate; f6p: fructose-6-phosphate; fdp: D-Fructose 1,6-bisphosphate; glp: glucose-1-phosphate; g3p: Glyceraldehyde 3-phosphate; g6p: glucose-6-phosphate; gal: galactose; gal1p: galactose-1-phosphate; glc__D: glucose; lac__D, lac__L: lactate; lcts: lactose; oaa: Oxaloacetate; pyr: pyruvate; ru5p__D: D-Ribulose 5-phosphate; tag6p__D: D-tagatose-6-phosphate; tagdp__D: D-Tagatose 1,6-biphosphate; udpg: UDPglucose; xu5p__D: D-Xylulose 5-phosphate. *Reaction abbreviations:* ACLD: Acetolactate decarboxylase; ACTD2: Acetoin dehydrogenase; ACLDC: Acetolactate decarboxylase; ACLS: Acetolactate synthase; BTDD_RR: R R butanediol dehydrogenase; D_LACt2: D lactate transport via proton symport; ENO: Enolase; FBA: Fructose-bisphosphate aldolase; FEDCabc: FED-Cabc; G6PDH2r: Glucose 6-phosphate dehydrogenase; GALkr: Galactokinase; GALUi: UTP-glucose-1-phosphate uridylyltransferase (irreversible); GAPD: Glyceraldehyde-3-phosphate dehydrogenase; GND: Phosphogluconate dehydrogenase; HEX1: Hexokinase (D-glucose:ATP); L_LACt2r: L lactate reversible transport via proton symport; LACZ: B-galactosidase; LCTSt3ipp: Lactose transport via proton aniport (periplasm); LDH__D: D-lactate dehydrogenase; LDH__L: L-lactate dehydrogenase; PC: Pyruvate carboxylase; PFL: Pyruvate formate lyase; PTA: Phosphotransacetylase; PGL: 6-phosphogluconolactonase; PFK: Phosphofructokinase; PGK: Phosphoglycerate kinase; PGI: Glucose-6-phosphate isomerase; PGM: Phosphoglycerate mutase; PGMT: Phosphoglucomutase; phosphoketolase: phosphoketolase; PYK: Pyruvate kinase; RPE: Ribulose 5-phosphate 3-epimerase; UGLT: UDPglucose-hexose-1-phosphate uridylyltransferase.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: *P. freudenreichii* GEM fitted on pure culture data. (a) FBA optimization with optimal parameters applied to the central carbon metabolism at the inoculation step. The color scale represents the reaction flux values predicted by the FBA and normalized by the highest flux value of the illustrated pathways. (b) Dynamics of *P. freudenreichii* population size in pure culture after parameter inference in the model (green line) and in the experimental data (orange lines, experiment replicates). (c) Lactate, acetate, propionate and succinate concentrations simulated (green points) and observed (orange points) at the final time point of a second experiment with different lactate inoculation (see Material and method) *Compound abbreviations*: ac: acetate; akg: 2-Oxoglutarate; fum: fumarate; glc\_D: glucose; lac\_D, lac\_L: lactate; lcts: lactose; mmcoa\_R: (R)-Methylmalonyl-CoA; mmcoa\_S: (S)-Methylmalonyl-CoA; oaa: Oxaloacetate; ppa: propionate.; ppap: Propanoyl phosphate; ppcoa: propionyl-coa; pyr: pyruvate; succ: succinate; succoa: succinyl-coa. *Reaction abbreviations*: 2131pyrpp: Methylmalonyl-CoA carboxyltransferase 5S subunit; CS: Citrate synthase; D\_LACt2: D lactate transport via proton symport; FRD2rpp: Fumarate reductase / succinate dehydrogenase (irreversible) (periplasmic, membrane potential dissipating); L\_LACt2r: L lactate reversible transport via proton symport; LACZ: B-galactosidase; LDH\_D: D-lactate dehydrogenase; LDH\_L: L-lactate dehydrogenase; MME: Methylmalonyl-CoA epimerase; MMM2: Methylmalonyl-CoA mutase; OOR3r: 2-oxoglutarate synthase (rev); PC: Pyruvate carboxylase; PPAKr: Propionate kinase; PPCSCT: Propanoyl-CoA: succinate CoA-transferase; PTA2: Phosphate acetyltransferase.</figcaption></figure>
## 2.3 Follow-up of metabolites production in co-culture experiments by the community-wide model.

Monitoring of growth and metabolic productions was performed during cheese production. *L. lactis*, *L. plantarum*, and *P. freudenreichii* reached respectively a culturability of 8.45 log₁₀ CFU g⁻¹, 8.47 log₁₀ CFU g⁻¹ and 8.59 log₁₀ CFU g⁻¹ in approximately 1200 hours. We built a community model by gathering the fitted individual models without further optimization, but the inclusion of the new plateau phase values in the logistic growth regulation in eq. 5, in order to reproduce *in silico* the metabolic behavior of the community (see eq. (3)-(4) and Material and methods §4.2). We note that, since no further fitting was performed, the co-culture experiments can be considered as a testing set involving unseen validation points.

The community model predicted the growth of each species quite accurately (Fig. 5(a)-(c), Supp. Table S4) despite a slight over-estimation of *P. freudenreichii* growth during the exponential phase. Metabolites dynamics were also correctly rendered. Lactose consumption is slower in the model compared to experiments, while citrate dynamics fits well with experimental data. Lactate production is slightly over-predicted during exponential phase, resulting in a slight over-acidification of the cheese during the molding and demolding steps. Acetate production is under-predicted during the early phases of cheese manufacture, but production rate during cheese maturation is correctly modeled. The main part of these shortcomings may be due to the
over-acidification of the medium during the first few hours after inoculation: the pH drops too sharply during the exponential phase of the LABs, which rapidly activates the pH regulation of their lactose uptake, inducing discrepancies in lactose curve until demoulding. However, during cheese maturation, these differences are quickly reduced and model outputs accurately predict the experimental data throughout ripening. Curves for propionate, and to a lesser extent succinate, are accurately predicted by the model, suggesting that production of organoleptic compounds by *P. freudenreichii* is correctly captured. Hence, the main characteristics shaping cheese flavor are correctly rendered by the model at the end of ripening (acidification, propionate, succinate, and to a less extent acetate). A global *R*² of 0.98 is obtained for this validation step (Supplementary data, Figure S1, right panel). In Figure S1 (right panel), three outliers can be observed for lactose and citrate depletion. They correspond to time points where the model predicts a depletion of the compounds while metabolite residuals can be measured in the experimental samples. They occur at C0w for citrate, indicating that the model predicts citrate depletion slightly too early, and at C4w and C7w for lactose, indicating little inaccuracies for small lactose amounts.

Furthermore, the metabolic pathways activated during mono-culture are also activated in co-culture (Supp. Fig. S3).
<figure><img src="image_6.png" /><figcaption>Figure 5: (a-c) Growth of respectively *L. lactis*, *P. freudenreichii* and *L. plantarum* species computed by the community model (green lines) or experimentally-observed (orange lines, experiment replicates). (d-j) Computational (green lines) and experimental (orange lines, experiment replicates) co-culture pH and metabolic profiles. Abbreviations: Linoc, Lactic acid bacteria inoculation (t = 0); Finoc, *P. freudenreichii* inoculation (t = 18 hours); Cm, molding (t = 19.5 hours); Cdm, demolding (t = 40 hours); C0w, start of ripening (t = 60 hours); C4w, Fourth week ripening (t = 732 hours); C7w, seventh week ripening (t = 1236 hours).</figcaption></figure>
## 2.4 Dynamical interactions within the community

Screening metabolite consumption and production during co-culture can reveal microbe-microbe interactions, such as cross-feeding or nutritional competition. We first surveyed the dynamics of metabolite exchange fluxes for each bacterium (Fig. 6, (a)) by computing, at each time step of the dynamics, the population-wide exchange flux $\mu_{i,j}b_i$ for metabolite j and bacterium i and the associated predicted metabolite concentrations. For citrate and lactose, *L. lactis* shows a strong uptake during its exponential phase, which rapidly vanishes after complete substrate depletion for citrate and activation of the pH-dependant regulation for lactose. Interestingly, *L. plantarum* fluxes are negligible for this substrate, suggesting that other metabolites support its growth, such as e.g. amino-acids. This may be due to the slower growth of *L. plantarum* during the first hours of cheese manufacture: citrate is depleted and acidification activates lactose consumption down-regulation when *L. plantarum* is still at low density. Lactose depletion is achieved by *P. freudenreichii* during the early phase of cheese ripening.

Lactate production reflects lactose consumption: lactate is strongly produced by *L. lactis* populations until $t = 250h$ (ripening), while only a slight production by *L. plantarum* can be observed. Lactate is however consumed by *P. freudenreichii* during ripening, which contributes to the pH increase of the system observed in Figure 5 d. Acetate is produced by the three bacteria: *L. lactis* shows a strong acetate production until $t = 250h$, while *L. plantarum*'s contribution increases slowly until a plateau phase and drops
at $t = 250h$. After reaching its plateau, *P. freudenreichii* keeps a constant rate of acetate production during cheese maturation. The same pattern is observed for butanediol production by the two LAB. Diacetyl is only produced by *L. lactis*, and its production stops after lactose depletion. Propionate and succinate are produced by *P. freudenreichii*, which keeps a constant release rate during ripening.

For each strain, we additionally computed the total net exchange fluxes of metabolites by integrating over time the FBA prediction $\int_0^T \mu_{i,j}(t)b_i(t)dt$. After normalization of this net exchange by the total exchange among the community (i.e. the sum of the individual net exchanges), we obtained the contribution of each strain to the metabolite fate (Fig. 6, (b)). These contributions confirm that succinate and propionate are solely produced by *P. freudenreichii*, while diacetyl (production) and citrate (consumption) are metabolized by *L. lactis*. *L. plantarum* contributes to butanediol production, although the main producer is predicted to be *L. lactis*, and to a smaller extent to acetate production, together with *L. lactis* and *P. freudenreichii*. Cross-feeding is observed for lactate between *L. lactis* (producer) and *P. freudenreichii* (consumer). Lactose consumption is shared between *L. lactis* and *P. freudenreichii*: interestingly, *P. freudenreichii* starts consuming lactose after *L. lactis* stops metabolizing this substrate, indicating a time segregation of lactose use, and no direct competition for this substrate.

The above model demonstrates that capturing the complexity of metabolic processes in the bacterial community during cheese making requires consid-
eration of the underlying dynamics occurring in the system. We wondered
whether applying *a priori* computational approaches relying on genome-
scale metabolic models could highlight additional putative interactions be-
tween the strains. To that end, we used two flux-based methods, namely
SMETANA (Zelezniak et al., 2015) and MICOM Diener et al. (2020a), and
one reasoning-based approach, Metage2Metabo (Belcour et al., 2020) in order
to suggest metabolic complementarity among the consortium.

MICOM and SMETANA highlighted respectively 14 and 25 putatively
exchanged metabolites (Supp. Table S6 and S7), while Metage2Metabo iden-
tified 11 metabolites that could not be produced by a species without inter-
actions among the community (Supp. Table S8). A first observation is that
a relatively high number of metabolites are common in the predictions of ex-
changeable compounds of SMETANA and MICOM: lactate, which was also
predicted in the dynamic model, phenylalanine, serine, malate, succinate,
xanthine, H2S, 2-ketoglutarate, glycollate and acetaldehyde. Other predicted
exchanges include mostly additional amino-acids (isoleucine, proline, glycine,
alanine). We used metatranscriptomic data in order to assess the validity of
the most plausible interactions according to the SMETANA score (see Meth-
ods): H2S (from L. lactis and L. plantarum to P. freudenreichii), ribose (from
L. lactis and plantarum to P. freudenreichii), glycerol (from L. lactis and
L. plantarum to P. freudenreichii) and phenylanaline (from L. plantarum to
P. freudenreichii). We verified the expression of genes associated to the pro-
duction of the metabolite in the donor species, and to the consumption of the
compound in the receiver species (see Methods and Supp. Fig. S2). Results
suggest that interactions for H₂S, ribose and glycerol are plausible at several
steps of cheese making. Conversely, whereas gene expression data shows that
phenylalanine consumption pathways are highly expressed in *P. freudenre-
ichii*, their production by *L. plantarum* is not highly expressed suggesting
that this interaction is less effective. Finally, the set of metabolites whose
production is predicted by Metage2Metabo to require interactions among
the community include various fatty-acids, as well as galactose-1-phosphate,
benzyl-CoA, glyceraldehyde and xantosine, indicative of metabolic comple-
mentarity between the strains for the related metabolic pathways.
<figure><img src="image_7.png" /><figcaption>Figure 6: (a) **Consumption and production flux dynamics.** We represent for each bacterium the dynamics of the consumption and production fluxes of the followed-up metabolites (colored plain lines). We also represent the dynamics of the concentration of the corresponding metabolite in the cheese (dashed black line, right axis for scale). (b) **Bacterial relative contribution to the metabolite fate.** We represent the relative contribution of each bacteria to the final metabolite concentration. Namely, we computed each microbial overall contribution by integrating in time its consumption or production flux ($\int_0^t \mu_i j b_i$ for metabolite j and bacteria i, see eq.4), and normalized the result by the sum over the three bacteria. The value 1 (or -1) then represents the net production (or consumption), i.e. the difference between the final and initial concentration of the metabolite (horizontal dashed lines). For each metabolite, the bacterial contribution is represented by a colored bar, with positive value for production, and negative value for consumption. Abbreviations: Linoc, Lactic acid bacteria inoculation (t = 0); Finoc, *P. freudenreichii* inoculation (t = 18 hours); Cm, molding (t = 19.5 hours); Cdm, demolding (t = 40 hours); C0w, start of ripening (t = 60 hours); C4w, Fourth week ripening (t = 732 hours); C7w, seventh week ripening (t = 1236 hours).</figcaption></figure>
# 3 Discussion

In this work we provide an accurate model of the bacterial metabolism occurring in milk during cheese making. Using metabolic models of the three inoculated bacterial strains and multi-omics data, we were able to accurately reproduce metabolite production patterns and the dynamics of the bacterial populations. We relied on pure culture experiments to calibrate the individual metabolic models with few parameters for each, thereby limiting overfitting during community simulations. An originality of the model is its ability to predict the dynamics over the entire cheese making process (seven weeks). Our work supplies hypotheses on the underlying metabolic pathways at stake in the milk environment, as well as the contribution of each strain in the consumption of nutrients and the producibility of organoleptic compounds.

Knowledge from microbiology experts, scientific literature, and multi-omics data were necessary to achieve our high-quality genome-scale metabolic models and subsequent dynamic simulations. The central carbon metabolism of LAB mainly produces lactic acid via glycolysis, tagatose and Leloir pathways from lactose degradation (Widyastuti et al., 2014; Van Rooijen et al., 1991; Kleerebezem et al., 2003); for the first time, citrate utilization and heterolactic fermentation were observed in the curated model of *L. plantarum* in milk. Ensuring the individual metabolic models can activate these pathways was therefore an important step of their validation, which enabled
reproducing the enzymatic activities observed in (Quatravaux et al., 2006; Carroll et al., 1999). The propionic bacterium *P. freudenreichii* converts lactose, and preferentially lactate, into propionate according to (Loux et al., 2015; Thierry et al., 2011). Following (Borghei et al., 2021), we added the propionyl-CoA:succinate CoA transferase (2.8.3.-) enzyme for completing the Wood-Werkman cycle and thus enabling propionate production.

In (Özcan et al., 2020), individual tuning parameters of all media com-
pounds are imposed in order to explain metabolomics at the species level, and
to predict it at the community level. Compared to this work, we implemented
a dynamic FBA based on (Mahadevan et al., 2002), that used strain specific
parameters for predicting growth, pH and metabolic concentrations. Namely,
we supplemented the classical dFBA framework with dynamic bounds on
the uptake and release of key compounds and with a population-sensitive
regulation of the growth rate. These improvements aimed to capture reg-
ulations that are not directly related to the availability of metabolites in
the growth media, and that would be missed with a classical dFBA model.
Rather, these regulations, that depend on the physico-chemical state of the
environment *i.e.* pH and bacterial density, can be modelled by tuning state-
dependent dynamic bounds for exchanged metabolites and growth rate. To
avoid over-fitting, we narrowed down the number of inferred parameters for
each bacterium, keeping only two parameters tuning pH regulation for the
LAB, and one parameter driving the growth for *L. plantarum* and *P. freuden-
reichii*. Additional calibration was conducted for *P. freudenreichii*, without
any additional mathematical optimization, by computing upper bounds for
metabolites from metabolomic data obtained in mono-culture.

The co-culture model gave insights into the community behaviour dur-
ing cheese manufacture, and suggested temporality and contribution of each
species to aroma coumpounds production. It showed that propionate and
succinate production could only be attributed to *P. freudenreichii* in line
with (Cao et al., 2021), and that diacetyl seems to be only produced by
*L. lactis*. According to the community model, butanediol was produced dur-
ing the moulding and at the beginning of the ripening by both *L. lactis* and
*L. plantarum*, as observed in the respective mono-culture models. Regard-
ing acetate production, the model infers an early production from moulding
to ripening steps by *L. lactis* and *L. plantarum* followed by the major in-
volvement of *P. freudenreichii* during ripening. Concerning carbon source
utilization, while *L. plantarum* is completely equipped for lactose and cit-
rate metabolisms, which are fully activated in mono-culture, it appears that
these metabolic pathways are strongly down-regulated in co-culture. When
cultivated with *L. lactis*, as *L. plantarum* growth is the slowest of both, it
achieves its maximal metabolic efficiency after *L. lactis*'s exponential phase.
As *L. lactis* mainly grows on lactose producing lactate, its exponential phase
is associated with a strong acidification of the environment, which activates
the pH-dependant repression of lactose metabolism in *L. plantarum*. *L. plan-
tarum* may use amino-acids or other metabolites not followed-up by the
dynamical model to support its growth. To further explore this hypothe-
sis, we used community-scale metabolic model exploration tools to identify
metabolic interactions outside the scope of the compounds we screened in
the dynamical co-culture model, highlighting other putative metabolic inter-
actions, of which a high proportion related to amino-acids.

In the co-culture model, *P. freudenreichii* keeps an active metabolism in
the later phase of cheese manufacturing, during ripening. This metabolism is
reflected by a lactate consumption all along the ripening phase, and a release
of organoleptic compounds at a constant rate (Fig. 6). This behaviour is
coherent with known metabolic capabilities of *P. freudenreichii* and with
the metabolomics data obtained during cheese processing. In the model,
*P. freudenreichii* also consumes the lactose remaining after *L. lactis* growth.
Indeed, during its exponential phase, *L. lactis* mainly metabolizes lactose,
strongly producing lactate and making the pH drop down. The pH reduction
inhibits the pH-regulated lactose uptake in the two LABs, hence stopping
lactose consumption by the *Lactobacilli*. Since lactose remains in the media
and *P. freudenreichii* is the unique bacteria still capable to use it, it completes
lactose depletion. In metatranscriptomics data, the LACZ gene is highly
activated in *P. freudenreichii* during ripening, supporting lactose degradation
by this strain after acidification (see Fig. S2).

This work highlights how the combination of expert knowledge, omic data
and metabolic modeling can provide precise and accurate predictions on the
mechanisms responsible for the dynamic behaviour of a microbial community.
On the other hand, it also points out that the amount of data and the efforts
necessary to create such high-quality models still remain an expensive price
to pay despite the improvement of simulation approaches. Methodological
development are yet to be proposed in order to automatize the calibration of
models with data, and ensure both the correctness of the predictions and scal-
ability to larger communities or communities of empirical composition. This
is a necessary step towards effective digital twins of bacterial metabolism for
cheese production, i.e. digital models accurate enough to produce actionable
knowledge. In the context of microbial fermentation during food processing,
digital twins could be used for instance to define environmental controls on
the microbial communities, such as maturation, temperature or salt. Those
could be screened and optimized with the model to achieve a specific target,
such as an organoleptic profile of the cheese.

# 4 Material and methods

## 4.1 Biological data

**Community composition.** Strains used in this project are described in
(Cao et al., 2021). Briefly, the controlled bacterial community is composed
of two lactic acid bacteria (LAB), *L. lactis* subsp. *lactis* biovar *diacetylac-
tis* CIRM-BIA1206 (termed *L. lactis*), *L. plantarum* CIRM-BIA465 (termed
*L. plantarum*), and one propionic bacterium, *P. freudenreichii* CIRM-BIA122
(termed *P. freudenreichii*) provided by the International Centre for Microbial
Resources-Food Associated Bacteria (CIRM-BIA<fnref n="1" />).

**Cheese manufacturing and sample collection.** The standard cheese manufacture protocol described in (Cao et al., 2021) is used. Briefly, after collection of fresh milk, pasteurization (76°C, 20 seconds), skimming, standardization (until 30 grams of fat and 36 mg of calcium per kg of milk, *L. lactis* and *L. plantarum* were inoculated at 5.7 and 5.2 log₁₀ CFU g⁻<fnref n="1" /> (Colony Forming Unit). After prematuration of the inoculated milk (at 14 °C for 18 h), *P. freudenreichii* was inoculated at 6.1 log₁₀ CFU g⁻<fnref n="1" /> and stirred for 20 min. The milk was then poured in vats and warmed at 33°C for 30 min. Commercial rennet was then added and coagulum was cut, stirred, washed, drained, prepressed and transferred into moulds for pressing. Demolding occurred on the third day of cheese making, and the cheeses were salted (immersion at 12°C for 10h in saturated brine), dried overnight and vacuum-packed in plastic bags on the fourth day of cheese making for ripening during 7 weeks at 13 °C.

After LAB bacteria inoculation (Linoc, t=0h) and *P. freudenreichii* inoculation (Finoc, t=18h), samples were collected at 5 different production stages for further analysis as described in (Cao et al., 2021): molding stage (Cm, t=19.5h), demolding stage (Cdm, t=40h), start of ripening (C0w, t=60h) and after 4 and 7 weeks of ripening (C4w, t=732h and C7w, t=1236h). Cheese processing was replicated four times with independent inocula and

<footnote n="1">¹https://collection-cirmbia.fr/</footnote>
milks at different days. Samples were collected at each production step for
numeration (4 replicates), biochemistry (2 replicates) and RNA sequencing
(2 replicates). All the data but RNA-seq and single strain experiments were
retrieved from Cao et al. (2021). We recall below the different data collection
protocols to ease the reading.

**Single strain culture experiments.** In order to calibrate metabolic models, a training dataset is built using single strain culture experiments. Single strain growth experiments of the three bacteria were performed in UHT full fat milk (Delisse). Dairy cultures were performed in duplicate in 15 ml Falcon tubes under microaerobic conditions, without shaking and using one tube per sampling time to avoid oxydative stresses, at 30°C. For *P. freudenreichii*, the milk was supplemented with lactate (20 g/L) and with casein peptone (10 g/L). A second growth experiment was performed for *P. freudenreichii* in order to measure metabolite levels (acetate, lactate, propionate and succinate) at the final time. *L. lactis* was inoculated at 1 × 10⁶ and 4 × 10⁶ CFU g⁻¹, *L. plantarum* at 3 × 10⁶ and 6 × 10⁶ CFU g⁻¹ and *P. freudenreichii* at 1.4 × 10⁷ and 9 × 10⁶ CFU g⁻¹. Pure culture growth data are available in Table S1.

**Plate-counting** Following (Cao et al., 2021), the growth kinetics were followed and expressed as Colony Forming Unit per milliliter (CFU mL⁻¹) using the previously described micromethod (Baron et al., 2006) in biological replicates and technical duplicate. All strains were grown at 30°C under mi
croaerophilic conditions in standard broths: the M17 broth containing 0.5% (w/v) lactose (Difco, Beckton Dickinson and Company, Sparks, USA) for L. lactis, Man, Rogosa, and Sharpe broth (MRS, pH 5.4, Biokar Diagnostics, Beauvais, France) for L. plantarum, and YEL broth containing 1% lactate (Malik et al., 1968) for P. freudenreichii until the visualization of colonies from one to six days. Results were expressed as mean of technical duplicates.

**Biochemistry** Milk samples were directly diluted 40-fold in the $H_2SO_4$ solution and (0.05M) and centrifuged in Vivaspin-2 (Sartorius) and frozen at -20°C. Cheese samples were first blended in deionized water at 40°C (1:5 w/w) in a filter bag and incubated at 40°C for 1 h. The suspensions were filtered on Whatmann 40 paper. The filtrates were diluted six-fold for Cm, four-fold for Cdm curds and C0w cheeses, and two-fold for ripened cheeses (C4w and C7w) in the $H_2SO_4$ solution (0.05 M) and frozen at -20°C. All samples were defrosted at 4°C and then centrifuged (8500 g, 20 min, 4°C). Supernatants were filtered on Chromafil Xtra PVDF-45/13, 0.45 µm pore size (Macherey-Nagel GmbH & Co. KG, Germany). The compounds were separated using a Thermo scientific Ultimate 300 (Thermo scientific, Les Ulis France) equipped with a Rezek ROA-Organic Acid H+ (300*7.8mm) (Phenomenex, California), using 5 mM $H_2SO_4$ as the mobile phase, and detected using both UV at 210 nm (UVD 170UV Dionex) and/or refractometry (RI2031 plus, Jasco). The compounds were identified and quantified using
standard solutions of lactate, pyruvate, acetate, citrate and propionate, lactose, galactose (concentrations ranging from 0.06 to 1 g L⁻¹). Data were acquired and processed by using the Chromeleon software.

Volatile compounds were analyzed using headspace (HS) trap extraction, gas chromatography and mass spectrometry (GC-MS) according to (Pogačić et al., 2015). The samples were injected in a random order, with standards (mixture of eleven volatiles: four esters (ethyl acetate, ethyl propanoate, ethyl butanoate, and ethyl hexanoate), two aldehydes (3-methylbutanal and benzaldehyde), two ketones (2-heptanone and 2-nonanone), 2,3-butanedione, dimethyl disulfide, and 3-methylbutanol) and blank samples (boiled deionized water) to monitor possible MS drift and carryover. Compounds were eluted on an Elite-WAX ETR (30m x 0.25mm ID x 0.25µm, PerkinElmer USA) column. They were identified by comparing their mass spectral data and linear retention indexes (LRI) calculated on a polar column with that of reference standard compounds and with data from Library NIST 2008 (Scientific Instrument Services, Ringoes, NJ, USA) and PubChem. Raw data files were converted to time- and mass-aligned chromatographic peaks areas, using the XCMS open-source package implemented with the R statistical software. The volatile compounds were semi quantified using the abundance of one selected mass fragment (m/z), in arbitrary units.

**Metatranscriptomic data.** Ten grams of cheese were dispersed in 90g of a 2% sodium citrate solution during 1 min at max speed followed by 1 min
at low speed in stomacher (Humeau, Treillière, France) at room tempera-
ture. Ten milliliters of the dispersion were centrifuged at 10 000 g at 4°C
during 5 min. Cell pellets were stored at –80°C for RNA extraction. After
defrost, the samples were centrifuged again for 5 min at 4°C, 10 000 g to dis-
card the supernatant. Cell pellets were suspended in 200 µL lysis buffer (50
mM Tris–HCl, 1 mM EDTA; pH 8.0) 20 mg/g lysozyme (MP Biomedicals,
Illkirch, France) and 50 U/g mutanolysin (Sigma, Saint Quentin Fallavier,
France) and incubated for 30 min at 24°C. Suspensions were then transferred
to 2 mL tubes containing 50 mg zirconium beads (diameter: 0.1 mm; BioSpec
Products, Bartlesville) and 50 µL SDS (20 %). Samples were then shaken
in a Precellys Evolution (Bertin, Montigny-le-Bretonneux, France) for two
cycles of 40 s at 6500 rpm. RNA was then extracted from the cell lysates us-
ing Qiazol and RNeasy Mini kit (Qiagen) according to (Falentin et al., 2010)
with two or three steps of subsequent DNase treatment (Dnase Rnase free,
Ambion) according to the manufacturer's instructions. RNA integrity (RIN)
was evaluated using an Agilent 2100 bioanalyzer (Agilent Technologies, Santa
Clara, CA). RNA concentrations were quantified using Qubit Fluorometric
quantitation (Thermofisher, France). Lack of DNA contamination was also
checked using Qbit and DNA quantitation protocol according to manufac-
turer instructions. RNA samples with a RIN value greater than 8, indicating
a good RNA integrity, were kept for further analysis.

Metatranscriptomic data were obtained at five time steps of cheese man-
ufacturing. Illumina libraries and sequencing were performed by GeneWiz
(Leipzig, Germany). Reads are available at ENA (Cambridge) under accession number ENA: PRJEB42478. Sequenced reads were trimmed using Trim Galore and mapped against reference genomes with Bowtie 2 (Langmead et al., 2009). No mismatch was allowed and only reads that mapped to unique sequences were further analysed. Reads mapping to each coding sequence were counted using HTSeq-count (`-stranded=reverse`, `-a 0`, `-t CDS`, `-i db_xref`, `-m union`), (Putri et al., 2022).

Statistical analysis were conducted using R statistical software version 4.0.3 (R Core Team, 2021). RNA-Seq raw count data were normalized in two steps. First, a species-specific scaling factor was applied to eliminate composition biases between libraries. For each of the three species, they were calculated using all the genes of the considered species and two others rows for the other species calculated as the sum of the counts attributed to each of these two other species, using the method TMM (Trimmed Mean of M-Values) as implemented in the package edgeR version 3.32.1 (Robinson et al., 2009). A supplementary within-sample normalization step was performed to correct for gene length and enable comparison. Finally replicates were averaged.

**Targeted metabolomics.** Methods are available in (Cao et al., 2021). Briefly, sugars and organic acids in the samples were quantified using high performance liquid chromatography (HPLC). Volatile compounds were analysed using headspace (HS) trap extraction coupled to gas chromatography-
mass spectrometry (GC-MS).

## 4.2 Metabolic modeling

**Genome sequencing** Cell pellets (equivalent to $1e^{10}$ CFU) were obtained by centrifugation for 10 min at 5000 g from strain cultures in broth medium. DNA was extracted using DNeasy midi kits (Qiagen) according to the protocol described in (Falentin et al., 2010). Genomes were sequenced by Illumina NovaSeq 2 x 150 pb and PacBio Sequel (Eurofins Genomics, Constance, Germany). Genomes were assembled using Unicycler (Wick et al., 2017a,b) on the Galaxy Genotoul bioinformatics platform (Toulouse, France). Genome sequences were integrated in the MicroScope platform hosted at Genoscope (CEA, Evry, France) for automatic annotation according to (Vallenet et al., 2019). Annotated genomes are available at ENA (Cambridge) under accession number ENA: PRJEB54980.

**Genome-scale metabolic model reconstruction.** Individual bacterial genomes were functionally annotated using the MicroScope platform (Vallenet et al., 2019). Genome-scale metabolic network were then reconstructed from these annotated genomes with CarveMe (Machado et al., 2018) version 1.4.1. Each genome scale metabolic model was polished using ModellPolisher (Römer et al., 2016) which improved annotations by adding information related to BIGG identifiers from BiGG Models knowledgebase (King et al., 2015). The consistency of each model was checked with MEM-
OTE (Lieven et al., 2020) https://memote.io at the end of the polishing test. The medium composition included the minimal media found in CarveMe and compounds in milk https://ciqual.anses.fr/#/aliments/19024/lait-entier-pasteurise, with amino-acids were considered in excess (Supplementary File 1).

**Flux balance analysis.** A genome-scale metabolic model (GEM) was constructed for each strain using the metabolic network and the medium composition, and analysed using Flux Balance Analysis (FBA) (Orth et al., 2010) FBA computes an optimal flux distribution satisfying stoichiometric and biophysical feasibility providing maximal biomass, by solving the linear programming problem

$$ \begin{array}{ll} \text{find } v^* \in \mathbb{R}^{N_r} \text{ such that } v^* = & \underset{S.v=0}{\arg\max} \quad v_{growth} \quad (1) \\ & v_{min}^{in} \le v \le v_{max}^{in} \\ & v_{min}^{ex} \le v \le v_{max}^{ex} \end{array} $$

where $v_{growth}$ is the flux of the model biomass reaction, $S$ is the stoichiometry matrix derived from the metabolic model, $v_{min}^{in}$ and $v_{max}^{in}$ (resp. $v_{min}^{ex}$ and $v_{max}^{ex}$) are minimal and maximal values for the intracellular (resp. exchange) fluxes defining flux bounds for the $\mathbb{R}^{N_r}$-dimensional flux vector $v$. The value of $(v_{min}^{in}, v_{max}^{in})$ and $(v_{min}^{ex}, v_{max}^{ex})$ set in the GEMs are available in Supplementary Material.
Equation (1) then defines, for each bacterium *i*, the mapping $\mu_i$

$$ \mu_i(v_{min}^{ex}, v_{max}^{ex}) = v^* \qquad (2) $$

between the constraint vectors on exchange reactions defining the nutritional environment and the optimal flux $v^*$ obtained in this nutritional context.

The FBA optimization problem (1) was solved using CobraPy (Ebrahim et al., 2013) version 0.17.1 using a defined lexicographic order (Gomez et al., 2014) to guarantee a unique solution.

**GEM manual refinement** A manual refinement step was performed for all models to ensure organoleptic compound production and the activity of metabolic functions of interest with respect to literature knowledge. The resulting modifications can be found in Supplementary Material Table S3. As the three species are micro-aerobic, an uptake bound of oxygen is defined as the minimal value which permitted the production of aroma compounds.

**Dynamic modeling** The dynamic behavior of each culture system was computed using dynamic flux balance analysis (dFBA) (Mahadevan et al., 2002) supplemented by additional population dynamics mechanisms.

We note $b = (b_i)$ for $i \in \mathcal{B} = \{L. lactis, L. plantarum, P. freudenreichii\}$ the bacterial population densities, and $m = (m_j)_{1 \le j \le N_m}$ the concentration of the $N_m$ metabolic compounds that are screened by the model. Note that these metabolites are taken from a short list of compounds that are uptaken
or released by at least one of the three micro-organisms and that are particularly relevant from a biological point of view (namely citrate, acetate, diacetyle, succinate, lactose, and lactate). These metabolites hence couple the metabolic models of the three bacteria.. We set the system, for $i \in \mathcal{B}$ and $1 \le j \le N_m$

$$ \partial_t b_i = \mathcal{R}_i(b_i)\mu_{i,i}((c_{min,i}^{ex}, c_{max,i}^{ex})(b,m))b_i \quad (3) $$

$$ \partial_t m_j = \sum_{i \in \mathcal{B}} \mu_{i,j}((c_{min,i}^{ex}, c_{max,i}^{ex})(b,m))b_i \quad (4) $$

where the term $\mu_{i,j}$ is the component corresponding to metabolite $j$ (or biomass $i$) of $\mu_i$, computed from mapping 2 given the set of constraints $(c_{min,i}^{ex}, c_{max,i}^{ex})(b,m)$. These constraints depend dynamically on the state variables $b$ and $m$.

The vector $\mathcal{R}_i$ models a population-sensitive regulation process on the population growth in a phenomenological way with the term

$$ \mathcal{R}_i = \lambda_i \left(1 - \frac{b_i}{\beta_i}\right), \quad \text{for } i \in \mathcal{B} \qquad (5) $$

where $\lambda_i$ is a weight and $\beta_i$ is a population carrying capacity which is the plateau phase value of the corresponding bacteria in the experiments.

For a given substrate $j$ and the metabolic model $i$, the usual bounds are

$$ c_{min,i,j}^{ex} = \max\left(-\frac{m_j}{\Delta t \sum_{i \in \mathcal{M}_j} b_i}, v_{i,j}^{\text{int}}\right), \quad c_{max,i,j}^{ex} = 0, \quad (6) $$
where $M_j$ is the subset of bacteria metabolizing $j$ and $\Delta t$ is an import characteristic time. These equations reflect uniform sharing of resources among bacteria during the time window $\Delta t$. When the substrate is in excess, the flux is limited by the intrinsic import limit $v_{i,j}^{int}$.

Classical dFBA did not allow by itself to capture known microbial physiology features during growth in milk such as pH regulation on lactose import. We then added additional regulations on some substrates and products. Following (Özcan et al., 2020), as the LAB bacteria metabolism of lactose is down-regulated by undissociated lactic acid with exponential decrease, for $j = lcts_e$ and $i \in \{L. lactis, plantarum\}$ we have

$$c_{min,i,j}^{ex} = \max\left(-\frac{m_{lcts_e}}{\Delta t * \sum_{l \in M(lcts_e)} b_l}, -\mu_{max,lcts}^i * 10^{(-k_{lac} * \phi_{undiss})} - \mu_{min,lcts}^i\right) \quad (7)$$

where $k_{lac}$ is an exponential decay, $\mu_{max,lcts}$ and $\mu_{min,lcts}$ are maximal and minimal values and $\phi_{undiss}$ is a function computing the concentration of undissociated lactic acid from lactate. The function $\phi_{undiss}$ is derived from Henderson-Hasselbalch equation as in (Özcan et al., 2020) and reads

$$\phi_{undiss}(m_{lac\_L_e}, m_{lac\_D_e}) = \frac{m_{lac\_L_e} + m_{lac\_D_e}}{1 + 10^{c_1*(m_{lac\_L_e} + m_{lac\_D_e}) + c_2}} \quad (8)$$

In this equation, the term $c_1 * (m_{lac\_L_e} + m_{lac\_D_e}) + c_2$ is a linear approximation of $(pH - pKa)$, so that parameters $c_1$ and $c_2$ are inferred by linear regression from data of the co-culture experiment. Lactate production
is regulated in a similar way, while acetate production is regulated according
to lactose availability. See Supplementary material sec. B.5 for a detailed
description and justifications of the dynamical bounds set on the exchange
reactions.

For *P. freudenreichii*, production upper bounds were evaluated by computing per capita production flux from the metabolite dosages and the growth data in the pure culture experiments, assuming a constant flux (see Supplementary material B.7).

To sum up, the present model differs from usual dFBA framework by
three additional features: 1) the population-sensitive regulation in eq. (5),
2) the pH regulation in eq. (7), and 3) the regulation of certain metabolite
exports. These developments were necessary to take into account physiolog-
ical characteristics that are not directly related to metabolism but rather to
how bacteria sense their environment.

The initial condition, the inoculum and nutritional environment com-
posed of milk compounds, amino acids, and co-factors, were defined as spec-
ified in (Cao et al., 2021), after unit conversion (See supplementary material
for a list of metabolites). A cell mass of 0.33 × 10⁻¹² g CFU⁻¹ was applied for
biomass unit conversion. To derive this number, we assumed that a microbial
cell has a volume of 1 µm³ and a density of 0.33g cm⁻³. *P. freudenreichii*
growth was permitted in the community model simulation after 18 hours in
agreement with the inoculation time of the species (*Finoc* time point).
**pH simulation.** Following (Özcan et al., 2020), we obtain pH value from lactate concentration using the relation

$$pH = pKa_{lac} + c_1 * (m_{lac\_L_e} + m_{lac\_D_e}) + c_2$$

where the term $c_1$ and $c_2$ are common to every species and are obtained by linear regression of $pH - pKa_{lac}$ and lactate concentrations in the co-culture experiments, where $pKa_{lac}$ is the acid dissociation constant of lactate.

We note that

$$pH - pKa_{lac} = c_1 * (m_{lac\_L_e} + m_{lac\_D_e}) + c_2$$

so that eq.(8) can be rewritten with

$$\phi_{undiss}(m_{lac\_L_e}, m_{lac\_D_e}) = \frac{m_{lac\_L_e} + m_{lac\_D_e}}{1 + 10^{pH - pKa_{lac}}}.$$

Hence, a direct link can be done between pH, lactate concentrations $m_{lac\_L_e}$ and $m_{lac\_D_e}$, and the quantity $\phi_{undiss}$ that is used to define the pH-related regulation of lactose uptake.

**Numerical implementation** The dynamical system is solved by a semi-implicit Euler scheme. At time step $n$, after computation of $F_j$ the total flux
of compound j as defined by the set of GEMs,

$$F_j = \sum_{i \in \mathcal{B}} \mu_{i,j} \left( (c_{min,i}^{ex}, c_{max,i}^{ex})(b^n, m^n) \right) b_i,$$

we computed

$$b_i^{n+1} = b_i^n + \Delta t * F_{bi}$$
$$m_j^{n+1} = \begin{cases} m_i^n + \Delta t * F_j & \text{if } F_j > 0 \text{(explicit case)} \\ m_j^n / (1 - \Delta t * F_j / m_j^n) & \text{else (implicit case)} \end{cases}$$

This scheme guarantees solution positivity at all times.

The dFBA problems (3)–(4) and parameter inference are solved using custom Python scripts and the scientific libraries ` numpy ` v.1.19.5(Harris et al., 2020), ` pandas ` v.0.25.3 (McKinney et al., 2010), and ` scipy ` v.1.5.3 (Wes McKinney, 2010). During dFBA, the FBA is solved using a defined lexicographic order (Gomez et al., 2014).

## 4.3 Model fitting

To fit the model parameters, we adapted equations (3)–(4) to mono-culture experiments. We then inferred the model parameter vector $\theta$ of the resulting dFBA model by maximising the log-likelihood of observed data in the mono-culture. Namely, for $i \in \{L. lactis, L. plantarum\}$, we minimised the cost
function

$$J(b_i, pH|\theta_i, b_{i,exp}, pH_{exp}) = \left\| \frac{\log_{10}(b_i) - \log_{10}(b_{i,exp})}{\sigma_{\log,i,exp}} \right\|^2 + \alpha \left\| \frac{pH - pH_{exp}}{\sigma_{pH,exp}} \right\|^2 \quad (9)$$

where $\sigma_i$ (resp. $\sigma_{\log,i}$) is the standard deviation (resp. the log-transformed standard deviation) of the corresponding data, and $\theta_i = (k_{lac}, \mu_{max,lcts})$ if $i = L. lactis$ and $\theta_i = (k_{lac}, \mu_{max,lcts}, \lambda_i)$ if $i = L. plantarum$. For $P. freudenreichii$, the cost function is

$$J(b, m|\theta_i, b_{exp}, m_{exp}) = \left\| \frac{\log_{10}(b) - \log_{10}(b_{exp})}{\sigma_{\log,b,exp}} \right\|^2 + \alpha \left\| \frac{m - m_{exp}}{\sigma_{m,exp}} \right\|^2 \quad (10)$$

where $m_{exp}$ is the final dosage of acetate, lactate, propionate and succinate, and $\theta_i = (\lambda_i)$. We kept the number of fitted parameters small (2 parameters for $L. lactis$, 3 parameters for $L. plantarum$ and 1 parameter for $P. freudenreichii$), to avoid over-fitting. Optimization results are presented in Table S9.

## 4.4 Inference of metabolic interactions

In order to predict metabolic exchanges likely to occur in the microbial consortium, MICOM v0.32.0 (Diener et al., 2020b) and SMETANA v1.2.0 (Zelezniak et al., 2015) were performed on the milk medium. Additionally, the prediction of metabolic complementarity between the three metabolic models grown in milk was performed with Metage2Metabo v.1.5.0 (Belcour
et al., 2020). All metabolites predicted by the tools were classified into higher
level categories: compounds identifiers were mapped to the Metacyc v.25.5
database (Caspi et al., 2018), and the ontology of the compounds was used.
The list of predicted cross-feeding interactions is available in Supplementary
Table S8.

A subset of putative bacterial interactions were selected among the ones
predicted by SMETANA, according to their score (SMETANA score), for
which a threshold was set at 0.5. The relevance of selected exchanged
metabolites was assessed using metatranscriptomic data (RPKM and repli-
cated means). For each metabolite, reactions producing (in the producer
species) and consuming (in the consumer species) the metabolite were ex-
tracted from gene-protein-reactions associations. In order to evaluate the
putative expression of these genes, the expression values were compared to
the overall expression of genes in each species, separated into quartiles.

4.5 Visualisation and statistical analyses

Plots were generated using `matplotlib` version 3.3.4 (Hunter, 2007). Metabolic
maps for pathways of interest were generated with Escher (King et al., 2015)
v.1.7.3 on Python v.3.6.9, after running an FBA on each model. Fluxes of
reactions of interest for each species were normalized by the maximal flux
value of the targeted pathways in the corresponding species.
## 4.6 Code and data availability

Code for metabolic simulation and metabolic models is available in https://forgemia.inra.fr/tango/tango_models. Genomes of the three bacterial strains and their annotations are available in ENA under the accession number PRJEB54980. Metatranscriptomic data is available in ENA under the accession number PRJEB42478. Annotated genomes, RNA reads and RPKM are deposited in https://entrepot.recherche.data.gouv.fr/privateurl.xhtml?token=4bf01466-29dd-4d71-b25e-b775d9cc39dc. Biochemical dosage data and metabolomic data are available in (Cao et al., 2021).

## 5 Acknowledgements

The authors acknowledge the Galaxy Genotoul bioinformatics platform (Toulouse, France), the MicroScope platform hosted at Genoscope (CEA, Evry, France), and the GenOuest bioinformatics core facility (https://www.genouest.org) for providing the computing infrastructure. CF, DJS and SL were supported by the French National Research Agency (ANR) France 2030 PEPR Agroécologie et Numérique MISTIC ANR-22-PEAE-0011. CF and SL were supported by the Inria Exploratory Action SLIMMEST. We are grateful to the INRAE STLO Dairy Platform (https://www6.rennes.inrae.fr/plateforme_lait_eng/, doi: 10.15454/1.557240525786479E12) for providing help and support, to the CNIEL for partial funding and fruitful discussion.
sions on the interpretation of the results (TANGO Project).

# 6 Conflicts of interest

The authors declare no conflict of interest.

# 7 Authors contributions

Conceptualization: CF, HF, ML, SL; Data curation: JA, HF, ML; Formal analysis: CF, DJS, HF, ML, JA, WC, SL; Funding acquisition: HF; Investigation: CF, HF, ML, SL; Methodology: CF, HF, JA, ML, SL; Project administration: CF, DJS, HF, SL; Resources: HF; Software: CF, ML, SL; Supervision: CF, DJS, HF, SL; Validation: HF, CF, ML, SL; Visualization: CF, ML, SL; Roles/Writing - original draft: CF, HF, ML, SL; Writing - review & editing: all.

# 8 Supplementary information

The manuscript includes Supplementary Figures, Supplementary Material and a Supplementary File entitled *supplementary_file1.xlsx*.
# 9 Copyright

A CC-BY public copyright license has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission, in accordance with the grant's open access conditions.

# 10 Bibliography

## References

Abedi, E. and Hashemi, S. M. B. (2020). Lactic acid production – producing microorganisms and substrates sources-state of art. *Heliyon*, 6(10):e04974.

Baron, F., Cochet, M.-F., Ablain, W., Grosset, N., Madec, M.-N., Gonnet, F., Jan, S., and Gautier, M. (2006). Rapid and cost-effective method for micro-organism enumeration based on miniaturization of the conventional plate-counting technique. *Le Lait*, 86(3):251-257.

Belcour, A., Frioux, C., Aite, M., Bretaudieu, A., Hildebrand, F., and Siegel, A. (2020). Metage2metabo, microbiota-scale metabolic complementarity for the identification of key species. *eLife*.

Blasche, S., Kim, Y., Mars, R. A. T., Machado, D., Maansson, M., Kafkia, E., Milanese, A., Zeller, G., Teusink, B., Nielsen, J., Benes, V., Neves,
R., Sauer, U., and Patil, K. R. (2021). Metabolic cooperation and spa-
tiotemporal niche partitioning in a kefir microbial community. *Nature
Microbiology*, 6(2):196–208.

Borghei, R., Mobasheri, M., and Sobati, T. (2021). Genome scale constraint-
based metabolic reconstruction of propionibacterium freudenreichii DSM
20271. Preprint from Research square, pages 1–29.

Cao, W., Aubert, J., Maillard, M.-B., Boissel, F., Leduc, A., Thomas, J.-L., Deutsch, S.-M., Camier, B., Kerjouh, A., Parayre, S., Harel-Oger, M., Garric, G., Thierry, A., and Falentin, H. (2021). Fine-tuning of process parameters modulates specific metabolic bacterial activities and aroma compound production in semi-hard cheese. *Journal of Agricultural and Food Chemistry*, 69(30):8511–8529.

Carroll, N. M., Ross, R. P., Kelly, S. M., Price, N. C., Sheehan, D., and
Cogan, T. M. (1999). Characterization of recombinant acetolactate syn-
thase from Leuconostoc lactis NCW1. *Enzyme and Microbial Technology*,
25(1-2):61–67.

Caspi, R., Billington, R., Fulcher, C. A., Keseler, I. M., Kothari, A., Krum-
menacker, M., Latendresse, M., Midford, P. E., Ong, Q., Ong, W. K., Pa-
ley, S., Subhraveti, P., and Karp, P. D. (2018). The MetaCyc database of
metabolic pathways and enzymes. *Nucleic Acids Research*, 46(D1):D633–
D639.
Dank, A., van Mastrigt, O., Boeren, S., Lillevang, S. K., Abee, T., and Smid, E. J. (2021). Propionibacterium freudenreichii thrives in microaerobic conditions by complete oxidation of lactate to CO2. *Environmental Microbiology*, 23(6):3116–3129.

de Roissart, H. (DL 1994). *Bactéries lactiques : aspects fondamentaux et technologiques* / H. de Roissart, F. M. Luquet, coordonnateurs. In *Bactéries lactiques : aspects fondamentaux et technologiques*. Lorica, Uriage.

Deborde, C. and Boyaval, P. (2000). Interactions between pyruvate and lactate metabolism in Propionibacterium freudenreichii subsp. shermanii: In vivo 13C nuclear magnetic resonance studies. *Applied and Environmental Microbiology*, 66(5):2012–2020.

Diener, C., Gibbons, S. M., and Resendis-Antonio, O. (2020a). MICOM: Metagenome-scale modeling to infer metabolic interactions in the gut microbiota. *mSystems*, 5(1).

Diener, C., Gibbons, S. M., and Resendis-Antonio, O. (2020b). MICOM: Metagenome-scale modeling to infer metabolic interactions in the gut microbiota. *mSystems*, 5(1).

Ebrahim, A., Lerman, J. A., Palsson, B. O., and Hyduke, D. R. (2013). COBRApy: COnstraints-Based reconstruction and analysis for python. *BMC Systems Biology*, 7(1):74.
Falentin, H., Deutsch, S.-M., Jan, G., Loux, V., Thierry, A., Parayre, S., Maillard, M.-B., Dherbécourt, J., Cousin, F. J., Jardin, J., Siguier, P., Couloux, A., Barbe, V., Vacherie, B., Wincker, P., Gibrat, J.-F., Gail-lardin, C., and Lortal, S. (2010). The complete genome of propionibacterium freudenreichii CIRM-BIA1T, a hardy actinobacterium with food and probiotic applications. *PLoS ONE*, 5(7):e11748.

Fang, X., Lloyd, C. J., and Palsson, B. O. (2020). Reconstructing organisms in silico: Genome-scale models and their emerging applications. *Nature Reviews Microbiology*, pages 1–13.

Flahaut, N. A. L., Wiersma, A., van de Bunt, B., Martens, D. E., Schaap, P. J., Sijtsma, L., dos Santos, V. A. M., and de Vos, W. M. (2013). Genome-scale metabolic model for Lactococcus lactis MG1363 and its application to the analysis of flavor formation. *Applied Microbiology and Biotechnology*, 97(19):8729–8739.

Galimberti, A., Bruno, A., Agostinetto, G., Casiraghi, M., Guzzetti, L., and Labra, M. (2021). Fermented food products in the era of globalization: Tradition meets biotechnology innovations. *Current Opinion in Biotechnology*, 70:36–41.

Gomez, J. A., Höffner, K., and Barton, P. I. (2014). DFBAlab: A fast and reliable MATLAB code for dynamic flux balance analysis. *BMC bioinformatics*, 15(1):1–10.
Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., del Río, J. F., Wiebe, M., Peterson, P., Gérard-Marchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., and Oliphant, T. E. (2020). Array programming with NumPy. *Nature*, 585(7825):357–362.

Hunter, J. D. (2007). Matplotlib: A 2D graphics environment. *Computing in Science & Engineering*, 9(3):90–95.

King, Z. A., Dräger, A., Ebrahim, A., Sonnenschein, N., Lewis, N. E., and Palsson, B. O. (2015). Escher: A web application for building, sharing, and embedding data-rich visualizations of biological pathways. *PLoS Computational Biology*.

Kleerebezem, M., Boekhorst, J., Van Kranenburg, R., Molenaar, D., Kuipers, O. P., Leer, R., Tarchini, R., Peters, S. A., Sandbrink, H. M., Fiers, M. W., Stiekema, W., Klein Lankhorst, R. M., Bron, P. A., Hoffer, S. M., Nierop Groot, M. N., Kerkhoven, R., De Vries, M., Ursing, B., De Vos, W. M., and Siezen, R. J. (2003). Complete genome sequence of Lactobacillus plantarum WCFS1. *Proceedings of the National Academy of Sciences of the United States of America*, 100(4):1990–1995.

Langmead, B., Trapnell, C., Pop, M., and Salzberg, S. (2009). Langmead B, trapnell C, pop M, salzberg SL. Ultrafast and memory-efficient alignment
of short DNA sequences to the human genome. Genome biol 10: R25.
Genome biology, 10:R25.

Lieven, C., Beber, M. E., Olivier, B. G., Bergmann, F. T., Ataman, M., Babaei, P., Bartell, J. A., Blank, L. M., Chauhan, S., Correia, K., Diener, C., Dräger, A., Ebert, B. E., Edirisinghe, J. N., Faria, J. P., Feist, A. M., Fengos, G., Fleming, R. M. T., García-Jiménez, B., Hatzimanikatis, V., van Helvoirt, W., Henry, C. S., Hermjakob, H., Herrgård, M. J., Kaafarani, A., Kim, H. U., King, Z., Klamt, S., Klipp, E., Koehorst, J. J., König, M., Lakshmanan, M., Lee, D.-Y., Lee, S. Y., Lee, S., Lewis, N. E., Liu, F., Ma, H., Machado, D., Mahadevan, R., Maia, P., Mardinoglu, A., Medlock, G. L., Monk, J. M., Nielsen, J., Nielsen, L. K., Nogales, J., Nookaew, I., Palsson, B. O., Papin, J. A., Patil, K. R., Poolman, M., Price, N. D., Resendis-Antonio, O., Richelle, A., Rocha, I., Sánchez, B. J., Schaap, P. J., Sheriff, R. S. M., Shoaie, S., Sonnenschein, N., Teusink, B., Vilaça, P., Vik, J. O., Wodke, J. A. H., Xavier, J. C., Yuan, Q., Zakhartsev, M., and Zhang, C. (2020). MEMOTE for standardized genome-scale metabolic model testing. *Nature Biotechnology*, 38(3):272-276.

Loux, V., Mariadassou, M., Almeida, S., Chiapello, H., Hammani, A., Buratti, J., Gendrault, A., Barbe, V., Aury, J. M., Deutsch, S. M., Parayre, S., Madec, M. N., Chuat, V., Jan, G., Peterlongo, P., Azevedo, V., Le Loir, Y., and Falentin, H. (2015). Mutations and genomic islands can explain the strain dependency of sugar utilization in 21 strains of Propionibacterium
freudenreichii. BMC Genomics, 16(1):1–19.

Machado, D., Andrejev, S., Tramontano, M., and Patil, K. R. (2018). Fast automated reconstruction of genome-scale metabolic models for microbial species and communities. bioRxiv : the preprint server for biology, page 223198.

Mahadevan, R., Edwards, J. S., and Doyle, F. J. (2002). Dynamic flux balance analysis of diauxic growth in escherichia coli. Biophysical Journal, 83(3):1331–1340.

Makhlouf, A. (2006). Méthodologie Pour l'optimisation Dynamique Multicritère d'un Procédé Discontinu Alimenté: Application à La Production Bactérienne d'arômes Laitiers. PhD thesis, Institut National Polytechnique de Lorraine.

Malik, A. C., Reinbold, G. W., and Vedamuthu, E. R. (1968). An evaluation of the taxonomy of Propionibacterium. Canadian Journal of Microbiology, 14(11):1185–1191.

Mannaa, M., Han, G., Seo, Y.-S., and Park, I. (2021). Evolution of food fermentation processes and the use of multi-omics in deciphering the roles of the microbiota. Foods (Basel, Switzerland), 10(11):2861.

McKinney, W. et al. (2010). Data structures for statistical computing in python. In Proceedings of the 9th Python in Science Conference, volume 445, pages 51–56. Austin, TX.
Miyoshi, A., Rochat, T., and Gratadoux, J.-j. (2003). Miyoshi. *Genetics and Molecular Research*, 2(4):348-359.

Noronha, A., Modamio, J., Jarosz, Y., Guerard, E., Sompairac, N., Preciat, G., Daníelsdóttir, A. D., Krecke, M., Merten, D., Haraldsdóttir, H. S., Heinken, A., Heirendt, L., Magnúsdóttir, S., Ravcheev, D. A., Sahoo, S., Gawron, P., Friscioni, L., Garcia, B., Prendergast, M., Puente, A., Rodriguez, M., Roy, A., Rouquaya, M., Wiltgen, L., Žagare, A., John, E., Krueger, M., Kuperstein, I., Zinovyev, A., Schneider, R., Fleming, R. M. T., and Thiele, I. (2018). The Virtual Metabolic Human database: Integrating human and gut microbiome metabolism with nutrition and disease. *Nucleic Acids Research*, 47(D1):D614-D624.

Orth, J. D., Thiele, I., and Palsson, B. Ø. (2010). What is flux balance analysis? *Nature Biotechnology*, 28(3):245-248.

Özcan, E., Seven, M., Şirin, B., Çakır, T., Nikerel, E., Teusink, B., and Öner, E. T. (2020). Dynamic co-culture metabolic models reveal the fermentation dynamics, metabolic capacities and interplays of cheese starter cultures. *Biotechnology and Bioengineering*, 118(1):223-237.

Palles, T., Beresford, T., Condon, S., and Cogan, T. M. (1998). Citrate metabolism in Lactobacillus casei and Lactobacillus plantarum. *Journal of Applied Microbiology*, 85(1):147-154.

Palsson, B. O. and Varma, A. (1994). Metabolic flux balancing: Basic
concepts, scientific and practical use. *Biotechnology* (Reading, Mass.),
12(October):994–998.

Pogačić, T., Maillard, M.-B., Leclerc, A., Hervé, C., Chuat, V., Yee, A. L., Valence, F., and Thierry, A. (2015). A methodological approach to screen diverse cheese-related bacteria for their ability to produce aroma compounds. *Food Microbiology*, 46:145–153.

Putri, G. H., Anders, S., Pyl, P. T., Pimanda, J. E., and Zanini, F. (2022). Analysing high-throughput sequencing data in Python with HTSeq 2.0. *Bioinformatics* (Oxford, England), 38(10):2943–2945.

Quatravaux, S., Remize, F., Bryckaert, E., Colavizza, D., and Guzzo, J.
(2006). Examination of Lactobacillus plantarum lactate metabolism side
effects in relation to the modulation of aeration parameters. *Journal of
Applied Microbiology*, 101(4):903–912.

R Core Team (2021). *R: A Language and Environment for Statistical Computing*. Vienna, Austria.

Rau, M. H. and Zeidan, A. A. (2018). Constraint-based modeling in microbial food biotechnology. *Biochemical Society Transactions*, 46(2):249–260.

Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2009). edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. *Bioinformatics* (Oxford, England), 26(1):139–140.
Römer, M., Eichner, J., Dräger, A., Wrzodek, C., Wrzodek, F., and Zell, A. (2016). ZBIT bioinformatics toolbox: A web-platform for systems biology and expression data analysis. *PLoS ONE*.

Smid, E. J. and Kleerebezem, M. (2014). Production of aroma compounds in lactic fermentations. *Annual Review of Food Science and Technology*, 5(1):313–326.

Somerville, V., Grigaitis, P., Battjes, J., Moro, F., and Teusink, B. (2021). Use and limitations of genome-scale metabolic models in food microbiology. *Current Opinion in Food Science*.

Swindell, S. R., Benson, K. H., Griffin, H. G., Renault, P., Ehrlich, S. D., and Gasson, M. J. (1996). Genetic manipulation of the pathway for diacetyl metabolism in Lactococcus lactis. *Applied and Environmental Microbiology*, 62(7):2641–2643.

Tamang, J. P., Shin, D.-H., Jung, S.-J., and Chae, S.-W. (2016a). Functional properties of microorganisms in fermented foods. *Frontiers in Microbiology*, 7:578.

Tamang, J. P., Watanabe, K., and Holzapfel, W. H. (2016b). Review: Diversity of microorganisms in global fermented foods and beverages. *Frontiers in Microbiology*, 7(MAR).

Teusink, B., Wiersma, A., Molenaar, D., Francke, C., de Vos, W. M., Siezen, R. J., and Smid, E. J. (2006). Analysis of growth of lactobacillus plantarum
WCFS1 on a complex medium using a genome-scale metabolic model. *Journal of Biological Chemistry*, 281(52):40041-40048.

Thierry, A., Deutsch, S. M., Falentin, H., Dalmasso, M., Cousin, F. J., and Jan, G. (2011). New insights into physiology and metabolism of Propionibacterium freudenreichii. *International Journal of Food Microbiology*, 149(1):19-27.

Turgay, M., Bachmann, H.-P., Irmler, S., von Ah, U., Frö Hlich-Wyder, M.-T., Falentin, H., Deutsch, S.-M., Jan, G., and Thierry, A. (2020). Propionibacterium spp. and Acidipropionibacterium spp. In *Reference Module in Food Science*, Reference Module in Food Science. Elsevier.

Vallenet, D., Calteau, A., Dubois, M., Amours, P., Bazin, A., Beuvin, M., Burlot, L., Bussell, X., Fouteau, S., Gautreau, G., Lajus, A., Langlois, J., Planel, R., Roche, D., Rollin, J., Rouy, Z., Sabatet, V., and Médigue, C. (2019). MicroScope: An integrated platform for the annotation and exploration of microbial gene functions through genomic, pangenomic and metabolic comparative analysis. *Nucleic Acids Research*.

Van Rooijen, R. J., Van Schalkwijk, S., and De Vos, W. M. (1991). Molecular cloning, characterization, and nucleotide sequence of the tagatose 6-phosphate pathway gene cluster of the lactose operon of Lactococcus lactis. *Journal of Biological Chemistry*, 266(11):7176-7181.

Vaughan L. Crow (1986). Metabolism of aspartate by propionibacterium
freudenreichii subsp. shermanii: Effect on lactate fermentation. *Applied and Environmental Microbiology*, 52(2):359-365.

Wes McKinney (2010). Data Structures for Statistical Computing in Python. In van der Walt, S. and Jarrod Millman, editors, *Proceedings of the 9th Python in Science Conference*, pages 56-61.

Wick, R. R., Judd, L. M., Gorrie, C. L., and Holt, K. E. (2017a). Completing bacterial genome assemblies with multiplex MinION sequencing. *Microbial Genomics*, 3(10):1-7.

Wick, R. R., Judd, L. M., Gorrie, C. L., and Holt, K. E. (2017b). Unicycler: Resolving bacterial genome assemblies from short and long sequencing reads. *PLoS Computational Biology*, 13(6):1-22.

Widyastuti, Y., R., and Febrisiantosa, A. (2014). The role of lactic acid bacteria in milk fermentation. *Food and Nutrition Sciences*, 05(04):435-442.

Zelezniak, A., Andrejev, S., Ponomarova, O., Mende, D. R., Bork, P., and Patil, K. R. (2015). Metabolic dependencies drive species co-occurrence in diverse microbial communities. *Proceedings of the National Academy of Sciences*, 112(20):6449-6454.
A Supplementary figures

<figure><img src="image_8.png" /><figcaption>Figure S1: Model goodness of fit. All the model outputs ($\hat{y}$) are plotted versus their respective data value (y) in the training (i.e. the mono-culture data, left panel) and in the testing set (i.e. the co-culture data, right panel). Each point is colored according to the data type and the bisector is plotted (red dotted line). A global coefficient of determination is computed with the formula $R^2 = 1 - \frac{\sum_i (y_i - \hat{y}_i)^2}{\sum_i (y_i - \bar{y})^2}$, where $\bar{y}$ is the average value of y. We can observe a very good agreement of the model in the training set, and a globally good agreement in the testing set, with 3 outliers on lactose and citrate. We added time labels for these outliers to indicate to which time point they correspond.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure S2: Gene expression with metatranscriptomic data.** This heatmap displays the decile of the RPKM of gene transcript counts at a given time point (columns) and in a given micro-organism (the colored areas indicate to which bacteria the gene belongs). Then, the decile were color coded from highest deciles (dark red, highest expression in this micro-organism at this time point) to lowest decile (light red, low expression). The displayed genes were manually selected from the list of possible interactions raised by SMETANA: for a selected interaction metabolite, genes involved in production (in the donor), and in consumption (in the receiver) were retained for analysis. Selected interaction metabolites were *H2S*, ribose and phenylalanine. Additionally, genes involved in lactose consumption were added to this list. Abbreviations: Cm, molding (t = 19.5 hours); Cdm, demolding (t = 40 hours); C0w, start of ripening (t = 60 hours); C4w, Fourth week ripening (t = 732 hours); C7w, seventh week ripening (t = 1236 hours).</figcaption></figure>
<figure><img src="image_10.png" /></figure>

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

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

Figure S3: **Metabolic pathway switches during cheese production.** For each bacterium we ran a dFBA simulation and reaction fluxes per each pathways of interest in Figure 1 were retrieved at the inoculum of LAB (Linoc, t=0 hours), the inoculum of *P. freudenreichii* (Finoc, t=18 hours the molding (Cm, t=19.5 hours), demolding (Cdm, t=40hours), start of ripening (C0w, t=60 hours), fourth week ripening (C4w, t=732 hours) and seventh week ripening state (C7w, t=1236 hours). We normalized fluxes by the inoculum flux value and the average flux for each pathways is therefore represented.
# B Supplementary Material

## B.1 Pure culture data

### B.1.1 Strain pure culture data

<table><caption>Table S1: Pure culture data used for the dynamic calibration of each strain</caption><thead><tr><th rowspan="2">GEM</th><th rowspan="2">Time</th><th colspan="2">Bacterial density (g L<sup>-1</sup>)</th><th colspan="2">pH</th></tr><tr><th>rep1</th><th>rep2</th><th>rep1</th><th>rep2</th></tr></thead><tbody><tr><td rowspan="5">P. freudenreichii</td><td>0</td><td>1.452e-06</td><td>9.9e-07</td><td></td><td></td></tr><tr><td>23</td><td>7.59e-05</td><td>0.0001221</td><td></td><td></td></tr><tr><td>40</td><td>0.000561</td><td>0.000561</td><td></td><td></td></tr><tr><td>48</td><td>0.000396</td><td>0.000429</td><td></td><td></td></tr><tr><td>122</td><td>0.000759</td><td>0.00066</td><td></td><td></td></tr><tr><td rowspan="7">L. plantarum</td><td>0</td><td>9.9e-07</td><td>1.98e-06</td><td>6.7</td><td>6.7</td></tr><tr><td>5</td><td>2.97e-06</td><td>4.62e-06</td><td>6.68</td><td>6.68</td></tr><tr><td>7</td><td>5.61e-06</td><td>4.95e-06</td><td>6.67</td><td>6.65</td></tr><tr><td>9</td><td>8.58e-06</td><td>8.25e-06</td><td>6.57</td><td>6.54</td></tr><tr><td>14</td><td>1.452e-05</td><td>1.32e-05</td><td>6.55</td><td>6.53</td></tr><tr><td>16</td><td>2.772e-05</td><td>1.947e-05</td><td>6.52</td><td>6.5</td></tr><tr><td>79</td><td>3.63e-05</td><td>3.861e-05</td><td>5.7</td><td>5.71</td></tr><tr><td rowspan="7">L. lactis</td><td>0</td><td>3.3e-07</td><td>1.32e-06</td><td>6.7</td><td>6.7</td></tr><tr><td>5</td><td>5.445e-05</td><td>5.016e-05</td><td>6.49</td><td>6.49</td></tr><tr><td>7</td><td>5.841e-05</td><td>7.821e-05</td><td>6.33</td><td>6.3</td></tr><tr><td>9</td><td>6.138e-05</td><td>7.161e-05</td><td>6.07</td><td>6.08</td></tr><tr><td>14</td><td>5.709e-05</td><td>7.227e-05</td><td>5.93</td><td>5.94</td></tr><tr><td>16</td><td>5.115e-05</td><td>4.356e-05</td><td>5.88</td><td>5.85</td></tr><tr><td>79</td><td>8.2005e-05</td><td>4.125e-05</td><td>5.12</td><td>5.11</td></tr></tbody></table>
B.1.2 Acids dosage for *P. freudenreichii*

<table><caption>Table S2: Acids dosage data in mM for *P. freudenreichii*.</caption><thead><tr><th>Compounds</th><th>Concentration T<sub>0</sub></th><th>Concentration T<sub>f</sub> = 89 h</th></tr></thead><tbody><tr><td>Lactate</td><td>0.183171</td><td>0.087478</td></tr><tr><td>Acetate</td><td>0</td><td>0.051122</td></tr><tr><td>Succinate</td><td>0</td><td>0.003142</td></tr><tr><td>Propionate</td><td>0</td><td>0.112176</td></tr></tbody></table>

B.2 Manual refinement of GEMs

To ensure the production of extracellular compounds, we added manual modifications on intracellular bounds. These modifications are recapitulated in Table S3.
<table><caption>Table S3: Manual refinement of the GEMs and corresponding justifications.</caption><thead><tr><th>GEM</th><th>Reaction</th><th>Bound</th><th>Value</th><th>Justification</th></tr></thead><tbody><tr><td rowspan="4">L. lactis</td><td>ACTD2</td><td>lower<br>upper</td><td>0<br>0</td><td>Allows the activation of the acetolactate pathway</td></tr><tr><td>ACLDC</td><td>lower<br>upper</td><td>0<br>2</td><td>Allows the activation of the acetolactate pathway</td></tr><tr><td>ACLS</td><td>lower<br>upper</td><td>2<br>1000</td><td>Allows the activation of the acetolactate pathway</td></tr><tr><td>LACZ</td><td>lower<br>upper</td><td>0.0006<br>2</td><td>Force the consumption of lactose</td></tr><tr><td rowspan="5">L. plantarum</td><td>PFL</td><td>lower<br>upper</td><td>0<br>0</td><td>Permit the production of both isomers of lactate from lactose.</td></tr><tr><td>LCTSt3ipp</td><td>lower<br>upper</td><td>0<br>0</td><td>Leads the flux of lactose into the Glycolysis pathway</td></tr><tr><td>ACTD2</td><td>lower<br>upper</td><td>0<br>0</td><td>Allows the activation of the acetolactate pathway</td></tr><tr><td>GALUi</td><td>lower<br>upper</td><td>-30<br>10</td><td>Permit the production of both isomers of lactate from lactose.</td></tr><tr><td>ACLS</td><td>lower<br>upper</td><td>2<br>1000</td><td>Allows the activation of the acetolactate pathway</td></tr><tr><td rowspan="8">P. freudenreichii</td><td>6PGALSZ</td><td>lower<br>upper</td><td>0<br>0</td><td>Tagatose pathway was used</td></tr><tr><td>XYLI2</td><td>lower<br>upper</td><td>0<br>0</td><td>Production of fructose from glucose instead of glucose-6-phosphate. Glycolysis was not activated</td></tr><tr><td>POX2</td><td>lower<br>upper</td><td>0<br>0</td><td>Regulate flux of acetate</td></tr><tr><td>PPAKr</td><td>lower<br>upper</td><td>0<br>0</td><td>Inhibit acetate over production</td></tr><tr><td>LACZ</td><td>lower<br>upper</td><td>0<br>5</td><td>Allow consumption of lactose</td></tr><tr><td>PTAr</td><td>lower<br>upper</td><td>-1000<br>8</td><td>Allow the appropriate ratio of acetate/propionate</td></tr><tr><td>2131pyrpp</td><td>lower<br>upper</td><td>8<br>1000</td><td>Allow the appropriate ratio of acetate/propionate</td></tr></tbody></table>
## B.3 Community experimental data

<table><caption>Table S4: Co-culture data used for testing community model predictions: numbering and pH</caption><thead><tr><th rowspan="2">GEM</th><th rowspan="2">steps</th><th rowspan="2">Time (h)</th><th colspan="4">Bacterial population (g L<sup>-1</sup>)</th></tr><tr><th>Rep1</th><th>Rep2</th><th>Rep3</th><th>Rep4</th></tr></thead><tbody><tr><td rowspan="5">P. freudenreichii</td><td>Cm</td><td>19.5</td><td>3.16E-06</td><td>1.91E-06</td><td>4.11E-06</td><td>5.46E-06</td></tr><tr><td>Cdm</td><td>40.0</td><td>8.65E-05</td><td>3.61E-05</td><td>3.32E-05</td><td>3.89E-05</td></tr><tr><td>C0w</td><td>60.0</td><td>7.18E-05</td><td>2.20E-05</td><td>4.01E-05</td><td>4.46E-05</td></tr><tr><td>C4w</td><td>732.0</td><td>2.83E-04</td><td>5.68E-05</td><td>1.34E-04</td><td>1.07E-04</td></tr><tr><td>C7w</td><td>1236.0</td><td>2.58E-04</td><td>5.23E-05</td><td>1.29E-04</td><td>1.27E-04</td></tr><tr><td rowspan="5">L. plantarum</td><td>Linoc</td><td>0.0</td><td>1.69E-07</td><td>3.04E-08</td><td>3.23E-08</td><td>4.95E-08</td></tr><tr><td>Cm</td><td>19.5</td><td>6.71E-07</td><td>5.94E-07</td><td>8.75E-07</td><td>1.17E-06</td></tr><tr><td>Cdm</td><td>40.0</td><td>1.59E-05</td><td>6.20E-06</td><td>1.05E-05</td><td>1.17E-05</td></tr><tr><td>C0w</td><td>60.0</td><td>6.08E-05</td><td>3.00E-05</td><td>2.99E-05</td><td>2.57E-05</td></tr><tr><td>C4w</td><td>732.0</td><td>9.79E-05</td><td>6.55E-05</td><td>1.02E-04</td><td>9.44E-05</td></tr><tr><td rowspan="5">L. lactis</td><td>Linoc</td><td>0.0</td><td>6.29E-08</td><td>3.30E-07</td><td>2.08E-07</td><td>2.31E-07</td></tr><tr><td>Cm</td><td>19.5</td><td>6.51E-05</td><td>1.43E-04</td><td>1.78E-04</td><td>2.21E-04</td></tr><tr><td>Cdm</td><td>40.0</td><td>5.27E-04</td><td>5.27E-04</td><td>3.96E-04</td><td>4.74E-04</td></tr><tr><td>C0w</td><td>60.0</td><td>1.27E-03</td><td>4.64E-04</td><td>3.81E-04</td><td>3.76E-04</td></tr><tr><td>C4w</td><td>732.0</td><td>4.80E-04</td><td>2.85E-04</td><td>2.95E-04</td><td>2.87E-04</td></tr><tr><td rowspan="5">pH</td><td>Linoc</td><td>0.0</td><td>6.72E+00</td><td>6.71E+00</td><td>6.72E+00</td><td>6.72E+00</td></tr><tr><td>Cm</td><td>19.5</td><td>6.40E+00</td><td>6.34E+00</td><td>6.41E+00</td><td>6.30E+00</td></tr><tr><td>Cdm</td><td>40.0</td><td>5.42E+00</td><td>5.48E+00</td><td>5.54E+00</td><td>5.37E+00</td></tr><tr><td>C0w</td><td>60.0</td><td>5.35E+00</td><td>5.44E+00</td><td>5.35E+00</td><td>5.30E+00</td></tr><tr><td>C4w</td><td>732.0</td><td>5.29E+00</td><td>5.14E+00</td><td>5.08E+00</td><td>5.15E+00</td></tr><tr><td></td><td>C7w</td><td>1236.0</td><td>5.26E+00</td><td>5.21E+00</td><td>5.13E+00</td><td>5.16E+00</td></tr></tbody></table>
<table><caption>Table S5: Co-culture data used for testing community model predictions: metabolites</caption><thead><tr><th rowspan="2">Compound</th><th rowspan="2">steps</th><th rowspan="2">Time (h)</th><th colspan="2">Concentration (mM)</th></tr><tr><th>Rep1</th><th>Rep2</th></tr></thead><tbody><tr><td rowspan="6">Lactose</td><td>Linoc</td><td>0.0</td><td>1.42E-01</td><td>1.41E-01</td></tr><tr><td>Cm</td><td>19.5</td><td>6.72E-02</td><td>6.24E-02</td></tr><tr><td>Cdm</td><td>40.0</td><td>2.44E-02</td><td>2.40E-02</td></tr><tr><td>C0w</td><td>60.0</td><td>2.05E-02</td><td>1.60E-02</td></tr><tr><td>C4w</td><td>732.0</td><td>3.15E-03</td><td>3.50E-03</td></tr><tr><td>C7w</td><td>1236.0</td><td>2.43E-03</td><td>2.19E-03</td></tr><tr><td rowspan="5">Lactate</td><td>Cm</td><td>19.5</td><td>4.74E-02</td><td>4.49E-02</td></tr><tr><td>Cdm</td><td>40.0</td><td>1.18E-01</td><td>1.12E-01</td></tr><tr><td>C0w</td><td>60.0</td><td>1.28E-01</td><td>1.31E-01</td></tr><tr><td>C4w</td><td>732.0</td><td>1.43E-01</td><td>1.46E-01</td></tr><tr><td>C7w</td><td>1236.0</td><td>1.12E-01</td><td>1.18E-01</td></tr><tr><td rowspan="3">propionate</td><td>C0w</td><td>60.0</td><td>1.62E-03</td><td>1.46E-03</td></tr><tr><td>C4w</td><td>732.0</td><td>4.46E-02</td><td>4.00E-02</td></tr><tr><td>C7w</td><td>1236.0</td><td>6.70E-02</td><td>6.74E-02</td></tr><tr><td rowspan="6">acetate</td><td>Linoc</td><td>0.0</td><td>2.30E-04</td><td>2.06E-04</td></tr><tr><td>Cm</td><td>19.5</td><td>7.30E-03</td><td>6.02E-03</td></tr><tr><td>Cdm</td><td>40.0</td><td>1.53E-02</td><td>1.34E-02</td></tr><tr><td>C0w</td><td>60.0</td><td>1.75E-02</td><td>1.62E-02</td></tr><tr><td>C4w</td><td>732.0</td><td>3.34E-02</td><td>2.95E-02</td></tr><tr><td>C7w</td><td>1236.0</td><td>4.31E-02</td><td>4.20E-02</td></tr><tr><td rowspan="6">succinate</td><td>Linoc</td><td>0.0</td><td>0.00E+00</td><td>0.00E+00</td></tr><tr><td>Cm</td><td>19.5</td><td>7.88E-04</td><td>6.65E-04</td></tr><tr><td>Cdm</td><td>40.0</td><td>1.62E-03</td><td>1.74E-03</td></tr><tr><td>C0w</td><td>60.0</td><td>1.18E-03</td><td>1.65E-03</td></tr><tr><td>C4w</td><td>732.0</td><td>2.66E-03</td><td>2.22E-03</td></tr><tr><td>C7w</td><td>1236.0</td><td>2.60E-03</td><td>3.37E-03</td></tr><tr><td rowspan="6">citrate</td><td>Linoc</td><td>0.0</td><td>1.18E-02</td><td>1.14E-02</td></tr><tr><td>Cm</td><td>19.5</td><td>6.54E-03</td><td>6.74E-03</td></tr><tr><td>Cdm</td><td>40.0</td><td>1.15E-03</td><td>1.45E-03</td></tr><tr><td>C0w</td><td>60.0</td><td>7.65E-04</td><td>8.04E-04</td></tr><tr><td>C4w</td><td>732.0</td><td>0.00E+00</td><td>0.00E+00</td></tr><tr><td>C7w</td><td>1236.0</td><td>0.00E+00</td><td>0.00E+00</td></tr></tbody></table>
## B.4 Exchanged metabolites highlighted by computational approaches

### B.4.1 MICOM

<table><caption>Table S6: Exchangeable metabolites highlighted by MICOM on the cheese bacterial community. Abbreviations: Ll, *L. lactis*; Lp, *L. plantarum*, Pf, *P. freudenreichii*</caption><thead><tr><th>Bigg ID</th><th>Metacyc ID</th><th>Name</th><th>Metacyc Ontology</th><th>Export fluxes</th><th>Import fluxes</th></tr></thead><tbody><tr><td>acald_e</td><td>ACETALD</td><td>acetaldehyde</td><td>Aldehydes-Or-Ketones, Aldehydes</td><td>Pf; Ll</td><td>Lp</td></tr><tr><td>akg_e</td><td>2-KETOGLUTARATE</td><td>2-oxoglutarate</td><td>Acids, Organic-Acids</td><td>Pf</td><td>Lp</td></tr><tr><td>ala_D_e</td><td>D-ALANINE</td><td>D-ALANINE</td><td>All-Amino-Acids, Acids, Amino-Acids, Organic-Acids</td><td>Pf; Lp</td><td>Ll</td></tr><tr><td>ala_L_e</td><td>L-ALPHA-ALANINE</td><td>L-ALPHA-ALANINE</td><td>All-Amino-Acids, Acids, Amino-Acids, Organic-Acids</td><td>Pf; Lp</td><td>Ll</td></tr><tr><td>co2_e</td><td>CARBON-DIOXIDE</td><td>CARBON-DIOXIDE</td><td>Others, Others</td><td>Pf; Ll</td><td>Lp</td></tr><tr><td>coa_e</td><td>CO-A</td><td>CO enzyme A</td><td>Groups, Others</td><td>Pf; Ll</td><td>Lp</td></tr><tr><td>fe2_e</td><td>FE+2</td><td>Fe2+</td><td>Ions, Inorganic-Ions, Cations</td><td>Lp</td><td>Pf; Ll</td></tr><tr><td>gly_e</td><td>GLY</td><td>Glycine</td><td>All-Amino-Acids, Acids, Amino-Acids, Organic-Acids</td><td>Lp</td><td>Pf; Ll</td></tr><tr><td>glyc1t_e</td><td>GLYCOLATE</td><td>GLYCOLATE</td><td>Acids, Organic-Acids</td><td>Pf</td><td>Lp</td></tr><tr><td>gthox_e</td><td>OXIDIZED-GLUTATHIONE</td><td>glutathione disulfide</td><td>ORGANOSULFUR, All-Glutathiones</td><td>Lp</td><td>Ll</td></tr><tr><td>gthrd_e</td><td>GLUTATHIONE</td><td>glutathione</td><td>ORGANOSULFUR, All-Glutathiones, Thiols</td><td>Ll</td><td>Lp</td></tr><tr><td>h2o2_e</td><td>HYDROGEN-PEROXIDE</td><td>hydrogen peroxide</td><td>Peroxides, Others</td><td>Ll</td><td>Lp</td></tr><tr><td>h2o_e</td><td>WATER</td><td>H2O</td><td>Pseudo-Compounds, Others</td><td>Pf; Ll</td><td>Lp</td></tr><tr><td>h2s_e</td><td>HS</td><td>hydrogen sulfide</td><td>Ions, Anions, Inorganic-Ions</td><td>Ll</td><td>Lp; Pf</td></tr><tr><td>ile_L_e</td><td>ILE</td><td>isoleucine</td><td>All-Amino-Acids, Acids, Amino-Acids, Organic-Acids</td><td>Lp</td><td>Pf; Ll</td></tr><tr><td>lac_D_e</td><td>D-LACTATE</td><td>d-lactate</td><td>Acids, Organic-Acids</td><td>Pf</td><td>Ll; Lp</td></tr><tr><td>lac_L_e</td><td>L-LACTATE</td><td>l-lactate</td><td>Acids, Organic-Acids</td><td>Ll</td><td>Lp; Pf</td></tr><tr><td>mal_L_e</td><td>MAL</td><td>malate</td><td>Acids, Organic-Acids</td><td>Pf; Lp</td><td>Ll</td></tr><tr><td>nh4_e</td><td>AMMONIUM</td><td>ammonium</td><td>Ions, Inorganic-Ions, Cations</td><td>Lp</td><td>Pf; Ll</td></tr><tr><td>phe_L_e</td><td>PHE</td><td>phenylalanine</td><td>Aromatics, All-Amino-Acids, Acids, Amino-Acids, Organic-Acids, Organic-aromatic-compounds</td><td>Pf</td><td>Lp</td></tr><tr><td>pro_L_e</td><td>PRO</td><td>proline</td><td>All-Amino-Acids, Acids, Amino-Acids, Organic-Acids</td><td>Ll</td><td>Pf</td></tr><tr><td>ser_D_e</td><td>D-SERINE</td><td>serine</td><td>All-Amino-Acids, Acids, Amino-Acids, Organic-Acids</td><td>Ll</td><td>Lp; Pf</td></tr><tr><td>ser_L_e</td><td>SER</td><td>serine</td><td>All-Amino-Acids, Acids, Amino-Acids, Organic-Acids</td><td>Ll</td><td>Lp; Pf</td></tr><tr><td>succ_e</td><td>SUC</td><td>succinate</td><td>Acids, Organic-Acids</td><td>Ll</td><td>Pf</td></tr><tr><td>xan_e</td><td>XANTHINE</td><td>xanthine</td><td>Others, Others</td><td>Ll</td><td>Lp; Pf</td></tr></tbody></table>
## B.4.2 SMETANA

<table><caption>Table S7: Exchangeable metabolites highlighted by SMETANA on the cheese bacterial community. Abbreviations: Ll, *L. lactis*; Lp, *L. plantarum*, Pf, *P. freudenreichii*</caption><thead><tr><th>Bigg ID</th><th>Metacyc ID</th><th>Name</th><th>Metacyc Ontology</th><th>Export fluxes</th><th>Import fluxes</th></tr></thead><tbody><tr><td rowspan="4">acald_e</td><td rowspan="4">ACETALD</td><td rowspan="4">acetaldéhyde</td><td rowspan="4">Aldehydes-Or-Ketones, Aldehydes</td><td>Pf</td><td>Ll</td></tr><tr><td>Pf</td><td>Lp</td></tr><tr><td>Lp</td><td>Ll</td></tr><tr><td>Ll</td><td>Lp</td></tr><tr><td rowspan="2">akg_e</td><td rowspan="2">2-KETOGLUTARATE</td><td rowspan="2">2-oxoglutarate</td><td rowspan="2">Acids, Organic-Acids</td><td>Lp</td><td>Pf</td></tr><tr><td>Pf</td><td>Lp</td></tr><tr><td>fum_e</td><td>FUM</td><td>fumarate</td><td>Acids, Organic-Acids</td><td>Pf</td><td>Lp</td></tr><tr><td rowspan="3">glyc_e</td><td rowspan="3">GLYCEROL</td><td rowspan="3">GLYCEROL</td><td rowspan="3">Alcohols, All-Carbohydrates, Sugar-alcohols, Carbohydrates</td><td>Lp</td><td>Pf</td></tr><tr><td>Ll</td><td>Pf</td></tr><tr><td>Ll</td><td>Lp</td></tr><tr><td>glyc1t_e</td><td>GLYCOLATE</td><td>GLYCOLATE</td><td>Acids, Organic-Acids</td><td>Lp</td><td>Pf</td></tr><tr><td rowspan="3">h2s_e</td><td rowspan="3">HS</td><td rowspan="3">hydrogen sulfide</td><td rowspan="3">Ions, Anions, Inorganic-Ions</td><td>Lp</td><td>Pf</td></tr><tr><td>Pf</td><td>Lp</td></tr><tr><td>Ll</td><td>Lp</td></tr><tr><td rowspan="2">lac__D_e</td><td rowspan="2">D-LACTATE</td><td rowspan="2">d-lactate</td><td rowspan="2">Acids, Organic-Acids</td><td>Pf</td><td>Lp</td></tr><tr><td>Ll</td><td>Lp</td></tr><tr><td rowspan="2">lac__L_e</td><td rowspan="2">L-LACTATE</td><td rowspan="2">l-lactate</td><td rowspan="2">Acids, Organic-Acids</td><td>Ll</td><td>Lp</td></tr><tr><td>Ll</td><td>Lp</td></tr><tr><td rowspan="4">mal__L_e</td><td rowspan="4">MAL</td><td rowspan="4">malate</td><td rowspan="4">Acids, Organic-Acids</td><td>Pf</td><td>Ll</td></tr><tr><td>Lp</td><td>Ll</td></tr><tr><td>Pf</td><td>Lp</td></tr><tr><td>Ll</td><td>Lp</td></tr><tr><td>phe__L_e</td><td>PHE</td><td>phenylalanine</td><td>Aromatics, All-Amino-Acids, Acids, Amino-Acids, Organic-Acids, Organic-aromatic-compounds</td><td>Lp</td><td>Pf</td></tr><tr><td rowspan="4">rib__D_e</td><td rowspan="4">D-Ribopyranose</td><td rowspan="4">D-Ribopyranose</td><td rowspan="4">Aldehydes-Or-Ketones, All-Carbohydrates, Aldehydes, Carbohydrates</td><td>Ll</td><td>Pf</td></tr><tr><td>Lp</td><td>Pf</td></tr><tr><td>Pf</td><td>Lp</td></tr><tr><td>Ll</td><td>Lp</td></tr><tr><td rowspan="6">ser__D_e</td><td rowspan="6">D-SERINE</td><td rowspan="6">serine</td><td rowspan="6">All-Amino-Acids, Acids, Amino-Acids, Organic-Acids</td><td>Ll</td><td>Pf</td></tr><tr><td>Lp</td><td>Pf</td></tr><tr><td>Pf</td><td>Ll</td></tr><tr><td>Lp</td><td>Ll</td></tr><tr><td>Lp</td><td>Ll</td></tr><tr><td>Ll</td><td>Lp</td></tr><tr><td rowspan="2">succ_e</td><td rowspan="2">SUC</td><td rowspan="2">succinate</td><td rowspan="2">Acids, Organic-Acids</td><td>Ll</td><td>Pf</td></tr><tr><td>Lp</td><td>Pf</td></tr><tr><td rowspan="2">xan_e</td><td rowspan="2">XANTHINE</td><td rowspan="2">xanthine</td><td rowspan="2">Others, Others</td><td>Ll</td><td>Pf</td></tr><tr><td>Lp</td><td>Pf</td></tr></tbody></table>
### B.4.3 Metage2Metabo

<table><caption>Table S8: Metabolites predicted to be producible through metabolic cross-feeding only by Metage2Metabo</caption><thead><tr><th>Bigg ID</th><th>Metacyc ID</th><th>Name</th><th>Metacyc Ontology</th></tr></thead><tbody><tr><td>galt1p</td><td>D-galactose-1-phosphate</td><td>a D-galactopyranose 1-phosphate</td><td>All-Carbohydrates, Carbohydrates</td></tr><tr><td>benzcoa</td><td>BENZOYLCOA</td><td>benzoyl-CoA</td><td>Esters, Thioesters</td></tr><tr><td>xtsn</td><td>XANTHOSINE</td><td>XANTHOSINE</td><td>Organic-heterocyclic-compound, All-Nucleosides, Nitrogen-Molecular-Entities, Organic-heteromonocyclic-compounds, Organic-Heteropolycyclic-Compounds, Organonitrogen-Compounds, Nucleosides</td></tr><tr><td>glyald</td><td>GLYCERALD</td><td>glyceraldehyde</td><td>Aldehydes-Or-Ketones, All-Carbohydrates, Aldehydes, Carbohydrates</td></tr><tr><td>3hocoa</td><td>ø</td><td>ø</td><td>ø</td></tr><tr><td>3hhdcoa</td><td>CPD0-2232</td><td>(S)-3-hydroxyhexadecanoyl-CoA</td><td>Thioesters, Esters</td></tr><tr><td>3hdcoa</td><td>CPD0-2244</td><td>(S)-3-hydroxydecanoyl-CoA</td><td>Thioesters, Esters</td></tr><tr><td>3odcoa</td><td>CPD0-2123</td><td>3-oxodecanoyl-CoA</td><td>Thioesters, Esters</td></tr><tr><td>3hhcoa</td><td>OH-HEXANOYL-COA</td><td>(S)-3-hydroxyhexanoyl-CoA</td><td>Thioesters, Esters</td></tr><tr><td>udcpp</td><td>UNDECAPRENYL-P</td><td>all-trans-undecaprenyl phosphate</td><td>Lipids, Polyisoprenoids</td></tr><tr><td>3oocoa</td><td>CPD0-2106</td><td>3-oxooctanoyl-CoA</td><td>Thioesters, Esters</td></tr></tbody></table>

## B.5 Dynamic bounds on exchange reactions

In this section, we recapitulate and justify the regulation functions modulating the metabolite exchanges. For each GEM, the exchange reactions are listed, and their lower and upper bounds are specified.

### *L. lactis*
• EX\_ucts\_e:

$$ \text{- lower: } \max\left(-\frac{m_{ucts_e}}{\Delta t * \sum_{i \in \mathcal{M}(ucts_e)} b_i}, -\mu_{max,ucts}^{lact} * 10^{(-k_{lac} * \phi_{undiss})} - \mu_{min,ucts}^{lact}\right) \quad (11) $$

- upper: 1000

Justification: When lactose is depleted, the first term is activated, modeling uniform sharing of available lactose among the set of bacteria $M(lac_L)$ metabolizing lactose. Otherwise, a regulation by the pH is added: the amount of undissociated lactic acid molecules is computed from lactate concentrations through the function $\phi_{undiss}$ (see (8) for the expression of this function); then, an exponential increase from the value $-\mu_{max,ucts}^{lact} - \mu_{min,ucts}^{lact}$ to the value $-\mu_{min,ucts}^{lact}$ when the undissociated lactic acid concentration increases, with exponential rate $k_{lac}$. The values of $\mu_{max,ucts}^{lact}$, $\mu_{min,ucts}^{lact}$ and $k_{lac}$ are inferred from co-culture experiments.

• EX\_lac\_L\_e, EX\_lac\_D\_e:

- lower: 0

$$ \text{- upper: } \min\left(\frac{m_{ucts_e}}{\Delta t * \sum_{i \in \mathcal{M}(ucts_e)} b_i}, \mu_{max,ucts} * 10^{(-k_{lac} * \phi_{undiss})} + \mu_{min,ucts} * 4\right) \quad (12) $$

Justification: As lactate production is linked to lactose availability, a dynamic modulation for lactate export mimics lactose import (see eq. (13)), with a stoichio-metric factor 4.0.

• EX\_ac\_e:

- lower: 0

$$ \text{- upper: } \max\left(-\frac{m_{ucts_e}}{\Delta t * \sum_{i \in \mathcal{M}(ucts_e)} b_i}, v_{i,ac_e}^{exp}\right) $$

Justification: Acetate production is regulated by lactose availability when lactose is depleted, and is otherwise exported according to an intrinsic physiological export limit $v_{i,ac_e}^{exp} = 0.5$
• EX\_diact\_e, EX\_btd\_RR\_e, EX\_cit\_e:

$$max(-\frac{m_j}{\Delta t * \sum_{i \in \mathcal{M}(j)} b_i}, v_{i,j}^{int})$$

- lower: for i = L. lactis and j ∈ {EX\_diact\_e, EX\_btd\_RR\_e, EX\_cit\_e, EX\_ac\_e}

- upper: 1000

Justification: Usual consumption dynamic limitation (see eq. (6)). When substrate is not limited, import is constrained by intrinsic physiological limit modeled by $v_{i,j}^{int} = -8$ for j = L. lactis and j ∈ {EX\_diact\_e, EX\_btd\_RR\_e} and $v_{i,j}^{int} = -1$ for j = EX\_cit\_e, and otherwise by nutrient availability. Available substrate is uniformly shared by consuming micro-organisms.

Additionnaly, the lower bounds of EX\_glc\_D\_e (glucose), EX\_coa\_e (coenzyme A) and EX\_starch1200\_e (starch) were set to 0, since those compounds do not appear in milk composition.

## L. plantarum

• EX\_lcts\_e:

$$\text{- lower: } max(-\frac{m_{lcts_e}}{\Delta t * \sum_{i \in \mathcal{M}(lcts_e)} b_i}, -\mu_{max,lcts}^{plant} * 10^{(-k_{lac} * \phi_{undiss})} - \mu_{min,lcts}^{plant}) \quad (13)$$

- upper: 1000

Justification: When lactose is depleted, the first term is activated, modeling uniform sharing of available lactose among the set of bacteria $\mathcal{M}(lac_L)$ metabolizing lactose. Otherwise, a regulation by the pH is added: the amount of undissociated lactic acid molecules is computed from lactate concentrations through the function $\phi_{undiss}$ (see (8) for the expression of this function); then, an exponential increase from the value $-\mu_{max,lcts}^{plant} - \mu_{min,lcts}^{plant}$ to the value $-\mu_{min,lcts}^{plant}$ when the undissociated lactic acid concentration increases, with exponential rate $k_{lac}$. The values of $\mu_{max,lcts}^{plant}, \mu_{min,lcts}^{plant}$ and $k_{lac}$ are inferred from co-culture experiments.
• EX\_lac\_L\_e, EX\_lac\_D\_e:

- lower: 0

$$ \text{- upper: } \min\left(\frac{m_{lcts_e}}{\Delta t * \sum_{i \in \mathcal{M}(lcts_e)} b_i}, \mu_{max,lcts} * 10^{(-k_{lac} * \phi_{undiss})} + \mu_{min,lcts}\right) * 4 \quad (14) $$

Justification: As lactate production is linked to lactose availability, a dynamic modulation for lactate export mimics lactose import (see eq. (13)), with a stoichiometric factor 4.0.

• EX\_ac\_e:

- lower: 0

$$ \text{- upper: } \max\left(-\frac{m_{lcts_e}}{\Delta t * \sum_{i \in \mathcal{M}(lcts_e)} b_i}, v_{i,ac_e}^{exp}\right) $$

Justification: Acetate production is regulated by lactose availability when lactose is depleted, and is otherwise exported according to an intrinsic physiological export limit $v_{i,ac_e}^{exp} = 1.0$

• EX\_btd\_RR\_e:

$$ \max\left(-\frac{m_j}{\Delta t * \sum_{i \in \mathcal{M}(j)} b_i}, v_{i,j}^{int}\right), $$

for $j = L. plantarum$ and $j \in \{EX\_btd\_RR\_e\}$

- lower:

- upper: 1000

Justification: Usual consumption dynamic limitation (see eq. (6)). When substrate is not limited, import is constrained by intrinsic physiological limits modeled by $v_{i,j}^{int} = -5$ for $i = L. plantarum$ and $j \in \{EX\_ac\_e, EX\_btd\_RR\_e\}$, and otherwise by nutrient availability. Available substrate is uniformly shared by consuming micro-organisms.

Additionally, the lower bounds of EX\_glc\_D\_e (glucose), EX\_gal\_e (galactose), EX\_coa\_e (coenzyme A) and EX\_dha\_e (dihydroxyacetone) were set to 0, since those compounds do not appear in milk composition.
P. freudenreichii

• EXlleLcE:

- lower: -12

- upper: 1000

Justification: Static consumption of serine

• EX_cite:

- lower: -12

- upper: 1000

Justification: Static consumption of citrate

• EX_gulleLcE:

- lower: -12

- upper: 1000

Justification: Static consumption of glutamate

• EX_alleLcE:

- lower: -12

- upper: 1000

Justification: Static consumption of alanine
* EX\_lac\_L\_e, EX\_lac\_D\_e:

$$ \text{- lower:} \qquad \max\left(-\frac{m_{lac_{L,e}}}{\Delta t * \sum_{i \in \mathcal{M}(lac_L)} b_i}, F_{lac} * K_{freud} * \frac{b_{freud}}{B_{p, freud}}\right) \qquad (15) $$

- upper: 1000

Justification: When lactate is depleted, the first term is activated, modeling uniform sharing of available lactate among the set of bacteria $M(lac_L)$ metabolizing L\_lactate (resp. D\_lactate). Otherwise, a flux $F_{lac}$, estimated from metabolomic data in pure culture (see section B.7), is applied after modulation by the fraction $\frac{b_{freud}}{B_{p, freud}}$, where $B_{p, freud}$ is the plateau biomass concentration value in pure culture and $b_{freud}$ is the current biomass density, and the factor $K_{freud}$, which is inferred from the pure culture growth data. This modulation models a different order of magnitude of metabolism yield in pure or co-culture.

* EX\_lcts\_e:

$$ \text{- lower:} \qquad \max\left(-\frac{m_{lcts_e}}{\Delta t * \sum_{i \in \mathcal{M}(lcts_e)} b_i}, F_{lcts}\right) $$

- upper: 1000

Justification: When lactose is depleted, the first term is activated, modeling uniform sharing of available lactose among the set of bacteria $M(lac_L)$ metabolizing lactose. Otherwise, a constant flux $F_{lcts}$, estimated from metabolomic data in pure culture (see section B.7), is applied.

* EX\_ac\_e, EX\_ppa\_e, EX\_succ\_e:

- lower: 0

$$ \text{- upper:} \qquad F_i * K_{freud} * \frac{b_{freud}}{B_{p, freud}}, \quad \text{for } i \in \{\text{EX\_ac\_e, EX\_ppa\_e, EX\_succ\_e}\} $$

Justification: The maximal production rate is defined by $F_i$, for $i \in \{\text{EX\_ac\_e, EX\_ppa\_e, EX\_succ\_e}\}$, which is estimated from metabolomic data in pure culture (see section B.7). The same modulation term than in 15 is applied, modeling a different order of magnitude of metabolism yield in pure or co-culture.

Additionnaly, the lower bounds of EX\_starch1200\_e (starch) and EX\_dha\_e (di-
hydroxyacetone) were set to 0, since those compounds do not appear in milk composition.

## B.6 Dynamic bounds on intracellular reactions

When regulating metabolites involved in dynamic bounds for exchange reactions are de- completed, the corresponding substrate import vanishes, inducing low metabolic fluxes across the whole GEM. These low metabolic fluxes can conflict with the fixed bounds on in- tracellular reactions introduced in Table S3, leading to infeasible FBA models. To avoid such infeasibility, we added additional dynamic bounds on the intracellular reactions. We stress that these dynamic bounds on intracellular reactions have a very different status than those on the exchange reactions: while the dynamic bounds on the exchange reactions drive the dynamics of the metabolic flux in the GEM according to substrate availability all along the simulation, the dynamics bounds on the intracellular reactions are only active when a substrate is depleted. They are "security" bounds that activate in case of vanishing substrate and weak metabolic activity, to avoid numerical infeasibility.

For each reaction contained in Table S3 indicating manually refined reactions, we add the following dynamic regulation

$$lower\_bound = min(\frac{m}{\Delta t \sum_{i \in M(m)} b_i}, manual\ lower\ bound)$$

Hence, as the dynamic regulation is always positive, this regulation is active for strictly positive values of *manual lower bound*, the value indicated in Table S3. In that case, this regulation becomes active only when the metabolite *m* is depleted.

## B.7 Estimating bounds on metabolite production for *P. freudenreichii*

The dFBA equation in the pure culture experiment for *i* = *P. freudenreichii* reads,
<table><caption>Table S9: list of optimized parameters obtained after optimization and used in the fitted models.</caption><thead><tr><th>parameter</th><th><i>L. lactis</i></th><th><i>L. plantarum</i></th><th><i>P. freudenreichii</i></th></tr></thead><tbody><tr><td>&lambda;</td><td>&empty;</td><td>0.351</td><td>1.496</td></tr><tr><td><i>k</i><sub>lac</sub></td><td>255.591</td><td>300</td><td>&empty;</td></tr><tr><td>oxygene &mu;<sub>max,lactose</sub></td><td>30.666</td><td>7.968</td><td>&empty;</td></tr></tbody></table>

$$ \partial_t b_i = \mathcal{R}_i(b_i) \mu_{i,i} \left( (c_{min,i}^{ex}, c_{max,i}^{ex})(b,m) \right) b_i \quad (16) $$

$$ \partial_t m_j = \mu_{i,j} \left( (c_{min,i}^{ex}, c_{max,i}^{ex})(b,m) \right) b_i. \quad (17) $$

Assuming a constant production flux $\mu_{i,j}$ and integrating eq. (17) in time, we get

$$ \mu_{i,j} = \alpha \frac{m_j(T) - m_j(0)}{\int_0^T b_i(t) dt} $$

where $\alpha$ is a coefficient accounting for approximation errors.

We then evaluated $\mu_{i,j}$ from initial and final metabolite concentrations, and numerical integration (trapeze scheme) of the biomass, taking $\alpha = 1.08$. The resulting value is taken as lower bound for consumption and upper bound for production.

## B.8 Parameter values

The optimized parameter values are recapitulated in Table S9. The remaining parameter values are recapitulated in Table S10.
<table><caption>Table S10: Parameter values used for simulations</caption><thead><tr><th>parameter</th><th>description</th><th>value [unit]</th></tr></thead><tbody><tr><td>c<sub>1</sub></td><td>pH linear dependence on lactate</td><td>-11.1937 [-]</td></tr><tr><td>c<sub>2</sub></td><td>intercept of pH dependence on lactate</td><td>2.8346 [-]</td></tr><tr><td>dt<sub>com</sub></td><td>integration time step in co-culture</td><td>0.5 [h]</td></tr><tr><td>Tf<sub>com</sub></td><td>Final experiment time in co-culture</td><td>1236 [h]</td></tr><tr><td>dt<sub>mono</sub></td><td>integration time step in co-culture</td><td>0.2 [h]</td></tr><tr><td>Tf<sub>mono</sub></td><td>Final experiment time in co-culture</td><td>122,80, 80[h]</td></tr><tr><td>citrate<sub>concentration</sub></td><td>Inoculum citrate</td><td>2.23 [g/kg]</td></tr><tr><td>acetate<sub>concentration</sub></td><td>Inoculum acetate</td><td>0.01 [g/kg]</td></tr><tr><td>diacetyle<sub>concentration</sub></td><td>Inoculum diacetyle</td><td>0.0 [g/kg]</td></tr><tr><td>succinate<sub>concentration</sub></td><td>Inoculum succinate</td><td>0.085750 [g/kg]</td></tr><tr><td>lactose<sub>concentration</sub></td><td>Inoculum lactose</td><td>48 [g/kg]</td></tr><tr><td>lactate<sub>concentration</sub></td><td>Inoculum lactate</td><td>0.0 [g/kg]</td></tr><tr><td>thresh<sub>qs, cocult</sub></td><td>Carrying capacity for P. freudenreichii, L. lactis and L. plantarum in co-culture</td><td>-3.85, -3.47, -4.04 [log<sub>10</sub> g kg<sup>-1</sup>]</td></tr><tr><td>thresh<sub>qs, monocult</sub></td><td>Carrying capacity for P. freudenreichii, L. lactis and L. plantarum in mono-culture</td><td>-3.15, -4.15, -4.50 [log<sub>10</sub> g kg<sup>-1</sup>]</td></tr><tr><td>inoculum_com</td><td>Bacteria inoculation for P. freudenreichii, L. lactis, L. plantarum</td><td>6.1, 5.7, 5.2 [log<sub>10</sub> CFU/g]</td></tr><tr><td>oxygene</td><td>Value of the oxygene for P. freudenreichii, L. lactis, L. plantarum</td><td>2, 7, 11 [mmol.gDW<sup>-1</sup>.hr<sup>-1</sup>]</td></tr><tr><td>intrinsuc<sub>flux</sub></td><td>Value of the intrinsic flux for P. freudenreichii, L. lactis, L. plantarum</td><td>-10, -8, -5 [mmol.gDW<sup>-1</sup>.hr<sup>-1</sup>]</td></tr></tbody></table>